% Problem 8.10
disp('Problem 8.10')
clear
K=4;
n=[K];
d=[1 2 0];

% Transfer function of the system
disp('(a) Uncompensated system:')
printsys(n,d,'s')

% Open-loop poles
disp('Open-loop poles:')
p=roots(d)

subplot(221)
rlocus(n,d)
axis([-2 0.1 -3 3])
grid
title('(a) Uncompensated system')

disp('Closed-loop uncompensated system:')
[ncl,dcl]=cloop(n,d,-1)
printsys(ncl,dcl,'s')

% Operating point sd
OS=0.15; ts=2;
xi=sqrt(log(OS)^2/(log(OS)^2+pi^2)); xi=0.52;
om_n=3/(ts*xi)
sd=-xi*om_n + i*om_n*sqrt(1-xi^2)
th_p1=angle(sd-p(1))
th_p2=angle(sd-p(2))
alpha_c=pi+th_p1+th_p2;
while alpha_c > 2*pi,
    alpha_c=alpha_c-2*pi;
end
alpha_c

% Zero of PD compensator
zc=om_n*(xi*tan(alpha_c)+sqrt(1-xi^2))/tan(alpha_c)

nc=[1 zc]; 

disp('System with a phase-lead compensator:')
n_c=conv(nc,n);
d_c=d;

% Static gain of compensated system
dp1=abs(sd-p(1))
dp2=abs(sd-p(2))
dzc=abs(sd+zc)
disp('Static gain:')
Kc=dp1*dp2/(K*dzc)
subplot(222)
rlocus(n_c,d_c)
axis([-20 0.1 -10 10])
grid
title('(b) Compensated system')

disp('Closed-loop compensated system:')
[n_ccl,d_ccl]=cloop(Kc*n_c,d_c,-1);
printsys(n_ccl,d_ccl,'s')

disp('Poles of closed-loop compensated system:')
roots(d_ccl)

subplot(212)
t=0:0.1:5;
[y,x]=step(ncl,dcl,t);
hold
[yc,xc]=step(n_ccl,d_ccl,t);
plot(t,y,'--',t,yc,'-')
grid
xlabel('time [s]')
ylabel('y,yc')