%  Example 10.6 (continuation).  Generate the step response and the 
%  Nyquist plot.
%  Define the aerodynamics model.
%    States: alpha, theta_dot
%    Inputs: delta
%    Outputs:  eta, d*theta_2dot, theta_dot
Za=-5.24;Zd=-0.4;Ma=-47;Mtd=-4.69;Md=-1134;V=1.82;d=8.55e-4;
Aan=[Za  1
     Ma Mtd];
Ban=[Zd
     Md];
Can=[V*Za   0
     Ma*d Mtd*d
     0     1  ];
Dan=[V*Zd
     Md*d
       0 ];
Aeron=pck(Aan,Ban,Can,Dan);
%  Define the servo model.
wn=189;z=0.7;
As=[   0      1
    -wn*wn -2*z*wn];
Bs=[  0
    wn*wn];
Cs=eye(2);
Ds=[0 0]';
Servo=pck(As,Bs,Cs,Ds);
%  Perform the system interconnection for the step response.
systemnames = 'Aeron Servo K';
inputvar = '[etac]';
outputvar = '[Aeron(1);Servo(1);Servo(2)]';
input_to_Aeron='[Servo(1)]';
input_to_Servo='[K]';
input_to_K='[etac-Aeron(1)-Aeron(2); Aeron(3)]';
sysoutname='Gcl';
cleanupsysic='yes';
sysic
%  Generate the closed loop state model.
[Acl,Bcl,Ccl,Dcl]=unpck(Gcl);
%  Generate the step response.
t=0:.0015:1;
y=step(Acl,Bcl,Ccl,Dcl,1,t);
%  Plot the step response.
set(0,'DefaultAxesFontName','times')
set(0,'DefaultAxesFontSize',16)
set(0,'DefaultTextFontName','times')
figure(1)
clf
subplot(311),plot(t,30*y(:,1))
ylabel('Acceleration (g)')
grid
subplot(312),plot(t,30*57.3*y(:,2))
ylabel('Fin angle (deg)')
grid
subplot(313),plot(t,30*57.3*y(:,3))
xlabel('Time (sec)')
ylabel('Fin rate (deg/sec)')
grid
%  Perform the system interconnection for the loop transfer function.
systemnames = 'Aeron Servo K';
inputvar = '[deltacin]';
outputvar = '[-K]';
input_to_Aeron='[Servo(1)]';
input_to_Servo='[deltacin]';
input_to_K='[-Aeron(1)-Aeron(2); Aeron(3)]';
sysoutname='Gloop';
cleanupsysic='yes';
sysic
%  Generate the state model of the loop transfer function.
[Al,Bl,Cl,Dl]=unpck(Gloop);
%  Generate the nyquist plot.
figure(2)
clf
[re,im]=nyquist(Al,Bl,Cl,Dl);
plot(re,im,'r-',re,-im,'r-',-1,0,'r+',re(62),im(62),...
   'r^',re(62),-im(62),'rv')
xlabel('Real axis')
ylabel('Imaginary axis')
axis('square')
axis([-2 2 -2 2])
grid
