Responses to Control Input for Learjet 24

The MATLAB script below takes in constants and aerodynamic data/coefficients from an airplane (in this case the Learjet 24) and plots the step/impulse responses to various control inputs such as thrust, elevator, rudder, and aileron. The constants and aerodynamic data can be swapped out for a different airplane in different conditions. Feel free to use and edit at will. If you have any input/suggestions, then please email me at gadouryjames@gmail.com.
clc
clear all
close all
%% Geometric Data
S = 230;
c_bar = 7;
b = 34;

%% Inertial Data
W = 13000;
Ix = 28000;
Iy = 18800;
Iz = 47000;
Ixz= 1300;
Ixdash = (Ix*Iz - Ixz^2)/Iz;
Izdash = (Ix*Iz - Ixz^2)/Ix;
Ixzdash= Ixz/(Ix*Iz - Ixz^2);

%% Flight Condition Data
alt = 40000;
M = 0.7;
u0 = 677;
Q = 134.6;
h = 0.32;
alpha_trim = 2.7*pi/180;

m = 404.05235;
rho = 5.87*10^-4;

%% Trim Condition
theta0 = 0;

%% Aerodynamic Coefficients for Learjet 24

%Steady State
CL1 = 0.41;
CD1 = 0.0335;
Cm1 = 0;
CTx1 = 0.0335;
CmT1 = 0;

%Stability Derivatives
CD0 = 0.0216;
CDu = 0.104;
CDalpha = 0.3;
CTxu = -0.07;
CL0 = 0.13;
CLu = 0.4;
CLalpha = 5.84;
CLalphadot = 2.2;
CLq = 4.7;
Cm0 = 0.05;
Cmu = 0.05;
Cmalpha = -0.64;
Cmalphadot = -6.7;
Cmq = -15.5;
CmTu = -0.003;
CmTalpha = 0;

%Control Derivatives
CDdeltae =0;
CDiH = 0;
CLdeltae = 0.46;
CLiH = 0.94;
Cmdeltae=-1.24;
CmiH = -2.5;

%% Non-Dimensional stability and control derivatives

Cxu = -(CDu +2*CD1);
Cxalpha = -(CDalpha - CL1);
Cxalphadot = 0;
Cxq = 0;
Cxdeltae = -CDdeltae;
Cxdeltap = CTxu + 2*CTx1;

Czu = -(CLu +2*CL1);
Czalpha = -(CLalpha +CD1);
Czalphadot = -CLalphadot;
Czq = -CLq;
Czdeltae = -CLdeltae;
Czdeltap = 0;

Cmu = Cmu +2*Cm1;
Cmdeltap = CmTu +2*CmT1;

Cw0 = W/(0.5*rho*u0^2*S);
%% Longitudinal Dimensional Derivatives

Xu = rho*u0*S*Cw0*sin(theta0) +0.5*rho*u0*S*Cxu;
Xw = 0.5*rho*u0*S*Cxalpha;
Xq = 0.25*rho*u0*c_bar*S*Cxq;
Xwdot = 0.25*rho*c_bar*S*Cxalphadot;

Zu = -rho*u0*S*Cw0*cos(theta0)+0.5*rho*u0*S*Czu;
Zw = 0.5*rho*u0*S*Czalpha;
Zq = 0.25*rho*u0*c_bar*S*Czq;
Zwdot = 0.25*rho*c_bar*S*Czalphadot;

Mu = 0.5*rho*u0*c_bar*S*Cmu;
Mw = 0.5*rho*u0*c_bar*S*Cmalpha;
Mq = 0.25*rho*u0*(c_bar^2)*S*Cmq;
Mwdot = 0.25*rho*(c_bar^2)*S*Cmalphadot;

%% Lateral Directional Derivatives

%Stability Derivatives
Clbeta = -0.11;
Clp = -0.45;
Clr = 0.16;
CYbeta = -0.73;
CYp = 0;
CYr = 0.4;
Cnbeta = 0.127;
CnTbeta= 0;
Cnp = -0.008;
Cnr = -0.2;

%Control Derivatives
Cldeltaa = 0.178;
Cldeltar = 0.019;
CYdeltaa = 0;
CYdeltar = 0.14;
Cndeltaa = -0.02;
Cndeltar = -0.074;

%% Latitudinal Directional Derivatives
Yv = 0.5*rho*u0*S*CYbeta;
Yp = 0.25*rho*u0*b*S*CYp;
Yr = 0.25*rho*u0*b*S*CYr;
Lv = 0.5*rho*u0*b*S*Clbeta;
Lp = 0.25*rho*u0*(b^2)*S*Clp;
Lr = 0.25*rho*u0*(b^2)*S*Clr;
Nv = 0.5*rho*u0*b*S*Cnbeta;
Np = 0.25*rho*u0*(b^2)*S*Cnp;
Nr = 0.25*rho*u0*(b^2)*S*Cnr;

%% Dimensional Control Derivatives
Xdeltae = Cxdeltae*0.5*rho*u0^2*S;
Xdeltap = Cxdeltap*0.5*rho*u0^2*S;
Zdeltae = Czdeltae*0.5*rho*u0^2*S;
Zdeltap = Czdeltap*0.5*rho*u0^2*S;
Mdeltae = Cmdeltae*0.5*rho*u0^2*S*c_bar;
Mdeltap = Cmdeltap*0.5*rho*u0^2*S*c_bar;

Ydeltaa = CYdeltaa*0.5*rho*u0^2*S;
Ydeltar = CYdeltar*0.5*rho*u0^2*S;
Ldeltaa = Cldeltaa*0.5*rho*u0^2*S*b;
Ldeltar = Cldeltar*0.5*rho*u0^2*S*b;
Ndeltaa = Cndeltaa*0.5*rho*u0^2*S*b;
Ndeltar = Cndeltar*0.5*rho*u0^2*S*b;
%% Calculation of A matrices
g=32.2;

Along(1,1) = Xu/m;
Along(1,2) = Xw/m;
Along(1,3) = 0;
Along(1,4) = -g*cos(theta0);
Along(1,5) = 0;
Along(2,1) = Zu/(m-Zwdot);
Along(2,2) = Zw/(m-Zwdot);
Along(2,3) = (Zq +m*u0)/(m-Zwdot);
Along(2,4) = (1/Iy)*(Mq+Mwdot*(Zq+m*u0)/(m-Zwdot));
Along(2,5) = 0;
Along(3,1) = (1/Iy)*(Mu+(Mwdot*Zu)/(m-Zwdot));
Along(3,2) = (1/Iy)*(Mw +(Mwdot*Zw)/(m-Zwdot));
Along(3,3) = (1/Iy)*(Mq+Mwdot*(Zq+m*u0)/(m-Zwdot));
Along(3,4) = -(Mwdot*m*g*sin(theta0))/(Iy*(m-Zwdot));
Along(3,5) = 0;
Along(4,1) = 0;
Along(4,2) = 0;
Along(4,3) = 1;
Along(4,4) = 0;
Along(4,5) = 0;
Along(5,1) = -sin(theta0);
Along(5,2) = cos(theta0);
Along(5,3) = 0;
Along(5,4) = -u0*cos(theta0);
Along(5,5) = 0;

Alat(1,1) = Yv/m;
Alat(1,2) = Yp/m;
Alat(1,3) = Yr/m - u0;
Alat(1,4) = g*cos(theta0);
Alat(1,5) = 0;
Alat(2,1) = Lv/Ixdash + Ixzdash*Nv;
Alat(2,2) = Lp/Ixdash + Ixzdash*Np;
Alat(2,3) = Lr/Ixdash + Ixzdash*Nr;
Alat(2,4) = 0;
Alat(2,5) = 0;
Alat(3,1) = Ixzdash*Lv + Nv/Izdash;
Alat(3,2) = Ixzdash*Lp + Np/Izdash;
Alat(3,3) = Ixzdash*Lr + Nr/Izdash;
Alat(3,4) = 0;
Alat(3,5) = 0;
Alat(4,1) = 0;
Alat(4,2) = 1;
Alat(4,3) = tan(theta0);
Alat(4,4) = 0;
Alat(4,5) = 0;
Alat(5,1) = 0;
Alat(5,2) = 0;
Alat(5,3) = sec(theta0);
Alat(5,4) = 0;
Alat(5,5) = 0;

Blong(1,1) = Xdeltae/m;
Blong(1,2) = Xdeltap/m;
Blong(2,1) = Zdeltae/(m-Zwdot);
Blong(2,2) = Zdeltap/(m-Zwdot);
Blong(3,1) = Mdeltae/Iy + Mwdot*Zdeltae/(Iy*(m-Zwdot));
Blong(3,2) = Mdeltap/Iy + Mwdot*Zdeltap/(Iy*(m-Zwdot));
Blong(4,:) = 0;
Blong(5,:) = 0;

Blat(1,1) = Ydeltaa/m;
Blat(1,2) = Ydeltar/m;
Blat(2,1) = Ldeltaa/Ixdash + Ixzdash*Ndeltaa;
Blat(2,2) = Ldeltar/Ixdash + Ixzdash*Ndeltar;
Blat(3,1) = Ndeltaa/Izdash + Ixzdash*Ldeltaa;
Blat(3,2) = Ndeltar/Izdash + Ixzdash*Ldeltar;
Blat(4,:) = 0;
Blat(5,:) = 0;

long_eig = eig(Along)
lat_eig = eig(Alat)
ts_phugoid = 4/real(abs(long_eig(3)))
ts_rollconvergence = 4/real(abs(lat_eig(3)))
% ts_roll = 4/real(abs(long_eig(1)))
% phugoid_app = tf([-2.699 -1.783 -17.36],[1.393 0.03092 0.01896])
% impulse(phugoid_app)
% step(phugoid_app)

linear_sys = ss(Along,Blong,eye(5), zeros(5,2));
C = [0 0 0 1 0];
elev_ptc = ss(Along,Blong(:,1),C,0);
C = [1 0 0 0 0];
thrust_speed = ss(Along,Blong(:,2),C,0);
C = [0 0 0 0 1];
rudder_yaw = ss(Alat,Blat(:,1),C,0);
C = [0 0 0 1 0];
ail_roll = ss(Alat,Blat(:,2),C,0);

figure,
subplot(1,2,1)
impulseplot(elev_ptc*(180/pi))
title('Impulse Response,Input: Elevator, Output:Pitch Angle')
ylabel('Amplitude (in Degrees)')
grid on

subplot(1,2,2)
stepplot(elev_ptc*(180/pi))
title('Step Response,Input: Elevator, Output:Pitch Angle')
ylabel('Amplitude (in Degrees)')
grid on

figure,
subplot(1,2,1)
impulseplot(thrust_speed)
title('Impulse Response, Input:Thrust,Output:Speed')
ylabel('Amplitude (in ft/s)')
grid on

subplot(1,2,2)
stepplot(thrust_speed)
title('Step Response, Input:Thrust,Output:Speed')
ylabel('Amplitude (in ft/s)')
grid on


figure,
subplot(1,2,1)
impulseplot(rudder_yaw*(180/pi))
title('Impulse Response, Input:Rudder, Output:Yaw angle')
ylabel('Amplitude (in Degrees)')
grid on
subplot(1,2,2)
stepplot(rudder_yaw*(180/pi))
title('Step Response, Input:Rudder, Output:Yaw angle')
ylabel('Amplitude (in Degrees)')
grid on


figure,
subplot(1,2,1)
impulseplot(ail_roll*(180/pi))
title('Impulse Response, Input:Aileron, Output:Roll Angle')
ylabel('Amplitude (in Degrees)')
grid on
subplot(1,2,2)
stepplot(ail_roll*(180/pi))
title('Step Response, Input:Aileron, Output:Roll Angle')
ylabel('Amplitude (in Degrees)')
grid on