% Problem 8.12
disp('Problem 8.12')
clear
n=[4];
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([-7 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.1; ts=2;
xi=sqrt(log(OS)^2/(log(OS)^2+pi^2)); xi=xi+0.005;
om_n=3/(xi*ts)
sd=-xi*om_n + i*om_n*sqrt(1-xi^2)

% Phase contribution
th_p1=angle(sd-p(1))
th_p2=angle(sd-p(2))
th_c=pi+th_p1+th_p2;
while th_c - 2*pi > 0,
    th_c=th_c-2*pi;
end
disp('Controller phase contribution [rad]:')
th_c

% Zero of the Phase-lead compensator
zc=3.5
phi=angle(sd+zc)
j=1;
pc=xi*om_n+om_n*sqrt(1-xi^2)*tan(th_c-phi+pi/2)
nc=[1 zc];
dc=[1 pc];
disp('Phase-lead compensator:')
printsys(nc,dc,'s')

disp('System with a phase-lead compensator:')
n_c=conv(nc,n);
d_c=conv(dc,d);
subplot(222)
rlocus(n_c,d_c)
axis([-7 1 -3 3])
title('(b) Compensated system')
grid

% Static gain of compensated system
dp1=abs(sd-p(1))
dp2=abs(sd-p(2))
dpc=abs(sd+pc)
dzc=abs(sd+zc)
disp('Static gain:')
Kc=dp1*dp2*dpc/(4*dzc)

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 the closed-loop compensated system:')
roots(d_ccl)

subplot(212)
t=0:0.1:10;
[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')