
function [constant]=Cessna182
% function [constant]=Cessna182
% Version 1.0, 2/14/01
% Make a vector of constants for the Cessna 182 aircraft
% for the cruise configuration. Ref Roskam Page 480
% The order of these constants is critical. Do not change the order
% or code in other simulation functions will be in error.
% Primary authors: Kacie Burton, Kevin Dahya, Mellisa Glaser, Kerem Koray
% Wael Nour, and Tanya Tunistra
% AAE 421 Design Project Spring 2001

disp('Input constants for the Cessna 182 at cruise, Vt=28 ft/sec, Hp=646 ft.')

%**********************************************************************
% This file contains the basic Boiler Xpress parameters and computes
% the control derivatives from the function files.
%**********************************************************************

%**********************************************************************
%
% Variable identification for the Control Derivative Derivations:
%
%**********************************************************************
global g		%ft/s^2
g=32.2;		%gravitational constant
global Vcruise		%ft/sec
Vcruise = 28;		%taken from 451 Boiler Xpress report
global q		%dynamic pressure
q = 0.5*.00231*Vcruise^2;	
global n	%load factor for manuevering
n=1;
global Fp	% assume no turning of flow by propeller
Fp=0;
global Fpa		
Fpa=0;
global B	
B=1;
global alpha			%assumed steady level flight
alpha=0;

%**********************************************************************
% 
% Geometric parameters of the aircraft
%
%**********************************************************************

global weight			%lb 
weight=10.5;			%taken from Boiler Xpress 451 final report
global lift				%lb
lift = weight;
global eo				%Oswald Efficiency factor
eo = .75;				%taken from Boiler Xpress 451 final report
global A					%Aspect ratio of wing
A=9;						%taken from Boiler Xpress 451 final report
global k1 				%constant used to find induced drag
k1 =  1/( pi*A*eo );	
%**********************************************************************
%
% Main Wing
%
%**********************************************************************
%The following geometric parameters taken from the Boiler Xpress final report.
global cbar
cbar=1.22;								%normalized chord
global croot							%root chord (ft)
croot=1.42;
global Sref								%Refrence area (ft^2)
Sref=13.5;
global b 								%wing span (ft)
b = (Sref*A)^.5;
global broot							%root span between pods - (ft)
broot=3.0;
global btip								%span of one tip panel (ft)
btip=0.5*(b-broot);
global ctip								%tip chord(ft)
ctip=11.5*(1/12);
global xcg
xcg=((croot+ctip)/2)*0.35;			%location of the center of gravity (ft)
global xbar_cg 
xbar_cg =xcg/cbar;					%cg location normalized by MAC
global xw
xw=0.1;
global Swet_wing      				%Total wetted area of wing (ft^2)
Swet_wing = Sref*2*1.04;
global Sexp
Sexp=13.3;
global lambda							%unswept wing
lambda=0;
global lambda_c4						%sweep at quarter chord
lambda_c4=0;
global taper							%taper ratio of tip section
taper=ctip/croot;
global c_wing							%mean chord of wing (ft)%== For S1210 airfoil
c_wing = Sref/b;
global x_c_wing         			%Chordwise location of maximum thickness point
x_c_wing = 0.25;
global xbar_acw						%S1210 aifoil
xbar_acw = 0.3;
global t_c_wing 		       		%Precent of thickness of the airfoil  
t_c_wing = 0.12;
global beta
beta=1;
global gamma							%dihedral angle of the wing in rads
gamma=.0349;
global theta							%no twist of the wing
theta=0;
global Cm0								%Selig data about c/4
Cm0=-0.25;
global Clmax							%Selig S1210 airfoil at Re=150,000
Clmax = 1.83;
global CLmax							%Raymer's estimate for high aspect ratio aircraft
CLmax = 0.85*Clmax;
global alpha_0L						%rad
alpha_0L = -8*pi/180;
global Cl_alpha 						%from airfoil cl vs alpha plot
Cl_alpha = 5.44;
global eta1								%sideslip of aircraft								 
eta1 = Cl_alpha/(2*pi);
global CL								%lift coefficient of aircraft
CL=weight/(q*Sref);

%**********************************************************************
%
% Tail horizontal & Vertical surfaces
%
%**********************************************************************
global xh									%distance from cg to AC for horizontal tail
xh=3.5;
global Vh									%vertical distance from a/c cg to AC for horizontal tail
Vh=.5;
global xbar_ach
xbar_ach = 2.3;
global elev_ratio							%chord of elevator/chord of h.stab
elev_ratio = 0.2;
global eta_h 								%tail dyamic pressure ratio -Raymer p. 412
eta_h = 0.85;
global b_horz								%span of horizontal stabilizer (ft)
b_horz = 3.0;
global S_horz								%Area of horizontal stabilizer (ft^2)
S_horz = 2.2;
global c_horz								%Chord of the horizontal stab (ft)
c_horz = S_horz/b_horz;
global x_c_horz							%Chordwise location of maximum thickness point
x_c_horz = 0.4;
global t_c_horz							%Precent of thickness of the airfoil            
t_c_horz = 0.06;
global Swet_horz							%wetted area of horiztonal tail (ft^2)
Swet_horz = S_horz*2*1.02;
global Q_horz 								%Component Interference factor
Q_horz = 1.44;
global c_ver_root							%Root Chord of the vertical stab.(ft)
c_ver_root = 1.0;
global c_ver_tip							%Tip chord of vertical stab (ft)
c_ver_tip = 0.6;
global c_ver								%average chord of vertical tail
c_ver = (c_ver_root+c_ver_tip)/2;
global S_ver 								%area of both vertical tails (ft^2)
S_ver = 1.75;
global h_ver 
h_ver = 0.5*S_ver/(0.5*(c_ver_root+c_ver_tip));
global Swet_ver							%wetted area of veritcal tail (ft^2) (BOTH)
Swet_ver = S_ver*2*1.02;
global x_c_ver                		%Chordwise location of maximum thickness point
x_c_ver = 0.4;
global t_c_ver 							%Precent of thickness of the airfoil 
t_c_ver = 0.06;
global A_horz								%aspect ratio of horizontal tail
A_horz=b_horz^2/S_horz;
global h_H									%height difference between the aero center of horz tail from the cg
h_H=0;
global l_H									%length difference between the aero center of horz tail from the cg
l_H=4.125;
global Vh_bar
Vh_bar=xh/cbar*S_horz/Sref;
global Z_v									%vertical height difference from aero center of vert tail from cg
Z_v=.31833;
global l_c									%length difference 
l_c=3.33;
global l_v
l_v=40;


global CL0h									%lift coeff. at zero AOA for horizontal tail
CL0h = 0;
global Cl_alphah 							%1/rad from thin symmetrical airfoil
Cl_alphah = 3.72;
global alph_0L_h							%symmetrical h.stab section
alph_0L_h = 0;

%**********************************************************************
%
% Fuselage
%
%**********************************************************************

global c_fuselage							%length in ft
c_fuselage = 3.3;
global d_fuselage							%Diameter in ft
d_fuselage = power( 4*.2/pi, .5 );
global Swet_fuselage						%Area of the fuselage in ft^2
Swet_fuselage = 0.7*pi*d_fuselage*c_fuselage;

%**********************************************************************
%
% pod ((Nacelles)
%
%**********************************************************************

global c_pods 									%length (ft)
c_pods = 1.75;
global d_pods 									%Diameter (ft)
d_pods = power( 4*.05/pi, .5 );
global Q_pods									%Component Interference factor
Q_pods = 1.5;
global Swet_pods								%Area of a pod (ft^2) (EACH) 
Swet_pods = 0.7*pi*d_pods*c_pods;

%**********************************************************************
%
% booms
%
%**********************************************************************

global c_boom									%length in ft
c_boom = 3;
global d_boom 									%Diameter in ft
d_boom = .0625;
global Swet_boom 								%Area of a boom (ft^2 ) EACH
Swet_boom = pi*d_boom*c_boom;

%**********************************************************************
%
% Stability equations Derived by Boiler Xpress 451 Team
%
%**********************************************************************

T=0;
global T 
Zbar_t = .15;
global Zbar_t									%vertical distance of thrust line from c.g.
dEu_da = 1.3;
global dEu_da									%upwash derivative with alpha from Raymer fig 16.11
dE_da = .33;
global dE_da									%downwash derivative from fig.16.12
del_h =  1-dE_da;
global del_h 									%change in horiz.stab angle of attack with wing angle of attack
F = 1.07*(1+.25/11.1)^2;
global F 										%fuselage lift factor Raymer 12.9
iw =3*pi/180;
global iw 										%wing incidence angle
ih = -6*pi/180;	
global ih 										%tail incidence angle
Cmw = Cm0*(A*(cos(lambda))^2/(A+2*cos(lambda)));
global Cmw 
delta_alpha_CLmax=2*pi/180;
global delta_alpha_CLmax					%From airfoil data
Cl_alpha_w= (2*pi*A)/(2 + (4+((A^2/eta1^2)*(1+ (tan(lambda))^2)))^.5)*Sexp/Sref*F;
global Cl_alpha_w								%wing only Raymer ch 12
alpha_CLmax = CLmax/Cl_alpha + alpha_0L + delta_alpha_CLmax;
global alpha_CLmax 
CLh_de=Cl_alphah/pi*(acos(1-2*elev_ratio) + 2*(elev_ratio*(1-elev_ratio))^.5);
global CLh_de	  								 % Mark Peters A.12
Cl_alpha = Cl_alpha_w + Cl_alphah*eta_h*(S_horz/Sref)*(1-dE_da)
global Cl_alpha 								% whole aircraft From Roskam eq 3.24
Cmfus=0;
global Cmfus 									% assume no pitching moment from fuselage
Cma_fus=0;
global Cma_fus 
xbar_p=0;
global xbar_p 
del_p=0;
global del_p
Cmq=-2.2*eta_h*(S_horz/Sref)*Cl_alphah*(xbar_ach-xbar_cg)^2;
global Cmq										%Pitch damping derivative (Raymer 16.52)
CL_de = CLh_de*(S_horz/Sref)
global CL_de 									%change in lift with elevator deflection
Cma=Cl_alpha_w*(xbar_cg - xbar_acw)+Cma_fus - eta_h*(S_horz/Sref)*Cl_alphah*(del_h)*(xbar_ach-xbar_cg)...
   +Fpa/(q*Sref)*del_p*(xbar_cg-xbar_p);
global Cma


%**********************************************************************
%
% Stability equations Derived by Boiler Xpress 421 Team
%
%**********************************************************************

%Lift Force
[Cl_alpha_w,Cl_alpha_dot,Clq,Cl_alpha_h,de_da_M,Cl_alpha_wb]=liftforce(Cl_alphah,eta_h,Vh,dE_da,xw,cbar,Cl_alpha_w,A,lambda,d_fuselage,b,beta,Cl_alpha,A_horz,h_H,b_horz,l_H,taper,xbar_acw,croot,CL)


%Pitching Moment
[Cm0,Cma_dot,Cm_de,xbar_ac,Cm_a]=pitchingmoment(Cl_alphah,eta_h,Vh,xh,cbar,dE_da,CL_de,l_H,Cl_alpha_h,Vh_bar,lambda,lambda_c4,A,B,xw,Cl_alpha_w,S_horz,Sref,xbar_ach,xbar_cg,xbar_ac,Cl_alpha,xbar_acw,Cl_alpha_wb)


%Side Force
[Cyo,Cy_Beta,Cy_da,Cy_dr,Cy_p,Cy_r,Cy_Beta_v,Cl_alpha_v,K_prime,Kb]=sideforce(gamma,Sref,S_ver,lambda,A,Z_v,b_horz,l_v,b)


%Rolling Force 
[CLo,Cl_Beta,Cl_delta_a,Cl_p,Cl_delta_r,Cl_r,ad_cl]=rollingforce(CL,gamma,Z_v,alpha,b,Cy_Beta_v,Cl_alpha_v,K_prime,Kb,S_ver,Sref,B,lambda_c4,theta,l_v,A)

%Yawing Force
[Cno,Cn_beta,Cn_delta_a,Cn_delta_r,Cnp,Cnr]=yawingforce(CL,Cl_delta_a,Cy_dr,alpha,b,Cl_alpha_w,beta,xbar_cg,xbar_ac,A,lambda,cbar,Cl_p,B,l_v,Cy_Beta_v,Z_v,lambda_c4,theta,ad_cl)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%							Reading in Constants
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Mass Related Inputs
constant(1)=weight;			%Weight
constant(2)=g;					%Gravity
constant(3)=constant(1)/constant(2);   %Mass, slugs
constant(4)=2.2176;			%Ixx slug*ft^2
constant(5)=18.7568;			%Iyy	"
constant(6)=20.7513;			%Izz   "
constant(7)=1.1805;			%%Ixz  "
constant(8)=.507;				%Total Efficiency, not just propeller
%constant(9)=;					%unassigned

%Derived constant from the inertia data
constant(10)=constant(4)*constant(6)-constant(7)*constant(7);											%gamma
constant(11)=((constant(5)-constant(6))*constant(6)-constant(7)*constant(7))/constant(10);	%cl
constant(12)=(constant(4)-constant(5)+constant(6))*constant(7)/constant(10);						%c2
constant(13)=constant(6)/constant(10);																			%c3
constant(14)=constant(7)/constant(10);																			%c4
constant(15)=(constant(6)-constant(4))/constant(5);														%c5
constant(16)=constant(7)/constant(5);																			%c6
constant(17)=1/constant(5);																						%c7
constant(18)=(constant(4)*(constant(4)-constant(5))+constant(7)*constant(7))/constant(10);	%c8
constant(19)=constant(4)/constant(10);																			%c9

%Aircraft Geometry
constant(20)=Sref;	%S, wing area, ft^2
constant(21)=cbar;	%cbar, mean geometric chord, ft
constant(22)=b;		%b, wing span, ft
constant(23)=0.08307;%phiT, thrust inclination angle, RADIANS
constant(24)=0;		%dT, thrust offset distance, ft

%Drag Polar
constant(25)=.017;  %CDmin    Value found from Fig 3.5 Final Report
constant(26)=k1;	 	%k
constant(27)=0.6;		%CLdm     Value found from Fig 3.5 Final Report

%Lift Force
constant(28)=0.95;	 %Cl0		 Fig 3.4 Final Report
constant(29)=Cl_alpha;		 %Cl_alpha
constant(30)=CL_de;    
constant(31)=Cl_alpha_dot;  
constant(32)=Clq;

%Side Force
constant(33)=Cyo;
constant(34)=Cy_Beta;
constant(35)=Cy_da;
constant(36)=Cy_dr;
constant(37)=Cy_p;
constant(38)=Cy_r;

%Rolling Moment
constant(39)=CLo;
constant(40)=Cl_Beta;
constant(41)=Cl_delta_a;
constant(42)=Cl_delta_r;
constant(43)=Cl_p;
constant(44)=Cl_r;

%Pitching Moment
constant(45)=Cm0;
constant(46)=Cm_a;
constant(47)=Cm_de;
constant(48)=Cma_dot;
constant(49)=Cmq;

%Yawing Moment
constant(50)=Cno;
constant(51)=Cn_beta;
constant(52)=Cn_delta_a;
constant(53)=Cn_delta_r;
constant(54)=Cnp;
constant(55)=Cnr;

%Reference Positions
constant(56)=xbar_cg;
constant(57)=xbar_cg;

%Trim Conditions
constant(58)=28;
constant(59)=646;
constant(60)=0;

%The constants below are used only by CessnaLongSC and CessnaLatSC
constant(61)=0;		%CLu=0
constant(62)=0;		%CDu=0
constant(63)=0;   	%CTxu
constant(64)=0;		%Cmu
constant(65)=0;		%CmTu
constant(66)=0;		%CmTalpha
constant(67)=0;		%CDdeltae