% Geometry of Boiler Express to be used as input for stability.m 
% The units are english
% Variable syntax:
% eg) CLa_w				eg) AR_h
%		a is alpha			_h is horizontal tail
%		_w is wing
%
% some variables originally used were wrong, and compensated for, the names
% were not necessarily changed.
%
% ALL code altered by Jake Grasmick, Jeremy Davis, Dan Squibb, Sean McElroy, Andy Wiratama
% and Amanda Allen

%
biplane=0;				%is it a biplane? 1=yes, 0=no
TVT = 1;					%if plane has twin vertical tails, TVT=1
%if not then TVT=0

                    
W=10.25; %weight of aircraft [lbs] (from report)
U1=28;	%flight speed [ft/s] (from report)
rho=0.002308; %density [slug/ft^3] (standard day conditions)
q=.5*rho*U1^2; %dynamic pressure	[(slug ft)/(s2 ft2) or lbf/ft2]
gam=1.4; %gamma (standard day air conditions)
R=1718;	%gas constant air (standard day conditions)
T=536.4; %air temperature R (standard day conditions)
a1=sqrt(gam*R*T); %speed of sound  [ft/s] 
M=U1/a1;			%Mach speed
Beta=sqrt(1-M^2); %compressibility correction factor
dCDdM=0;		%change in coefficient of drag w.r.t. change in Mach number @ const. CL 
 				%from Fig. 4.1 in Roskam "Methods for Estimating Stability and Control Derivatives
             %of Conventional Subsonic Airplanes" 1977.
% below said zero due to it is described as the change in position which tend to occur in the high subsonic range, 
%Boiler Expresse's is only 28 ft/s  (txt pg 137)
dXacdm=0;	%change in aerodynamic center location w.r.t. Mach number @ const. CL (same ref. as above)

a=0*pi/180;	%wing angle of attack at cruise steady level flight [rad]
d=.3 ;    %fuselage diam or width at the wing root cord - measured (ft)
deda=.33;   %from report

% fuselage
xLE=1.25;            %distance from LE to nose of fuselage [ft] (measured)
xTOT=3;              %total length of fuselage [ft] (measured)
w=.29167;			   %average width of fuselage [ft] (report)
h=0.5;					%average height of fuselage [ft] (report)
l=3;                 %average length of fuselage [ft] (measured)
avgAREA=4*l*h+2*w*h; %average area of fuselage [ft^2]
dboom=0.04583;       %average diameter (width) of tail boom [ft] (report)
hboom=2.833;         %average length of tail boom [ft] (report)
tboom=0.00267;       %average thickness of tail boom [ft] (report)
avgAREAboom=2*(pi*(dboom/2)*hboom+2*(pi*(dboom/2)^2)); %average area of tail booms [ft^2] 
							%boom area negligble compared to fuselage area                     
                     
lb=xTOT;					%fig 7.19
KN=.0033;				%fig 7.19
KRl=1.2;					%fig 7.20 (ask Andrisani)

% wing
b_w=12.84;					%wing span	[ft] (measured)
cr_w=1.395833;				%root chord of the wing [ft] (measured)
ct_w=0.95833;				%tip chord of the wing [ft] (measured)
c_w=(cr_w + ct_w)/2;	   %mean aerodynamic chord of wing [ft]
Gam_w=0;					   %dihedral of wing (rad)
Xcg=0.35*c_w;           %CG position [ft] 

SBS=(h*xLE-(h-.2)*xLE)+h*(c_w-xLE)+h/2*(lb-(xLE+c_w));  
							%body side area, different eqn for different aircraft ft2
							%couldn't be determined before knowing c_w
							%SBS from fig 7.19



S_w=c_w*b_w;			%wing area (ft2)
%S_w=13.5
if biplane ==1
	AR_w=2*b_w^2/S_w;		%aspect Ratio of wing
	gap=1;					%gap between biplane wings (ft)
	stag=0;					%stagger (ft)
	Y=4/7*gap;				%vertical position from bottom plane of biplane ft
	Z_w=h/2-Y;			%vertical position of wingroot quartercord point to 
   						%fuselage centerline, positive downward
   						%chapter 7.1
else 
   [AR_w]=b_w^2/S_w;
   Y=0.6;					%vertical position from bottom of airplane [ft]
   Z_w=h/2-Y;				
end

CL_w=W/(q*S_w);					%3d wing coefficient of lift 
%Cla_w=.0986*180/pi;	%2d d(Cl)/d(alpha) (from report) 1/rad

%CLa_w calculated toooo many times, comment this one out and use
% calculated one

%CLa_w=0.044*180/pi;	%3d d(CL)/d(alpha) of wing fig 3.3 for monoplane not biplane

%CLa=  ;					%3d d(Cl)/d(alpha) (from report) 1/rad - CALCULATED NOT IN BE REPORT
CDo=0.0225;				%from plot 3.5, in report
Cmadot_w=0;				%from chapter 6.3 NEGLIGBLE
CnB_w=0;					%from chapter 7.3 - NEGLECT UNLESS HIGH AOA


CL_wb=CL_w;				%from chapter 7.2
ClB_CL2=0;				%fig 7.11 - converges at zero
KML=1;					%fig 7.12 - M SMALL
Kf=.935;					%fig 7.13 - estimated for A = approx 10
ClB_CL_A=-.0005;		%fig 7.14
ClB_Gam=0.000264;		%fig 7.15
KMG=1;					%fig 7.16 - est with SMALL M
d_l=sqrt(avgAREA/.7854); %eqn 7.11
dClB_Gam=-.0005*sqrt(AR_w)*(d_l/b_w)^2*(180/pi)^2; %eqn 7.10
dClB_Zw=-1.2*sqrt(AR_w)*Z_w/b_w*2*d_l/b_w; %eqn 7.12
theta_w=0;					%wing twist between root and tip sections, 
							%negative for washout shown in fig 7.17
dClB_thetatanL4=0;	%wing twist correction factor obtained from fig 7.17
BClp_Kap=-.575;		%from fig 8.1   - approximated                



lam_w=ct_w/cr_w;		%tip-to-chord ratio of wing
Lam_w=0; 				%LE sweep of wing     rad
Lamq_w=0; 				%quarter-chord sweep of wing    rad
B=sqrt(1-M^2*(cos(Lamq_w))^2); %eqn from symbol definition
Lam2_w=0;				%half-chord sweep rad from fig 2.5

Cla_w = 5.44 ;       %from report and Cl vs alpha plot for S1210  airfoil
Kap=Cla_w/(2*pi);
Xac_w=.3*c_w;			%from report


dCnpTheta=0.0002; 	%fig8.2
Theta=0;					%fig 8.2 - no twist
dCnp_add=0;     		%fig 8.3 - no flaps - in paper
delf=0;					%geometry - no flaps
adelf=0;					%section 10.1

ClrCL_M0=.26;			%fig 9.1
dClrTheta=-.0175;		%fig 9.2
Gam=0;					%geometry
Theta=0;					%geometry
dClradelfdelf=0;		%fig 9.3  - NO FLAPS
CnrCL2=-.025;			%fig 9.4 - A= approx 10, lam=.69
CnrCDo=-0.3;	   	%fig 9.5 - with cg to ac from report= ac at .46 and cg at .35

%The BoilerExpress is reported to have no flaps, thus all flap dependent values 
%  are considered to be zero.

adCLadCl=1.05;	   		%fig 10.2 - NO FLAPS - used aileron geometry
adCl=-0.65 ; 				%fig 10.2 - NO FLAPS - used aileron geometry
Kb=0.6;						%fig 10.3 - NO FLAPS - aileron geometry
Kprime=1;				%fig 10.7 - NO FLAPS - aileron geometry
Cldth=0;				   %fig 10.5 - NO FLAPS
CldCldth=0;  			%fig 10.6 - NO FLAPS

BCldPrKap=0.015;		%figure 11.1 - approximated 
CldR=0;					%right aileron deflection coeff
CldL=0;					%left aileron deflection coeff
dr=0*pi/180;			%right aileron deflection [rad]
dl=0*pi/180;			%left aileron deflection [rad]

K=-.1;					%empirical factor obrained from figure 11.3

% horizontal tail
b_h=2.71;				%h-tail span (ft)
cr_h=0.69;				%root chord (ft)
ct_h=0.69;				%tip chord (ft)
c_h=(cr_h+ct_h)/2;	%mean aero chord of h-t   (ft)
Gam_h =0;				%dihedral of h-tail (rad)
Y=0.6		;				%vertical position from bottom of airplane ft
Z_h=h/2-Y;				%vertical position of h-tailroot quartercord point to 
   						%fuselage centerline, positive downward
   						%chapter 7.1



S_h=b_h*c_h;					%h-tail area (ft2)
AR_h=b_h^2/S_h;		%Aspect Ratio of h-tail
eta_h=0.85;			   	%ratio of dynamic pressure at the tail to that of
							%the free stream  (from report)
hh=0;				   	%distance define in fig. 3.7
CLa_h=0.086;			%fig 3.3

%airfoil of horiz tail unknown, left as NACA 0006
Cla_h=.0871*180/pi;  %2d d(Cl)/d(alpha) for NACA 0006
CL_h=2*pi*Cla_h;		%linear approximation of h-tail CL
CL_hb=CL_h;				%another approximation

l_h=3.5-xLE;			%distance defined in fig 3.7

CLhb=CL_h;				%from chapter 7.2
ClB_CL2_h=0;				%fig 7.11
KML_h=1;					%fig 7.12
Kf_h=1;						%fig 7.13
ClB_CL_A_h=-.0005;		%fig 7.14
ClB_Gam_h=-0.00014;				%fig 7.15
KMG_h=1;					%fig 7.16
dClB_Gam_h=-.0005*sqrt(AR_h)*(d_l/b_h)^2*(180/pi)^2; %eqn 7.10
dClB_Zh=-1.2*sqrt(AR_h)*Z_h/b_h*2*d_l/b_h; %eqn 7.12
theta_h=0;				%wing twist between root and tip sections, 
							%negative for washout shown in fig 7.17
dClB_thetatanL4_h=0;	%wing twist correction factor obtained from fig 7.17
BClp_Kap_h=-.39;		%from fig 8.1  with AR approx 10, K .8658 and Beta .9996




lam_h=ct_h/cr_h;		%h-tail tip-to-chord ratio
Lam_h=0; 	%LE sweep of h-tail	rad
Lamq_h=0; 	%quarter-chord sweep of h-tail	rad
Lam2_h=0; 		%half-chord sweep rad from fig 2.5
Kap2=1;






% vertical tail
a_v=0;					%geometric angle of attack
b_v=0.875;
c_v=0.69;
lam_v=0.465;
Lam_v=21*pi/180;     %rad
Lamq_v=10.9*pi/180;  %rad
%S_v=1.18;
S_v=2*1.18;				% 2 times due to two tails - report gives total S_v is 2.0 ft^2
AR_v=b_v^2/S_v;		%from geometry
eta_v=1;
l_v=4.1;					%from geometry
Z_v=b_v/2;				%from geometry

%calculated in Stability, Jeremy also found a value
CLa_v=2*pi;

k_v=1; 			   	%from fig 7.3


AvbAv=1.05;				%the ratio of aspect ratio of the vertical panel in 
							%the presence of the body to that of the isolated panel. 
                     %fig 7.5
                     
AvhbAvb=1.25;			%the ratio of the vertical panel aspect ratio in the
                     %presence of the horizontal tail and body to that of the
                     %panel in the presence of the body alone fig7.6
                     
KH=1.05; 					%a factor accopunting for 
							%the relative size of the horizontal and vertical tails
							%fig 7.7
AR_v_eff=AvbAv*AR_v*(1+KH*(AvhbAvb-1)); %eqn 7.6
Lam2_v=1.105;
CyB_v_eff= 2.45;		%for use with twin vert tails fig 7.9
CyB_v_wbh_eff= 0.94;%for use with twin vert tails fig 7.10

                     