disp('*** This code is written for MATLAB 5.2 ***')
echo on
%
% Longitudinal 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
aircraft='Learjet 24';
configuration='approach';
%
% 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;	% Angle of attack, degree
gamma1=0	% Steady state flight path angle, deg, see figure page 66
%					angle between horizontal and Xs axis (Wind at trim)
%
% 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
%
% Steady State (trim) Coefficients
CL1=1.64;	% Lift coefficient at trim, nondimensional
CD1=.256;	% Drag coefficient at trim, nondimensional
CTx1=.256;	% Thrust coefficient at trim, nondimensional
Cm1=0;		% Aerodynamic Pitching Moment coefficient at trim, nondimensional
CmT1=0;		% Thrust Pitching Moment coefficient at trim, nondimensional
%
%Longitudinal Coefficients and Stability Derivatives
CD0=.043;       % paracite drag, nondimensional
CDu=0;          % drag variation with speed, nondimensional
CDalpha=1.06;   % drag variation with alpha, nondimensional
CTxu=-.6;       % thrust coefficient vatiation with speed, nondimensional
CL0=1.2;        % lift coefficient at zero alpha, nondimensional
CLu=.04;        % lift coefficient variation with speed, nondimensional
CLalpha=5.04;	% Lift curve slope, nondimensional
CLalphadot=1.6; % lift coefficient variation with alphadot, nondimensional
CLq=4.1;        % lift coefficient variation with pitch rate, nondimensional
CM0=.047;       % pitching moment coef. at zero alpha, nondimensional
Cmu=-.01;       % pitch moment coef. variation with pitch rate, nondimensional
Cmalpha=-.66;	% Static longitudinal stabiity derivative, nondimensional
Cmalphadot=-5;	% Lag of downwash stability derivative, nondimensional
Cmq=-13.5;		% Swish effect in pitch, nondimensional
CmTu=.006;      % pitch moment due to thrust variation with speed, nondimensional
CmTalpha=0;     % pitch moment due to thrust variation with alpha, nondimensional
%
%Longitudinal Control Derivatives
CDdeltae=0;		% Drag from elevator, nondimensional
CLdeltae=.4;	% Lift from elevator, nondimensional
Cmdeltae=-.98;	% Pitching moment from elevator, nondimensional
CDih=0;			% Drag from incidence of the horizontal tail, nondimensional
CLih=.85;		% Lift from incidence of the horizontal tail, nondimensional
Cmih=-2.1;		% Pitlching moment from incidence of the horizontal tail, nondimensional
%
% Miscellaneous inputs
g=32.17;		% Acceleration of gravity, ft/sec*sec

%
% Computation of Dimensional Stability and Control Derivatives
%
% Preliminary calculations
% Stability axis inertia data, see page 346, eqn 5.94
Iyys=Iyyb;	% These are the same since that rotation is angle alpha1 
%				about the Y-axis
mass=W/g	% mass, slugs
theta1=gamma1+alpha1;	% Steady state 
d2r=pi/180;
r2d=180/pi;
gamma1rad=gamma1*d2r;

% See page 319
Xu=-qbar*S*(CDu+2*CD1)/(mass*U1)	% 1/sec
XTu=qbar*S*(CTxu+2*CTx1)/(mass*U1)	% 1/sec
Xalpha=-qbar*S*(CDalpha-CL1)/mass	% ft/sec*sec
Xdeltae=-qbar*S*CDdeltae/mass		% ft/sec*sec
%
Zu=-qbar*S*(CLu+2*CL1)/(mass*U1)	% 1/sec
Zalpha=-qbar*S*(CLalpha+CD1)/mass	% ft/sec*sec
Zalphadot=-qbar*S*cbar*CLalphadot/(2*mass*U1)	% ft/sec
Zq=-qbar*S*cbar*CLq/(2*mass*U1)		% ft/sec
Zdeltae=-qbar*S*CLdeltae/mass		% ft/sec*sec
%
Mu=qbar*S*cbar*(Cmu+2*Cm1)/(Iyys*U1)	% 1/ft/sec
MTu=qbar*S*cbar*(CmTu+2*CmT1)/(Iyys*U1)	% 1/ft/sec
Malpha=qbar*S*cbar*Cmalpha/Iyys			% 1/sec*sec
MTalpha=qbar*S*cbar*CmTalpha/Iyys		% 1/sec*sec
Malphadot=qbar*S*cbar*cbar*Cmalphadot/(2*Iyys*U1)	% 1/sec
Mq=qbar*S*cbar*cbar*Cmq/(2*Iyys*U1)		% 1/sec
Mdeltae=qbar*S*cbar*Cmdeltae/Iyys		% 1/sec*sec

%
% Computation of Primed Dimensional Stability and Control Derivatives
%	
XuP=Xu+XTu;
XalphaP=Xalpha;
XthetaP=-g*cos(gamma1rad);
XdeltaeP=Xdeltae;
%
ZuP=Zu/(U1-Zalphadot);
ZalphaP=Zalpha/(U1-Zalphadot);
ZqP=(Zq+U1)/(U1-Zalphadot);
ZthetaP=-g*sin(gamma1rad)/(U1-Zalphadot);
ZdeltaeP=Zdeltae/(U1-Zalphadot);
%
MuP=Mu+MTu+Malphadot*Zu/(U1-Zalphadot)
MalphaP=Malpha+MTalpha+Malphadot*Zalpha/(U1-Zalphadot)
MqP=Mq+Malphadot*(Zq+U1)/(U1-Zalphadot)
MthetaP=-Malphadot*g*sin(gamma1rad)/(U1-Zalphadot)
MdeltaeP=Mdeltae+Malphadot*Zdeltae/(U1-Zalphadot)
%
% Assemble the A and B matrices
% X=[u(ft/sec) alpha(rad) q(rad/sec) theta(rad)]'
%
A=[XuP XalphaP   0 XthetaP;
   ZuP ZalphaP ZqP ZthetaP;
   MuP MalphaP MqP MthetaP;
     0       0   1       0];
B=[XdeltaeP ZdeltaeP MdeltaeP 0]';
C=[1  0     0   0;
   0  r2d   0   0;
   0  0    r2d  0;
   0  0     0  r2d];
D=[0 0 0 0]';
%
% Analyze the linear equations of motion
%
sys=ss(A,B,C,D);
set(sys,'statename',{'u(f/s)' 'alpha(r)' 'q(r/s)' 'theta(r)'})
set(sys,'inputname','deltae(r)')
set(sys,'outputname',{'u(f/s)' 'alpha(d)' 'q(d/s)' 'theta(d)'})
sys
[Wn,Z,Poles]=damp(sys)
tfsys=tf(sys)
zpksys=zpk(sys)
