%  Figure 6.62      Feedback Control of Dynamic Systems, 8e
%                   Franklin, Powell, Emami
%   

clear all
close all
clf
 
w = logspace(-2,1,500);

num = 3;

den = conv([2 1],[1 1]);
den = conv(den,[0.5 1]);

[mag,phas] = bode(num,den,w);
[OLGM,OLPM,OLWcg,OLWcp] = margin(mag,phas,w);
data0 = [OLGM,OLPM,OLWcg,OLWcp]


%Lag compensator 

numl = 3*[5 1];
denl = [15 1];

numc = conv(num,numl);
denc = conv(den,denl);

[magc,phasc] = bode(numc,denc,w);
[D1GM,D1PM,D1Wcg,D1Wcp] = margin(magc,phasc,w);
data1 = [D1GM,D1PM,D1Wcg,D1Wcp]

[data0; data1]

dencl = denc+[0 0 0 numc];

t = 0:.1:20;
y = step(numc,dencl,t);

figure(1)

subplot(2,1,1)

loglog(w,mag,'-',w,ones(500,1),'-',w,ones(500,1),'-','LineWidth',2);

text( 2*0.01, 4*1, "K = 3")

axis([.01 10 .1 10])
grid on;
%xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Bode Plot for lag-compensation design for Ex 6.18 (a) magnitude');

subplot(2,1,2)

semilogx(w,phas,'-',w,-180*ones(500,1),w,-180*ones(500,1),'LineWidth',2);

axis([.01 10 -250 -50])
grid on;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase');


figure(2)

subplot(2,1,1)

loglog(w,ones(500,1),'-',w,magc,'--',w,ones(500,1),'-','LineWidth',2);

text( 2*0.1, 7*1, "D_c = 3*[5 1]/[15 1]")

axis([.01 10 .1 10])
grid on;
%xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Bode Plot for lag-compensation design for Ex 6.18 (a) magnitude');

subplot(2,1,2)

semilogx(w,-180*ones(500,1),w,phasc,'--',w,-180*ones(500,1),'LineWidth',2);

text( 2*0.1, 7*1, "D_c = 3*[5 1]/[15 1]")

axis([.01 10 -250 -50])
grid on;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase');


figure(3)

subplot(2,1,1)

loglog(w,mag,'-',w,magc,'--',w,ones(500,1),'-','LineWidth',2);

text( 2*0.01, 4*1, "K = 3")
text( 2*0.1, 7*1, "D_c = 3*[5 1]/[15 1]")

axis([.01 10 .1 10])
grid on;
%xlabel('\omega (rad/sec)');
ylabel('Magnitude');
title('Bode Plot for lag-compensation design for Ex 6.18 (a) magnitude');

subplot(2,1,2)

semilogx(w,phas,'-',w,phasc,'--',w,-180*ones(500,1),'LineWidth',2);

axis([.01 10 -250 -50])
grid on;
xlabel('\omega (rad/sec)');
ylabel('Phase (deg)');
title('(b) phase');

