disp(' ')
disp(' ')
disp('*** This code is written for MATLAB 5.2 ***')
disp(' ')
echo on
%
% Lateral-Directional Airplane Stability and Control Analysis
%
% Ref: Airplane Flight Dynamics and Automatic Flight Controls
%      Part 1, 1994
%      by Jan Roskam
%
% See Data input format of Appendix B.
% Sample data for Cessna 182, cruise configuration, pages 480-482.

echo off
constant=BoilEX;
disp('Boiler Express, cruise configuration')

% To change to a different aircraft, only the code above needs to be changed.
%
% Reference Geometry
S=constant(20)		    % Wing area, ft*ft
cbar=constant(21)		% wing mean geometric chord, ft
span=constant(22)		% Wing span, ft
%
% Steady State (Trim) Flight Condition
h=constant(59);		% Altitude, ft
U1=constant(58)		% True airspeed, ft/sec
rho=0.00237691267925741*(1-6.87558563248308e-06*h)^(4.25591641274834) %slug/ft3, troposphere
qbar=.5*rho*U1*U1     %Dynamic pressure, lbf/ft^2
Xcg=constant(57)		% X-location of center of gravity, fraction of cbar
alpha1=constant(60)	% Trim angle of attack, degree
gamma1=0   % Trim flight path angle, degree (This is usually so. Check it!)<<<

%
% Mass data
W=constant(1)		% Weight, pounds (lbf)
Ixxb=constant(4)	% X-direction moment of inertia, slug*ft*ft
Iyyb=constant(5)	% Y-direction moment of inertia, slug*ft*ft
Izzb=constant(6)	% Z-direction moment of inertia, slug*ft*ft
Ixzb=constant(7)	% XZ-direction product of inertia, slug*ft*ft
%
% Lateral-Directional Coefficients and Stability Derivatives
Clb=constant(40)
Clp=constant(43)
Clr=constant(44)
Cyb=constant(34)
Cyp=constant(37)
Cyr=constant(38)
Cnb=constant(51)
Cntb=0            %(This is usually so. Check it!) <<<<<
Cnp=constant(54)
Cnr=constant(55)
%
% Lateral-Directional Control Derivatives
Clda=constant(41)
Cldr=constant(42)
Cyda=constant(35)
Cydr=constant(36)
Cnda=constant(52)
Cndr=constant(53)
% Miscellaneous inputs
g=32.17;		% Acceleration of gravity, ft/sec*sec
d2r=pi/180;
r2d=180/pi;
%
% Computation of Dimensional Stability and Control Derivatives
%
% Preliminary calculations
%
% Stability axis inertia data, see page 346, eqn 5.94
% Rotate axes through an angle alpha1 (see fig 5.15 page 347)
alpha1rad=alpha1*d2r;
csq=cos(alpha1rad)^2;
ssq=sin(alpha1rad)^2;
c2a=cos(2.*alpha1rad);
s2a=sin(2.*alpha1rad);
Ixxs=Ixxb*csq+Izzb*ssq-Ixzb*s2a
Izzs=Izzb*csq+Ixxb*ssq+Ixzb*s2a
Ixzs=0.5*(Ixxb-Izzb)*s2a+Ixzb*c2a

Abar=Ixzs/Ixxs;   % See page 349
Bbar=Ixzs/Izzs;

mass=W/g	% mass, slugs

% See page 348
qsdm=qbar*S/mass;
qsbd2mu=qsdm*span/(2*U1);
qsbdixx=qbar*S*span/Ixxs;
qsbbd2ixxu=qsbdixx*span/(2*U1);
qsbdizz=qbar*S*span/Izzs;
qsbbd2izzu=qsbdizz*span/(2*U1);
theta1rad=(alpha1+gamma1)*d2r;


Yb =qsdm*Cyb
Yp =qsbd2mu*Cyp
Yr =qsbd2mu*Cyr
Yda=qsdm*Cyda
Ydr=qsdm*Cydr

Lb =qsbdixx*Clb
Lp =qsbbd2ixxu*Clp
Lr =qsbbd2ixxu*Clr
Lda=qsbdixx*Clda
Ldr=qsbdixx*Cldr

Nb =qsbdizz*Cnb
Ntb=qsbdizz*Cntb
Np =qsbbd2izzu*Cnp
Nr =qsbbd2izzu*Cnr
Nda=qsbdizz*Cnda
Ndr=qsbdizz*Cndr

%
% Computation of system matrices
%	

%
% Assemble the A and B matrices. See page 349.
% x=[beta(rad) p(rad/sec) r(rad/sec) phi(rad) psi(rad)]'
% u=[aileron(rad) rudder(rad)]'
% E*xdot=a*x + b*u
a=[Yb Yp Yr-U1 g*cos(theta1rad) 0
   Lb Lp Lr       0             0
   Nb Np Nr       0             0
    0  1  0       0             0
	0  0  1       0             0];
b=[Yda Ydr
   Lda Ldr
   Nda Ndr
    0   0
	0   0];
E=[U1    0      0   0  0
   0     1    -Abar 0  0
   0   -Bbar    1   0  0
   0     0      0   1  0
   0     0      0   0  1];

Einv=inv(E);
A=Einv*a;
B=Einv*b;
% Outputs y=[beta(deg) p(deg/sec) r(deg/sec) phi(deg) psi(deg)]
C=r2d*eye(5);
D=zeros(5,2);


%
% Analyze the linear equations of motion
%
sys=ss(A,B,C,D);
set(sys,'statename',{'beta(r)' 'p(r/s)' 'r(r/s)' 'phi(r)' 'psi(r)'})
set(sys,'inputname',{'deltaA(r)' 'deltaR(r)'})
set(sys,'outputname',{'beta(d)' 'p(d/s)' 'r(d/s)' 'phi(d)' 'psi(d)'})
sys
[Wn,Z,Poles]=damp(sys)
tfsys=tf(sys)
zpksys=zpk(sys)