disp('*** This code is written for MATLAB 5.2 ***')
%echo o
echo off
hold off
clf
%
% 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
Xref=.32;	% X-location of the reference point for the pitching moments
%           % (may serve as the center of gravity) fraction of cbar
alpha1=5;	% Angle of attack, 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
%
% Steady State Coefficients
CL1=1.64;	% Lift coefficient, nondimensional
CD1=.256;	% Drag coefficient, nondimensional
CTx1=.256;	% Thrust coefficient, nondimensional
Cm1=0;		% Aerodynamic Pitching Moment coefficient, nondimensional
CmT1=0;		% Thrust Pitching Moment coefficient, nondimensional
%
%Longitudinal Coefficients and Stability Derivatives
CD0=.043;
CDu=0;
CDalpha=1.06;
CTxu=-.6;
CL0=1.2;
CLu=.04;
CLalpha=5.04;	% Lift curve slope
CLalphadot=1.6;
CLq=4.1;
CM0REF=.047;
CmuREF=-.01;
CmalphaREF=-.66;	% Static longitudianl stabiity derivative
CmalphadotREF=-5;	% Lag of downwash stability derivative
CmqREF=-13.5;		% Swish effect in pitch
CmTu=.006;
CmTalpha=0;
%
%Longitudinal Control Derivatives
CDdeltae=0;		% Drag from elevator
CLdeltae=.4;	% Lift from elevator
Cmdeltae=-.98;	% Pitching moment from elevator
CDih=0;			% Drag from incidence of the horizontal tail
CLih=.85;		% Lift from incidence of the horizontal tail
Cmih=-2.1;		% Pitlching moment from incidence of the horizontal tail
%
% Miscellaneous inputs
g=32.17;		% Acceleration of gravity, ft/sec*sec
gamma1=0		% Steady state flight path angle, deg, see figure page 66
%					angle between horizontal and Xs axis (Wind at trim)




%
% 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;


XcgMat=0.4:.01:.5;
jj=length(XcgMat);
RPolesMat=zeros(4,jj);
IPolesMat=zeros(4,jj);
% >>>>>> start loop
for kk=1:jj
Xcg=XcgMat(kk)
% Move c.g. position fom Xref to Xcg
% Xcg=.32         % fraction of cbar
Dx=Xref-Xcg;
CM0=CM0REF-Dx*CL0;
Cmu=CmuREF-Dx*CLu;
Cmalpha=CmalphaREF-Dx*CLalpha;	% Static longitudianl stabiity derivative
Cmalphadot=CmalphadotREF-Dx*CLalphadot;	% Lag of downwash stability derivative
Cmq=CmqREF-Dx*CLq;		% Swish effect in pitch
Xac=Xref-CmalphaREF/CLalpha  % fraction of cbar
static_margin=Xac-Xcg % fraction of cbar


% 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)
RPolesMat(:,kk)=real(Poles);
IPolesMat(:,kk)=imag(Poles);
tfsys=tf(sys)
zpksys=zpk(sys)
plot(real(Poles),imag(Poles),'x')
hold on
end
% >>>>>> end loop

xlabel('Real(poles)')
ylabel('Imag(poles)')
title('LearJet (approach) Poles as a function of Xcg for Xcg=0.4:.01:.5')
strr=['Xac= ',num2str(Xac),' cbar'];
text2(.1,.8,strr)
text2(.1,.75,'A pole crosses into the LHP when Xcg is between .45 and .46')
XcgMat
RPolesMat
IPolesMat
hold off
