% Problem 9.9
clear
%  System model
num=[0 0 0 40];
den=conv([1 4 0],[1 10]);

% (a,b) Bode diagram and stability margins
w=logspace(-1,1,200);
[mag,phase]=bode(num,den,w);

% Gain and phase margins
margin(mag,phase,w)
subplot(211)
grid
subplot(212)
grid

% (c) Gm=10dB
[Gm,Pm,Wcp,Wcg]=margin(mag,phase,w);
fprintf('\n(b) Phase margin (for K=1)    = %f degrees\n\n',Pm)
fprintf('    Gain margin (for K=1)     = %f dB\n\n',20*log10(Gm))
fprintf('    Phase crossover frequency = %f rad/s\n\n',Wcp)
fprintf('    Gain crossover frequency  = %f rad/s\n\n',Wcg)
K=10^((-10+20*log10(Gm))/20);
Essp=0;
Kv=40*K/(4*10); Essv=1/Kv;
Ka=0; Essa=1/Ka;
fprintf('\n(c) Static gain K (for Gm=10dB)  = %f \n',K)
fprintf('    Unit step steady state error = 0\n')
fprintf('    Unit ramp steady state error = %f \n',Essv)
fprintf('    Unit par. steady state error = %f \n',Essa)

% (d) Pm=45 degrees
j=1;
while phase(j) > -135,
   j=j+1;
end
K=1/mag(j);
fprintf('(d) Static gain K (at Pm=45 degrees) = %f \n',K)
fprintf('    new crossover frequency wcp = %f \n',w(j))

% Closed-loop system poles
[nc,dc]=cloop(K*num,den,-1);
roots(dc)