disp(' ')
disp(' ')
disp(' ')
disp(' Start Here')
echo on
%
% Lateral-Directional Airplane Stability and Control Analysis
%
% Ref: Airplane Flight Dyanmics and Automatic Flight Controls
%      Part 1, 1994
%      by Jan Roskam
%
% See Data input format of Appendix B
%
% Sample data for LearJet, approach configuration, pages 522-523.
%
% Aircraft and configuration
echo off
aircraft='Learjet 24';
configuration='approach';
disp([aircraft '  ' configuration])
%
% Reference Geometry
S=230		% Wing area, ft*ft
cbar=7		% wing mean geometric chord, ft
span=34		% Wing span, ft
%
% Steady State (Trim) Flight Condition
h=0;		% Altitude, ft
M=.152		% Mach number, nondimensional
U1=170		% True airspeed, ft/sec
qbar=34.3	% Dynamic pressure, lbf/ft*ft
Xcg=.32		% X-location of center of gravity, fraction of cbar
alpha1=5	% Trim angle of attack, degree
gamma1=0   % Trim flight path angle, degree

%
% Mass data
W=13000		% Weight, pounds (lbf)
Ixxb=28000	% X-direction moment of inertia, slug*ft*ft
Iyyb=18800	% Y-direction moment of inertia, slug*ft*ft
Izzb=47000	% Z-direction moment of inertia, slug*ft*ft
Ixzb=1300	% XZ-direction product of inertia, slug*ft*ft
%
% Lateral-Directional Coefficients and Stability Derivatives,
%            nondimensional
Clb=-.173
Clp=-.39
Clr=.45
Cyb=-.73
Cyp=0
Cyr=.4
Cnb=.15
Cntb=0
Cnp=-.13
%Cnr=-.26
%Cnr=-.50
%Cnr=-.40
%Cnr=-.30
%Cnr=-.20
%Cnr=-.70
%Cnr=-.90
%Cnr=-1.50
Cnr=0
%
% Lateral-Directional Control Derivatives
Clda=.149
Cldr=.014
Cyda=0
Cydr=.14
Cnda=-.05
Cndr=-.074
% 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)
