% You should run LatSC.m before running this script.

d2r=pi/180;
echo off
ifig=0;
t1=0:.01:5;					% time scale to expose roll mode
mag=20*d2r;
uvector=doublet(t1,2,4,mag);	% deltaA doublet of ten degree
zero=zeros(length(t1),1);
u=[uvector zero];
y=lsim(sys,u,t1);


ifig=ifig+1;
figure(ifig)
clf
subplot(3,1,1)
ts=['Aileron response, ',aircraft,', ',configuration,' configuration'];
plot(t1,y(:,2))
title(ts)
xlabel('time(sec)')
ylabel('p(deg/sec)')
text2(.1,.15,'Note Da controls roll rate through the roll mode.')
subplot(3,1,2)
plot(t1,y(:,3))
xlabel('time(sec)')
ylabel('r(deg/sec)')
subplot(3,1,3)
plot(t1,u(:,1))
xlabel('time(sec)')
ylabel('Da(rad)')
text2(.1,.15,'The roll mode involves primarily motion variables p and phi.')




ifig=ifig+1;
figure(ifig)
clf
subplot(4,1,1)
plot(t1,y(:,1))
xlabel('time(sec)')
ylabel('beta(deg)')
ts=['Aileron response, ',aircraft,', ',configuration,' configuration'];
title(ts)
subplot(4,1,2)
plot(t1,y(:,4))
xlabel('time(sec)')
ylabel('phi(deg)')
text2(.1,.15,'Note phi is the integral of roll rate.')
kk=find(y(:,4)>=30)
% time to 30 degrees is
T30=t1(kk(1))
ts2=['Time to 30 degrees= ',num2str(T30),' sec.'];
text2(.1,.25,ts2)
hold on; plot(t1(kk(1)),y(kk(1),4),'rx')

subplot(4,1,3)
plot(t1,y(:,5))
xlabel('time(sec)')
ylabel('psi(deg)')
subplot(4,1,4)
plot(t1,u(:,1))
xlabel('time(sec)')
ylabel('Da(rad)')
text2(.1,.15,'The roll mode involves primarily motion variables p and phi.')
%%%%%%%%%%%%

% Rudder doublet input
t2=0:.05:20;
mag=10*d2r;
uvector=doublet(t2,.2,.4,mag);	% deltaR doublet of ten degree
zero=zeros(length(t2),1);
u=[zero uvector];
y=lsim(sys,u,t2);

ifig=ifig+1;
figure(ifig)
clf


subplot(3,1,1)
plot(t2,y(:,2))
ts=['Rudder doublet response, ',aircraft,', ',configuration,' configuration'];
title(ts)
xlabel('time(sec)')
ylabel('p(deg/sec)')
text2(.1,.15,'Note the slightly unstable dutch roll mode.')
subplot(3,1,2)
plot(t2,y(:,3))
xlabel('time(sec)')
ylabel('r(deg/sec)')
text2(.1,.15,'The dutch rol mode involves all motion variables.')
subplot(313)
plot(t2,u(:,2))
xlabel('time(sec)')
ylabel('Dr(rad)')
text2(.1,.15,'Rudder doublets excite the dutch roll mode.')


ifig=ifig+1;
figure(ifig)
clf
subplot(4,1,1)
plot(t2,y(:,1))
xlabel('time(sec)')
ylabel('beta(deg)')
ts=['Rudder doublet response, ',aircraft,', ',configuration,' configuration'];
title(ts)
text2(.1,.15,'Note the slightly unstable dutch roll mode.')
subplot(4,1,2)
plot(t2,y(:,4))
xlabel('time(sec)')
ylabel('phi(deg)')
text2(.1,.15,'The dutch rol mode involves all motion variables.')
subplot(4,1,3)
plot(t2,y(:,5))
xlabel('time(sec)')
ylabel('psi(deg)')
subplot(414)
plot(t2,u(:,2))
xlabel('time(sec)')
ylabel('Dr(rad)')
text2(.1,.15,'Rudder doublets excite the dutch roll mode.')



%%%%%%%%%%%%%%%


% initial condition response to expose spiral mode
t3=0:.1:100;
zero=zeros(length(t3),1);
u=[zero zero];
%[V,D] = EIG(A,B)
[V,D]=eig(A);
x0=V(:,3);
y=lsim(sys,u,t3,x0);
ifig=ifig+1;
figure(ifig)
clf
subplot(5,1,1)
ts=['Initial condition response to show spiral mode, ',aircraft,', ',configuration,' configuration'];
plot(t3,y(:,2))
title(ts)
xlabel('time(sec)')
ylabel('p(deg/sec)')
text2(.1,.15,'Note the spiral mode is unstable.')

subplot(5,1,2)
plot(t3,y(:,3))
xlabel('time(sec)')
ylabel('r(deg/sec)')
subplot(5,1,3)
plot(t3,y(:,1))
xlabel('time(sec)')
ylabel('beta(deg)')
subplot(5,1,4)
plot(t3,y(:,4))
xlabel('time(sec)')
ylabel('phi(deg)')
text2(.1,.15,'Note phi is the integral of roll rate.')
text2(.1,.35,'Note the spiraling motion, sometimes called the death spiral.')
subplot(5,1,5)
plot(t3,y(:,5))
xlabel('time(sec)')
ylabel('psi(deg)')

%%%%%%%%%%%%
