% Script to analyse the aircraft system assuming flat earth.
% Version 1.1, 2/22/01. Perturbed inputs are now written to the workspace
% by the Simulink simulation and used as small pert. inputs
% in the linear simulation.

disp(' ');disp('Start here');disp(' ');
disp('Compute a quick trim')
Vt=220.1% ft/sec
Hp=5000
constant=Cessna182;
[x1,u1]=QuickTrim(Vt,Hp,constant)
flag=3; t=0;
[y1,x0,str,ts] = aircraft(t,x1,u1,flag,constant,x1);
[a,b,c,d]=linmod('FlatEarth2',x1,u1); %Find linear model of nonlinear pendulum
[Wn,Z]=damp(a) % Find properties of the poles of linearized system
poles=eig(a)
PSYS=ss(a,b,c,d); %Create a linear state space system
disp(' ')
disp('Small Perturbation Initial condition for linear simulation.')
disp('Note that this is X(0)-XR, where X(0) is the same I.C. as')
disp('used in the nonlinear simulation.')
X0=x1;
XR=x1;
YR=y1;
dx0=X0-XR; %Use this small perturbation initial condition for linear simulation
input='Time response comparison for elevator input. Flat Earth Model, Cessna 182';
tstart=taircraft(1)
tstop=taircraft(length(taircraft))
t=tstart:.02:tstop;
u=zeros(length(t),4);
%YI = INTERP1(X,Y,XI)
u(:,1)=interp1(taircraft,deltaE,t');
u(:,2)=interp1(taircraft,deltaA,t');
u(:,3)=interp1(taircraft,deltaR,t');
u(:,4)=interp1(taircraft,deltabhp,t');
[y,t]=lsim(PSYS,u,t,dx0); %Do linear simulation
%Plot the linear results and the nonlinear results
%  which have been stored in tnlsim and ynlsim.
figure(1) 
disp('Take care to add the reference state xR to the small perturbation')
disp('solution in order to get the total approximate solution.')
disp('Only the total approximate solution should be compared to the ')
disp('results from the nonlinear simulation.')

figure(1)
subplot(311)
plot(t,y(:,5)+YR(5),'-',taircraft,yaircraft(:,5),':')
xlabel('time (sec)')
ylabel('Q (rad/sec)')
legend('linear sim','nonlinear sim')
title(input)
axis([0 5 -0.05 0.05])

subplot(312)
plot(t,y(:,8)+YR(8),'-',taircraft,yaircraft(:,8),':')
xlabel('time (sec)')
ylabel('theta (rad)')
legend('linear sim','nonlinear sim')
axis([0 5 -0.05 0.05])

subplot(313)
plot(t,y(:,14)+YR(14),'-',taircraft,yaircraft(:,14),':')
xlabel('time (sec)')
ylabel('alpha (rad)')
legend('linear sim','nonlinear sim')
axis([0 5 -0.05 0.05])

figure(2)
subplot(311)
plot(t,y(:,6)+YR(6),'-',taircraft,yaircraft(:,6),':')
xlabel('time (sec)')
ylabel('R (rad/sec)')
legend('linear sim','nonlinear sim')
title(input)

subplot(312)
plot(t,y(:,9)+YR(9),'-',taircraft,yaircraft(:,9),':')
xlabel('time (sec)')
ylabel('psi (rad)')
legend('linear sim','nonlinear sim')

subplot(313)
plot(t,y(:,15)+YR(15),'-',taircraft,yaircraft(:,15),':')
xlabel('time (sec)')
ylabel('beta (rad)')
legend('linear sim','nonlinear sim')

figure(3)
subplot(311)
plot(t,y(:,4)+YR(4),'-',taircraft,yaircraft(:,4),':')
xlabel('time (sec)')
ylabel('P (rad/sec)')
legend('linear sim','nonlinear sim')
title(input)

subplot(312)
plot(t,y(:,7)+YR(7),'-',taircraft,yaircraft(:,7),':')
xlabel('time (sec)')
ylabel('phi (rad)')
legend('linear sim','nonlinear sim')

subplot(313)
plot(t,y(:,11)+YR(11),'-',taircraft,yaircraft(:,11),':')
xlabel('time (sec)')
ylabel('Y (ft)')
legend('linear sim','nonlinear sim')

poles=poles'
Z=Z'
abel=1
TranFunc=zpk(Z,poles,abel)
