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

% Transfer function of the system
disp('(a) Uncompensated system:')
printsys(K*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(K*n,d,-1);

% Operating point sd
OS_d=0.1; ts_d=2;
xi_d=sqrt(log(OS_d)^2/(log(OS_d)^2+pi^2)); xi_d=xi_d+0.005;
om_nd=3/(xi_d*ts_d)
s_d=-xi_d*om_nd + om_nd*sqrt(1-xi_d^2)*i

% Controller phase contribution
th_p1=angle(s_d-p(1))
th_p2=angle(s_d-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_nd*(xi_d*tan(alpha_c)+sqrt(1-xi_d^2))/tan(alpha_c)
if zc>0
   nc=[1 zc]
else
   disp('System cannot be compensated by a single compensator!!!')
   disp('Introducing double PD compensator.')
   alpha_c=alpha_c/2;
   zc=om_nd*(xi_d*tan(alpha_c)+sqrt(1-xi_d^2))/tan(alpha_c)
   nc=conv([1 zc],[1 zc])
end

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

% Static gain of compensated system
dp1=abs(s_d-p(1))
dp2=abs(s_d-p(2))
dzc=abs(s_d+zc)
if length(nc)>2
   dzc=dzc*dzc;
end
disp('Static gain:')
Kc=dp1*dp2/dzc

subplot(222)
rlocus(n_c,d_c)
axis([-15 0.1 -6 6])
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')
