clear
close all

% Flight conditions and atmospheric parameters (sea level)
a = 1116.45;        % Speed of sound (ft/s)
rho = 0.0023769;    % Density (slugs/ft^3)
g = 32.2;           % Gravitational acceleration (ft/s^2)

Ma = 0.158;
u0 = Ma*a;
theta0 = 0;

% Moments of inertia
Ix = 1048;
Iy = 3000;
Iz = 3530;
Ixz = 0;

S = 184;    % Wing planform (ft^2)
b = 33.4;   % Wing span (ft)
c = 5.7;    % Wing mean aerodynamic chord (ft)

% Aircraft parameters
W = 2750;   % Weight (lbs)
m = W/g;    % Mass (slugs)
Pdyn = (1/2)*rho*u0^2;
CW = W/(Pdyn*S);

% Lateral-directional stability derivatives

CYbeta = -0.564;
Clbeta = -0.074;
Cnbeta = 0.071;
CYp = 0;
Clp = -0.410;
Cnp = -0.0575;
CYr = 0;
Clr = 0.107;
Cnr = -0.125;
CYdeltaa = 0;
Cldeltaa = 0.134;
Cndeltaa = -0.0035;
CYdeltar = 0.157;
Cldeltar = 0.107;
Cndeltar = -0.072;

% Lateral-directional dimensional stability derivatives

Yv = (1/2)*rho*u0*S*CYbeta;
Yp = (1/4)*rho*u0*b*S*CYp;
Yr = (1/4)*rho*u0*b*S*CYr;
Lv = (1/2)*rho*u0*b*S*Clbeta;
Lp = (1/4)*rho*u0*b^2*S*Clp;
Lr = (1/4)*rho*u0*b^2*S*Clr;
Nv = (1/2)*rho*u0*b*S*Cnbeta;
Np = (1/4)*rho*u0*b^2*S*Cnp;
Nr = (1/4)*rho*u0*b^2*S*Cnr;

xi = Ix*Iz-Ixz^2;

ScriptLv = (Ix/xi)*(Iz*Lv + Ixz*Nv);
ScriptLp = (Ix/xi)*(Iz*Lp + Ixz*Np);
ScriptLr = (Ix/xi)*(Iz*Lr + Ixz*Nr);
ScriptNv = (Iz/xi)*(Ixz*Lv + Ix*Nv);
ScriptNp = (Iz/xi)*(Ixz*Lp+ Ix*Np);
ScriptNr = (Iz/xi)*(Ixz*Lr + Ix*Nr);

ALD = [Yv/m, Yp/m, (Yr/m - u0), g*cos(theta0);
    (1/Ix)*ScriptLv, (1/Ix)*ScriptLp, (1/Ix)*ScriptLr, 0;
    (1/Iz)*ScriptNv, (1/Iz)*ScriptNp, (1/Iz)*ScriptNr, 0;
    0, 1, tan(theta0), 0];

[VLD,DLD] = eig(ALD)