clc 
clear all

g = 32.174
theta1 = 0
%Cessna design parameters
S = 174 %Planform Area (ft^2)
c = 4.9 %Chord length (ft)
b = 36.0 %Wing span (ft)

%Climb Flight Condition
altitude = 0
M = 0.120
U_1 = 133.5
q_bar = 21.2
c_frac = 26.4
a_1 = 5.4

%Mass Data in Climb, Cruise, and Approach Conditions
W = 2650 %Weight (lbs)
m = W/g %Mass
I_xx_B = 948 %Moment of Inertia along the X-axis (slug-ft^2)
I_yy_B = 1346 %Moment of Inertia along the Y-axis (slug-ft^2)
I_zz_B = 1967 %Moment of Inertia along the Z-axis (slug-ft^2)
I_xz_B = 0 %Moment of Inertia along the XZ plane (slug-ft^2)

%Steady State Coefficients for Climb
CL_1 = 0.719 %Lift Coefficient
CD_1 = 0.057 %Drag Coefficient
CT_1 = 0.057 %Coefficient of Thrust
Cm_1 = 0 %Moment Coefficient
CmT_1 = 0 %Moment Coefficient due to Thrust

%Longitudinal Coefficients and Stability Derivatives
CD_0 = 0.0270
CD_u = 0
CD_a = 0.380
CT_xu = -0.171
CL_0 = 0.307
CL_u = 0
CL_a = 4.41
CL_adot = 1.7
CL_q = 3.9
Cm_0 = 0.04
Cm_u = 0
Cm_a = -0.650
Cm_adot = -5.57
Cm_q = -15.2
Cm_Tu = 0
Cm_Ta = 0

%Longitudinal Control and Hinge Moment Derivatives
CD_dele = 0
CL_dele = 0.43
Cm_dele = -1.369
Ch_a = -0.0545
Ch_del_e = -0.594

%Definition of Longitudinal, Dimensional Stability Derivatives
X_u = (-q_bar*S*(CD_u+2*CD_1))/(m*U_1)
X_a = (-q_bar*S*(CD_a-CL_1))/m
X_dele = (-q_bar*S*CD_dele)/m
Z_u = (-q_bar*S*(CL_u+2*CL_1))/(m*U_1)
Z_a = (-q_bar*S*(CL_a+CD_1))/m
Z_adot = (-q_bar*S*c*CL_adot)/(2*m*U_1)
Z_q = (-q_bar*S*c*CL_q)/(2*m*U_1)
Z_dele = (-q_bar*S*CL_dele)/m

M_u = (q_bar*S*c*(Cm_u+2*Cm_1))/(I_yy_B*U_1)
M_a = (q_bar*S*c*Cm_a)/I_yy_B
M_adot = (q_bar*S*(c^2)*Cm_adot)/(2*I_yy_B*U_1)
M_q = (q_bar*S*(c^2)*Cm_q)/(2*I_yy_B*U_1)
M_dele = (q_bar*S*c*Cm_dele)/I_yy_B

%Denominator Polynomial Approximation
D_1 = U_1*[1 -(M_q+Z_a/U_1+M_adot) (Z_a*M_q/U_1-M_a)]

%AOA to Elevetor Deflection Numerator Polynomial Approximation
Na = [Z_dele (M_dele*U_1-M_q*Z_dele)]

%Transfer Function
sys = tf(Na,D_1)
damp(sys)