%Aerodynamics Calculations 
%Team Orion-Fall 2000
%Bridget Fitzpatrick
%calculations based on a bi-plane design

clear all
close all

%input givens
span=6.6; %span-ft
chord=1.5; %chord-ft
df=(0.417+.5)/2; %equivalent fuselage diameter-ft
lf=5.85; %fuselage length-ft 
SHT=3.33; %area of the horizontal tail-ft
SVT=1.3; %area of the vertical tail-ft  
rho=0.002308; %this is the density at sea level-slug/ft^3
Sref=19.8; %ft^2
W=10.2; %lb

%defining velocities
Vstall=20; %this is the stall Velocity-ft/s
Vcruise=25;
Vmax=30;

%defining elevator deflections
deltaE1=-15*(pi/180); %elevator deflection-rad
deltaE2=0; %elevator deflection-rad
deltaE3=15*(pi/180); %elevator deflection-rad
Ee=.4; %this is the flap to chord ratio of the elevator
eta_h=0.9; %this is the dynamic pressure ratio of the horizontal tail and the wing
tau_e=0.6; %this is the efficiency of the elevator. Calculated from Roskams 2-D graph p.61

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculations using Warner's Aerodynamic Design Performance paper
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

span_l=span; %this portion will set the upper and lower wing span
span_u=span; %equal to each other.
Gap=.15*span;

%Warner states that the optimum distance between the the two wings is .15*span of upper wing
x=Gap/span;
sigma=63.112*(x^4) - 62.595*(x^3) + 23.759*(x^2) - 5.2979*x + 0.9984; %this value is found from Airplane Performance Design paper

%excel was used to create a fourth order curve fit
Ru=span/chord; %AR of upper wing
Rl=span/chord; %AR of lower wing. these two values are =
R=Ru; %setting the AR = to the AR from one wing
chord_l=chord; %chord of the lower wing
chord_u=chord; %chord of the upper wing. these two values are =
Su=Ru*(chord_u^2);
S=R*((chord_u^2)+(chord_l^2));

%equations apply only when there is no declage
n=1+ (2*sigma*chord_u*chord_l)/((chord_u^2)+(chord_l^2));
ncheck=1+sigma; %this should be equal to n
y=Gap/chord;
beta=0.072*(y^2)-0.226*y+0.208; %taken from text of Airplane Performance

%excel was used to create a second order curve fit
%lift curve slope
alpha=[-6.94 -5.95 -4.91 -3.9 -2.86 -1.85 -0.84 0.2 1.22 2.24 3.27 4.29 5.31 6.34 7.36 8.37 9.39 10.41 11.41 12.42 13.42 14.43].*(pi/180); %needs to be set before can figure Clmax. this is the experimental value for stall given
alphadeg=alpha*(180/pi);
dClslope=0.0986*(180/pi);%this is a theoretical value taken from lift curve slope. Re of stall was used.-rad
dCl_slopeb=1/((1/dClslope)+(n/(pi*R))+(beta/2));%taken from airplane performance design paper

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Values using Roskam-Volume 6, Chapter 8
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Calculating Alpha Zero Lift of the Wing
value=-0.42; %found using Figure 8.41. A taper ratio of 1 and 0 quarter-chord sweep angle
%delta_ao/eps_t=value;  
i_w=(2.2)*(pi/180); %this is found from finding Cl cruise value of .7 and using curve fit to find the angle that correlates-rad
alpha_wing=alpha-i_w; %this is the airplane aoa- the wing incidence angle-Eqn(8.20)
alpha_ol=(-4.856)*(pi/180); %Alpha Zero Lift. This value is taken from graph of Cl vs alpha for RE~200000-rad
eps_t=value*alpha_ol;
ratio= 1;%this is the ratio of alpha_ol/alpha_ol2. This had to be assumed
%where alpha_ol2 is alpha zero lift at M=0.3.
alpha_olw=(alpha_ol+(value)*eps_t)*(ratio); %this is the wing zero lift aoa-Eqn(8.21)
%delta_ao/eps_t is the change in wing zero lift aoa per degree of linear wing twist found using Figure 8.41 
%eps_t is the wing twist angle found using Figure 8.41
%alpha_ol/alpaha_ol2 is found from Figure 8.42
%Calculation of Wing Lift Curve Slope: CL_alphaw
beta=1; %assuming this b/c beta=(1-(M^2))^.5; and  M~0;
deltac_2=0; %no sweep
Cl_alpha_wt=0.0986*(180/pi); %this is the theoretical airfoil lift curve slope
%taken from curve fit of Cl vs alpha and throwing out points around stall
kw=Cl_alpha_wt/(2*pi); %this is found using M=2*pi/beta2 :248 -Eqn(8.25)
ARW=4.4; %this is the aspect ratio of the wing
CL_alphaw=2*pi*ARW/(2+(((ARW^2)/kw^2)+4)^.5); 
%this is the wing lift slope curve for conventional, straight, tapered wings w/ moderate sweep angles 
%Calculation of the Wing maximum lift coefficient
%Clw_stall ; %this is the spanwise value of the local section lift coefficient
%CLmax_w=(1/S)*int(b*c*Clw_stall,0,1);
%Airplane zero-angle-of-attack lift coefficient: CLo
%the horizontal tail was approximated by a NACA 0006. Theoretically a flat plate.
Cl_alpha_ht=(0.0871)*(180/pi);%this was obtained same way as wing CL_alpha
h_h=(7.5)/12; %distance between the horizontal tail centerline and wing center line-ft
l_v=4.1; %distance from the wing c/4 to the vertical tail c/4-ft
l_h=4.1; %distance from the wing c/4 to the horizontal tail c/4-ft
kh=(1-h_h/span)*((2*l_h/span)^(1/3));
%Cl_alpha_ht/(2*pi); %this is the theoretical ratio of the airfoil lift curve slope
ARH= 3; %this is the aspect ratio of the horizontal tail
CL_alpha_h= 2*pi*ARH/(2+(((ARH^2)/kh^2)+4)^.5);%this is the horizontal tail lift curve slope found from 8.1.3.2
Kwf=1+0.025*(df/span)-0.25*(df/span)^2; %wing-fuselage interference factor
CL_alpha_wf=Kwf*CL_alphaw; %wing-fuselage lift curve slope found from 8.1.5.3 
CLo_wf=(i_w-alpha_olw)*CL_alpha_wf; %Eqn(8.33)     
qstall=.5*rho*Vstall^2; %this is the dynamic pressure of the entire aircraft
qcruise=.5*rho*Vcruise^2; %this is the dynamic pressure of the entire aircraft
qmax=.5*rho*Vmax^2; %this is the dynamic pressure of the entire aircraft
ada_h=.9; %this is the ratio of the dynamic pressure between the horizontal tail and total dynamic pressure
%this is subject to change.  Check out p272
i_h=0; %incidence of the horizontal tail 
eps_oh=0; %is the horizontal tail downwash angle for zero airplane angle of attack
CLo= CLo_wf+ CL_alpha_h*ada_h*(SHT/S)*(i_h-eps_oh); %Eqn(8.32)
CLo_h=CL_alpha_h*ada_h*(SHT/S)*(i_h-eps_oh);
CL_deltae_h=(CL_alpha_h/pi)*((acos(1-2*Ee))+2*sqrt(Ee*(1-Ee)));
CL_h1=CLo_h + CL_alpha_h*alpha + CL_deltae_h*deltaE1; %coefficient of lift for the horizontal tail at E1
CL_h2=CLo_h + CL_alpha_h*alpha + CL_deltae_h*deltaE2; %coefficient of lift for the horizontal tail at E2
CL_h3=CLo_h + CL_alpha_h*alpha + CL_deltae_h*deltaE3;  %coefficient of lift for the horizontal tail at E3

%Calculation of airplane lift curve slope: CLalpha
taperratio=1;
Ka=(1/ARH)-1/(1+ARH^1.7);
Klambda=(10-3*taperratio)/7;
Khoriz=(1-h_h/span)/(2*l_h/span)^(1/3);
downwashtail=4.44*(Ka*Klambda*Khoriz*(cos(0)^.5)^1.19);
CLalpha= CL_alpha_wf+ CL_alpha_h*ada_h*(SHT/S)*(1-downwashtail);
CLdeltae=CL_alpha_h*eta_h*(SHT/S).*tau_e;

%overall coefficient of lift and lift calculations
CL21=CLalpha*alpha + CLo+ CLdeltae*deltaE1;
CL22=CLalpha*alpha + CLo+ CLdeltae*deltaE2;
CL23=CLalpha*alpha + CLo+ CLdeltae*deltaE3;
CLmax21=max(CL21)
CLmax22=max(CL22)
CLmax23=max(CL23)
L21stall=.5*rho*CL21*S*(Vstall^2);
L22stall=.5*rho*CL22*S*(Vstall^2);
L23stall=.5*rho*CL23*S*(Vstall^2);
L21cruise=.5*rho*CL21*S*(Vcruise^2);
L22cruise=.5*rho*CL22*S*(Vcruise^2);
L23cruise=.5*rho*CL23*S*(Vcruise^2);
L21max=.5*rho*CL21*S*(Vmax^2);
L22max=.5*rho*CL22*S*(Vmax^2);
L23max=.5*rho*CL23*S*(Vmax^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculation of CDo using Raymer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%skin friction drag coefficient
%wing
Re_perlengthstall=.002308*Vstall./3.717E-7; %Reynolds number divided by length
Re_wingstall=Re_perlengthstall.*chord;
Cf_wingstall=(0.455./((log10(Re_wingstall)).^2.58));
Re_perlengthcruise=.002308*Vcruise./3.717E-7; %Reynolds number divided by length
Re_wingcruise=Re_perlengthcruise.*chord;
Cf_wingcruise=(0.455./((log10(Re_wingcruise)).^2.58));
Re_perlengthmax=.002308*Vmax./3.717E-7; %Reynolds number divided by length
Re_wingmax=Re_perlengthmax.*chord;
Cf_wingmax=(0.455./((log10(Re_wingmax)).^2.58));

%fuselage
Re_fuselagestall=Re_perlengthstall.*lf;
Cf_fuselagestall=(0.455./((log10(Re_fuselagestall)).^2.58));
Re_fuselagecruise=Re_perlengthcruise.*lf;
Cf_fuselagecruise=(0.455./((log10(Re_fuselagecruise)).^2.58));
Re_fuselagemax=Re_perlengthmax.*lf;
Cf_fuselagemax=(0.455./((log10(Re_fuselagemax)).^2.58));
%lamdaf=(fuselage_length./fuselage_height); %fineness ratio
%Kf=1+(60./(lamdaf.^3))+(lamdaf./400);

%vertical tail
%Vertical Tail Paramters
ARv=1.3; %from Raymer 185
lamdav=0.6;
b_vtail=sqrt(ARv*SVT);
Croot_vtail=2*SVT/(b_vtail*(1+lamdav));
Ctip_vtail=lamdav*Croot_vtail;
Cbar_vtail=(2/3)*Croot_vtail*(1+lamdav+(lamdav^2))/(1+lamdav);
Re_vtailstall=Re_perlengthstall.*Cbar_vtail;
Re_vtailcruise=Re_perlengthcruise.*Cbar_vtail;
Re_vtailmax=Re_perlengthmax.*Cbar_vtail;
Cf_vtailstall=(0.455./((log10(Re_vtailstall)).^2.58));
Cf_vtailcruise=(0.455./((log10(Re_vtailcruise)).^2.58));
Cf_vtailmax=(0.455./((log10(Re_vtailmax)).^2.58));
%Kvtail=(2.*cos(sweep_qcv));

%horizontal tail
%Horizontal Tail Parameters
lamdah=0.6; %found from Raymer p85
ARh=3; %found from Raymer p85
b_htail=sqrt(ARh*SHT);
Croot_htail=2*SHT/(b_htail*(1+lamdah));
Ctip_htail=lamdah*Croot_htail;
Cbar_htail=(2/3)*Croot_htail*(1+lamdah+(lamdah^2))/(1+lamdah);
Re_htailstall=Re_perlengthstall.*Cbar_htail;
Cf_htailstall=(0.455./((log10(Re_htailstall)).^2.58));
Re_htailcruise=Re_perlengthcruise.*Cbar_htail;
Cf_htailcruise=(0.455./((log10(Re_htailcruise)).^2.58));
Re_htailmax=Re_perlengthmax.*Cbar_htail;
Cf_htailmax=(0.455./((log10(Re_htailmax)).^2.58));

%struts
Frontal_Area_St=11.*(Gap.*8.33E-3);
Factor_St=0.4;	%D/q / Frontal Area for rectangular strut estimated from Raymer p. 349
Cd_St=(Frontal_Area_St.*Factor_St)./S;
Cd_St=1.2.*Cd_St; %Raymer p. 347, estimate that accounts for mutual interference

%Landing Gear
d_tire=.208; %ft-diamter of the tire
t_tire=.04167; %ft-thickness of the tire
Frontal_Area_Lg=2.*(d_tire.*t_tire);
Factor_Lg=0.25; %D/q / Frontal Area for regular tire from Raymer p. 349
Cd_Lg=(Frontal_Area_Lg.*Factor_Lg)./S;
Cd_Lg=1.2.*Cd_Lg; %%Raymer p. 347, stimate that accounts for mutual interference

%fix this section
%form factor
ffc_fu=1.1;
ffc_w=2;
%ffc_w=(1+(.6/(.25))*(.14)+100*(.14^4));
ffc_ht=(1+(.6/(.25))*(.14)+100*(.14^4));
ffc_vt=(1+(.6/(.25))*(.14)+100*(.14^4));

%fix this section
%Interference effects on component drag
Qc_fu=1;
Qc_w=1.5;
Qc_ht=1;
Qc_vt=1;

%wetted areas
Swet_fu=11.15; %ft^2
Swet_w= 39.89; %ft^2
Swet_ht=9.2805; %ft^2
Swet_vt=9.2805; %ft^2
Swet_st=0.8387; %ft^2
Swet_lg=0.0028; %ft^2

%fuselage
Cd_fustall=Cf_fuselagestall*ffc_fu*Qc_fu*Swet_fu/Sref;
Cd_fucruise=Cf_fuselagecruise*ffc_fu*Qc_fu*Swet_fu/Sref;
Cd_fumax=Cf_fuselagemax*ffc_fu*Qc_fu*Swet_fu/Sref;

%wing
Cd_wstall=Cf_wingstall*ffc_w*Qc_w*Swet_w/Sref;
Cd_wcruise=Cf_wingcruise*ffc_w*Qc_w*Swet_w/Sref;
Cd_wmax=Cf_wingmax*ffc_w*Qc_w*Swet_w/Sref;

%horizontal tail
Cd_htailstall=Cf_htailstall*ffc_ht*Qc_ht*Swet_ht/Sref;
Cd_htailcruise=Cf_htailcruise*ffc_ht*Qc_ht*Swet_ht/Sref;
Cd_htailmax=Cf_htailmax*ffc_ht*Qc_ht*Swet_ht/Sref;

%vertical tail
Cd_vtstall=Cf_vtailstall*ffc_vt*Qc_vt*Swet_vt/Sref;
Cd_vtcruise=Cf_vtailcruise*ffc_vt*Qc_vt*Swet_vt/Sref;
Cd_vtmax=Cf_vtailmax*ffc_vt*Qc_vt*Swet_vt/Sref;

%struts
%Cd_ststall=Cfc_st*ffc_st*Qc_st*Swet_st/Sref;
%Cd_stcruise=Cfc_st*ffc_st*Qc_st*Swet_st/Sref;
%Cd_stmax=Cfc_st*ffc_st*Qc_st*Swet_st/Sref;

%landing gear
%Cd_lgstall=Cfc_lg*ffc_lg*Qc_lg*Swet_lg/Sref;
%Cd_lgcruise=Cfc_lg*ffc_lg*Qc_lg*Swet_lg/Sref;
%Cd_lgmax=Cfc_lg*ffc_lg*Qc_lg*Swet_lg/Sref;
Cdnewstall=Cd_fustall+Cd_wstall+Cd_htailstall+Cd_vtstall+Cd_St;
Cdnewcruise=Cd_fucruise+Cd_wcruise+Cd_htailcruise+Cd_vtcruise+Cd_St;
Cdnewmax=Cd_fumax+Cd_wmax+Cd_htailmax+Cd_vtmax+Cd_St;
CDostall=Cdnewstall*1.18; %margin factor added for unknowns
CDocruise=Cdnewcruise*1.18; %margin factor added for unknowns
CDomax=Cdnewmax*1.18; %margin factor added for unknowns
r=1;
mu=1;
e=((mu^2)*(1+r^2))/((mu^2)+(2*sigma*mu*r)+(r^2)); %Oswald efficiency factor for a bi-plane
CDi21=(CL21.^2)./(pi*ARW*e); %induced drag 
CDi22=(CL22.^2)./(pi*ARW*e); %induced drag 
CDi23=(CL23.^2)./(pi*ARW*e); %induced drag 
CD21stall=CDostall+CDi21;
CD21cruise=CDocruise+CDi21;
CD21max=CDomax+CDi21;
CD22stall=CDostall+CDi22;
CD22cruise=CDocruise+CDi22;
CD22max=CDomax+CDi22;
CD23stall=CDostall+CDi23;
CD23cruise=CDocruise+CDi23;
CD23max=CDomax+CDi23;
D21stall=.5*rho*S*CD21stall*(Vstall^2);
D22stall=.5*rho*S*CD22stall*(Vstall^2);
D23stall=.5*rho*S*CD23stall*(Vstall^2);
D21cruise=.5*rho*S*CD21cruise*(Vcruise^2);
D22cruise=.5*rho*S*CD22cruise*(Vcruise^2);
D23cruise=.5*rho*S*CD23cruise*(Vcruise^2);
D21max=.5*rho*S*CD21max*(Vmax^2);
D22max=.5*rho*S*CD22max*(Vmax^2);
D23max=.5*rho*S*CD23max*(Vmax^2);
k=1/(pi*ARW*e);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Wing Coefficient of Drag Prediction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Wing Zero-Lift Drag Prediction: CDo
%thicktochord=.14; %this is the thickness to chord ratio for the airfoil
%Lprime=2;%this is the airfoil thickness location parameter. Figure 4.4
%Lprime =2 for (t/c)max at Xt<30c
%Swet_w=39.89; %wetted area of the wing in ft^2
%Cf_w=.007; %this is the turbulent flat plate friction coefficient of the wing. Figure 4.3
%Problem with this b/c M=0 never intersects Re=1.5*10^5
%R_LS=0.1; %this is the lifting surface correction factor. Figure 4.2
%Problem with this b/c sweep=0 is not on the plot
%R_wf=1.05; %this is the wing fuselage interference factor. Figure 4.1
%Problem with this b/c Re of 1.5*10^5 is not on the graph
%CDo_w= R_wf*R_LS*Cf_w*(1+Lprime*(thicktochord)+ 100*(thicktochord)^4)*(Swet_w/S); %(4.2.1.1)
%Wing Drag due to Lift: CD_Lw
%CL_h=CL_alpha_h*alpha; %check this out
%CL_w=CL2-CL_h*S_h/S; %this is the wing lift coefficient
%w=1; %this is zero-lift drag factor due to linear twist. Figure 4.10 
%this is a guess. need to figure out what this should be
%Vind=.0017; %this is the induced drag factor due to linear twist. Figure 4.9
%This needs to be checked!!!!!!!!!!!!
%CD_Lw=((CL_w).^2)/(pi*ARW*e)+2*pi*CL_w*eps_t*Vind+4*(pi.^2)*(eps_t.^2)*w; %Eq(4.2.1.2)
%CD_w=CDo_w+CD_Lw; %Eq(4.5)
%CDw_wave=   ; %this is the wave drag coefficient. Figure 4.11
%Subsonic Fuselage Drag Coefficient: CDLfus
%!CDLfus=  ; %this is the fuselage drag coefficient due to lift. (4.3.1.2)
%Cf_fus=0.007; %this is the skin friction coefficient of the fuselage Figure4.3
%Swetfus=11.15; %this is the wetted area of the fuselage 
%Sbfus=
%db=sqrt((4/pi)*Sbfus); %this is the fuselage base diameter. Figure(4.17)
%CDofus-base=
%Sfus=
%CDbfus=.029*((db/df)^3)/((CDofus-base*(S/Sfus))^.5)*(Sfus/S); %this is the fuselage base drag coefficient. Eq(4.32)
%CDofus=R_wf*Cf_fus*(1+60/((lf/df)^3)+.0025*(lf/df))*Swetfus/S% +CDbfus; %this is the fuselage zero-lift drag coefficient. (4.3.1.1)

%Fuselage Drag Coefficient due to lift: CDLfus
%ada=   ; %drag ratio. Figure 
%CDLfus=ada*Cdc*(alpha^3)*Splf_fus/S;
%CDfus=CDofus+CDLfus; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% More from NACA paper
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%coefficient of lift and lift calculations
CL11stall= dCl_slopeb*alpha + CLo +(CLdeltae*deltaE1);
CL11cruise= dCl_slopeb*alpha + CLo +(CLdeltae*deltaE1);
CL11max= dCl_slopeb*alpha + CLo +(CLdeltae*deltaE1);
CL12stall=dCl_slopeb*alpha + CLo +(CLdeltae*deltaE2);
CL12cruise=dCl_slopeb*alpha + CLo +(CLdeltae*deltaE2);
CL12max=dCl_slopeb*alpha + CLo +(CLdeltae*deltaE2);
CL13stall=dCl_slopeb*alpha + CLo +(CLdeltae*deltaE2);
CL13cruise=dCl_slopeb*alpha + CLo +(CLdeltae*deltaE2);
CL13max=dCl_slopeb*alpha + CLo +(CLdeltae*deltaE2);
CLmax11stall=max(CL11stall);
CLmax11cruise=max(CL11cruise);
CLmax11max=max(CL11max);
CLmax12stall=max(CL12stall);
CLmax12cruise=max(CL12cruise);
CLmax12max=max(CL12max);
CLmax13stall=max(CL13stall);
CLmax13cruise=max(CL13cruise);
CLmax13max=max(CL13max);
L11stall=.5*rho*CL11stall*S*Vstall^2;
L12stall=.5*rho*CL12stall*S*Vstall^2;
L13stall=.5*rho*CL13stall*S*Vstall^2;
L11cruise=.5*rho*CL11cruise*S*Vcruise^2;
L12cruise=.5*rho*CL12cruise*S*Vcruise^2;
L13cruise=.5*rho*CL13cruise*S*Vcruise^2;
L11max=.5*rho*CLmax11max*S*Vmax^2;
L12max=.5*rho*CLmax12max*S*Vmax^2;
L13max=.5*rho*CLmax13max*S*Vmax^2;
CLavg1stall=(CL11stall+CL21)/2;
CLavg2stall=(CL12stall+CL22)/2;
CLavg3stall=(CL13stall+CL23)/2;
CLavg1cruise=(CL11cruise+CL21)/2;
CLavg2cruise=(CL12cruise+CL22)/2;
CLavg3cruise=(CL13cruise+CL23)/2;
CLavg1max=(CL11max+CL21)/2;
CLavg2max=(CL12max+CL22)/2;
CLavg3max=(CL13max+CL23)/2;

%coefficient of drag and drag calculations
CDi11stall=((CL11stall.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi11cruise=((CL11cruise.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi11max=((CL11max.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi12stall=((CL12stall.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi12cruise=((CL12cruise.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi12max=((CL12max.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi13stall=((CL13stall.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi13cruise=((CL13cruise.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CDi13max=((CL13max.^2)./(pi.*S)).*((chord_u.^2)+(chord_l.^2)+2.*sigma.*chord_l.*chord_u);
CD11stall=CDi11stall+CDostall;
CD11cruise=CDi11cruise+CDocruise;
CD11max=CDi11max+CDomax;
CD12stall=CDi12stall+CDostall;
CD12cruise=CDi12cruise+CDocruise;
CD12max=CDi12max+CDomax;
CD13stall=CDi13stall+CDostall;
CD13cruise=CDi13cruise+CDocruise;
CD13max=CDi13max+CDomax;
D11stall=.5*rho*CD11stall*S*Vstall^2;
D12stall=.5*rho*CD12stall*S*Vstall^2;
D13stall=.5*rho*CD13stall*S*Vstall^2;
D11cruise=.5*rho*CD11cruise*S*Vcruise^2;
D12cruise=.5*rho*CD12cruise*S*Vcruise^2;
D13cruise=.5*rho*CD13cruise*S*Vcruise^2;
D11max=.5*rho*CD11max*S*Vmax^2;
D12max=.5*rho*CD12max*S*Vmax^2;
D13max=.5*rho*CD13max*S*Vmax^2;
CLavg1stall=(CL11stall+CL21)/2;
CLavg2stall=(CL12stall+CL22)/2;
CLavg3stall=(CL13stall+CL23)/2;
CLavg1cruise=(CL11cruise+CL21)/2;
CLavg2cruise=(CL12cruise+CL22)/2;
CLavg3cruise=(CL13cruise+CL23)/2;
CLavg1max=(CL11max+CL21)/2;
CLavg2max=(CL12max+CL22)/2;
CLavg3max=(CL13max+CL23)/2;
CDavg1cruise=(CD11cruise+CD21cruise)/2;
CDavg2cruise=(CD12cruise+CD22cruise)/2;
CDavg3cruise=(CD13cruise+CD23cruise)/2;
CDavg1max=(CD11max+CD21max)/2;
CDavg2max=(CD12max+CD22max)/2;
CDavg3max=(CD13max+CD23max)/2;
CDavg1stall=(CD11stall+CD21stall)/2;
CDavg2stall=(CD12stall+CD22stall)/2;
CDavg3stall=(CD13stall+CD23stall)/2;
Davg1stall=(D11stall+D21stall)/2;
Davg1cruise=(D11cruise+D21cruise)/2;
Davg1max=(D11max+D21max)/2;
Davg2stall=(D12stall+D22stall)/2;
Davg2cruise=(D12cruise+D22cruise)/2;
Davg2max=(D12max+D22max)/2;
Davg3stall=(D13stall+D23stall)/2;
Davg3cruise=(D13cruise+D23cruise)/2;
Davg3max=(D13max+D23max)/2;
Lavg1stall=(L11stall+L21stall)/2;
Lavg1cruise=(L11cruise+L21cruise)/2;
Lavg1max=(L11max+L21max)/2;
Lavg2stall=(L12stall+L22stall)/2;
Lavg2cruise=(L12cruise+L22cruise)/2;
Lavg2max=(L12max+L22max)/2;
Lavg3stall=(L13stall+L23stall)/2;
Lavg3cruise=(L13cruise+L23cruise)/2;
Lavg3max=(L13max+L23max)/2;
CLmaxavg1stall=(CLmax11stall+CLmax21)/2;
CLmaxavg2stall=(CLmax12stall+CLmax22)/2;
CLmaxavg3stall=(CLmax13stall+CLmax23)/2;
CLmaxavg1cruise=(CLmax11cruise+CLmax21)/2;
CLmaxavg2cruise=(CLmax12cruise+CLmax22)/2;
CLmaxavg3cruise=(CLmax13cruise+CLmax23)/2;
CLmaxavg1max=(CLmax11max+CLmax21)/2;
CLmaxavg2max=(CLmax12max+CLmax22)/2;
CLmaxavg3max=(CLmax13max+CLmax23)/2;
LtoD1stall=Lavg1stall./Davg1stall;
LtoD1cruise=Lavg1cruise./Davg1cruise;
LtoD1max=Lavg1max./Davg1max;
LtoD2stall=Lavg2stall./Davg2stall;
LtoD2cruise=Lavg2cruise./Davg2cruise;
LtoD2max=Lavg2max./Davg2max;
LtoD3stall=Lavg3stall./Davg3stall;
LtoD3cruise=Lavg3cruise./Davg3cruise;
LtoD3max=Lavg3max./Davg3max;
%percent_difference1=((CLmax21-CLmax11)/CLmax21)*100;
%percent_difference2=((CLmax22-CLmax12)/CLmax22)*100;
%percent_difference1=((CLmax23-CLmax13)/CLmax23)*100;

%Calculating airfoil zero-lift pitching moment coefficient: CMo
x_ac=.23*chord; %Roskam assumes that the aerodynamic center of a biplane is located at

Cmw=[0.0222
0.0181
0.0151
0.0145
0.0129
0.0121
0.0117
0.0113
0.0115
0.0123
0.013
0.0139
0.0155
0.0171
0.0191
0.0207
0.0227
0.0261
0.0298
0.04
0.0525
0.0602];

%these are the 2-D moments for the SD 7062 taken from the nasg website
c_h=1.055;
S_w=20.091;
Cmh=0;
epso=3*(pi/180);
xac_w=.833;
addh=.9;
x_ach=3.566;
CLoh=0;
Kfimp=.005; %this is the imperical pitching moment factor
Wf=.417; %width of the fuselage in ft
x=(xac_w*(CL_alpha_h/CL_alphaw)*addh*(SHT/S_w)*x_ach*(1-(downwashtail)))/(1+(CL_alpha_h/CL_alphaw)*addh*(SHT/S_w)*x_ach*(1-(downwashtail)));
SM=0.18; %static margin
xcg=x-SM; %location of the center of gravity 
Cmh_deltae=(CL_alpha_h/pi)*(1-Ee)*sqrt((1-Ee)*Ee);
Cm_deltaE=((SHT*c_h)/(S_w*chord))*Cmh_deltae-((l_h*SHT)/(chord*S_w))*CL_alpha_h*(1-downwashtail);
Cmo=Cmw+(c_h*SHT)/(chord*S_w)*Cmh+(x/chord)*CL_alphaw-(l_h*SHT)/(chord*S_w)*(CLoh+CL_alpha_h*(i_h-epso));
Cm_alpha=(x/chord)*CL_alphaw-((l_h*SHT)/(chord*S_w))*CL_alpha_h*(1-downwashtail)
Cm_fus=((Kfimp*(Wf^2)*lf)/chord*S_w)*alphadeg; %Raymer eqn 16.24
Cm1=Cmo'+(Cm_alpha*alpha)+Cm_deltaE*deltaE1;
Cm2=Cmo'+(Cm_alpha*alpha)+Cm_deltaE*deltaE2;
Cm3=Cmo'+(Cm_alpha*alpha)+Cm_deltaE*deltaE3;

%Calculating the Cmcg using Raymer
Cmcg1stall=CLavg1stall'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h1)'*(x_ach-xcg);
Cmcg1cruise=CLavg1cruise'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h1)'*(x_ach-xcg);
Cmcg1max=CLavg1max'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h1)'*(x_ach-xcg);
Cmcg2stall=CLavg2stall'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h2)'*(x_ach-xcg);
Cmcg2cruise=CLavg2cruise'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h2)'*(x_ach-xcg);
Cmcg2max=CLavg2max'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h2)'*(x_ach-xcg);
Cmcg3stall=CLavg3stall'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h3)'*(x_ach-xcg);
Cmcg3cruise=CLavg3cruise'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h3)'*(x_ach-xcg);
Cmcg3max=CLavg3max'*((xcg-xac_w)/chord)+Cmw+Cm_fus'-addh*(SHT/S_w)*(CL_h3)'*(x_ach-xcg);

%Calculating Endurance and Loiter Velocity
endu1stall=(CLavg1stall.^(3/2)./CDavg1stall);
number1stall= max(endu1stall);
term11stall=(CLavg1stall.^(3/2))-number1stall.*k.*(CLavg1stall.^2);
term21stall=number1stall*CDostall;
CL_end1stall=CLavg1stall(:,15);
V_endstall=sqrt((2*(W/S))/(rho*CL_end1stall));
istall=find(endu1stall==max(endu1stall));
alpha_max_end1stall=alpha(istall);
CL_max_end1stall=CLavg1stall(istall);
CD_max_end1stall=CDavg1stall(istall);
endu1cruise=(CLavg1cruise.^(3/2)./CDavg1cruise);
number1cruise= max(endu1cruise);
term11cruise=(CLavg1cruise.^(3/2))-number1cruise.*k.*(CLavg1cruise.^2);
term21cruise=number1cruise*CDocruise;
CL_end1cruise=CLavg1cruise(:,15);
V_endcruise=sqrt((2*(W/S))/(rho*CL_end1cruise));
icruise=find(endu1cruise==max(endu1cruise));
alpha_max_end1cruise=alpha(icruise);
CL_max_end1cruise=CLavg1cruise(icruise);
CD_max_end1cruise=CDavg1cruise(icruise);
endu1max=(CLavg1max.^(3/2)./CDavg1max);
number1max= max(endu1max);
term11max=(CLavg1max.^(3/2))-number1max.*k.*(CLavg1max.^2);
term21max=number1max*CDomax;
CL_end1max=CLavg1max(:,15);
V_endmax=sqrt((2*(W/S))/(rho*CL_end1max));
imax=find(endu1max==max(endu1max));
alpha_max_end1max=alpha(imax);
CL_max_end1max=CLavg1max(imax);
CD_max_end1max=CDavg1max(imax);
LoverD_max1stall=max(LtoD1stall);
LoverD_max1max=max(LtoD1max);
LoverD_max1cruise=max(LtoD1cruise);
endu2stall=(CLavg2stall.^(3/2)./CDavg2stall);
number2stall= max(endu2stall);
term12stall=(CLavg2stall.^(3/2))-number2stall.*k.*(CLavg2stall.^2);
term22stall=number2stall*CDostall;
CL_end2stall=CLavg2stall(:,15);
V_endstall=sqrt((2*(W/S))/(rho*CL_end2stall));
istall=find(endu2stall==max(endu2stall));
alpha_max_end2stall=alpha(istall);
CL_max_end2stall=CLavg2stall(istall);
CD_max_end2stall=CDavg2stall(istall);
LoverD_max2stall=max(LtoD2stall);
endu2cruise=(CLavg2cruise.^(3/2)./CDavg2cruise);
number2cruise= max(endu2cruise);
term12cruise=(CLavg2cruise.^(3/2))-number2cruise.*k.*(CLavg2cruise.^2);
term22cruise=number2cruise*CDocruise;
CL_end2cruise=CLavg2cruise(:,15);
V_endcruise=sqrt((2*(W/S))/(rho*CL_end2cruise));
icruise=find(endu2cruise==max(endu2cruise));
alpha_max_end2cruise=alpha(icruise);
CL_max_end2cruise=CLavg2cruise(icruise);
CD_max_end2cruise=CDavg2cruise(icruise);
LoverD_max2cruise=max(LtoD2cruise);
endu2max=(CLavg2max.^(3/2)./CDavg2max);
number2max= max(endu2max);
term12max=(CLavg2max.^(3/2))-number2max.*k.*(CLavg2max.^2);
term22max=number2max*CDomax;
CL_end2max=CLavg2max(:,15);
V_endmax=sqrt((2*(W/S))/(rho*CL_end2max));
imax=find(endu2max==max(endu2max));
alpha_max_end2max=alpha(imax);
CL_max_end2max=CLavg2max(imax);
CD_max_end2max=CDavg2max(imax);
LoverD_max2max=max(LtoD2max);
endu3stall=(CLavg3stall.^(3/2)./CDavg3stall);
number3stall= max(endu3stall);
term13stall=(CLavg3stall.^(3/2))-number3stall.*k.*(CLavg3stall.^2);
term23stall=number3stall*CDostall;
CL_end3stall=CLavg3stall(:,15);
V_endstall=sqrt((2*(W/S))/(rho*CL_end3stall));
istall=find(endu3stall==max(endu3stall));
alpha_max_end3stall=alpha(istall);
CL_max_end3stall=CLavg3stall(istall);
CD_max_end3stall=CDavg3stall(istall);
LoverD_max3stallstall=max(LtoD3stall);
endu3cruise=(CLavg3cruise.^(3/2)./CDavg3cruise);
number3cruise= max(endu3cruise);
term13cruise=(CLavg3cruise.^(3/2))-number3cruise.*k.*(CLavg3cruise.^2);
term23cruise=number3cruise*CDocruise;
CL_end3cruise=CLavg3cruise(:,15);
V_endcruise=sqrt((2*(W/S))/(rho*CL_end3cruise));
icruise=find(endu3cruise==max(endu3cruise));
alpha_max_end3cruise=alpha(icruise);
CL_max_end3cruise=CLavg3cruise(icruise);
CD_max_end3cruise=CDavg3cruise(icruise);
LoverD_max3cruise=max(LtoD3cruise);
endu3max=(CLavg3max.^(3/2)./CDavg3max);
number3max= max(endu3max);
term13max=(CLavg3max.^(3/2))-number3max.*k.*(CLavg3max.^2);
term23max=number3max*CDomax;
CL_end3max=CLavg3max(:,15);
V_endmax=sqrt((2*(W/S))/(rho*CL_end3max));
imax=find(endu3max==max(endu3max));
alpha_max_end3max=alpha(imax);
CL_max_end3max=CLavg3max(imax);
CD_max_end3max=CDavg3max(imax);
LoverD_max3max=max(LtoD3max);

%%%%%%%%%%%%%%%%%%%%%%%%
%plotting values
%%%%%%%%%%%%%%%%%%%%%%%%

%stall

figure(1)
subplot(2,2,1)
plot(alphadeg,CL11stall,alphadeg,CL21,alphadeg,CLavg1stall)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at -15 deg elevator deflection at Stall Speed')
legend('Warner Method', 'Roskam Method','Average')

subplot(2,2,2)
plot(alphadeg,CL12stall,alphadeg,CL22,alphadeg,CLavg2stall)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at 0 deg elevator deflection at Stall Speed')
legend('Warner Method', 'Roskam Method','Average')

subplot(2,2,3)
plot(alphadeg,CL13stall,alphadeg,CL23,alphadeg,CLavg3stall)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at 15 deg elevator deflection at Stall Speed')
legend('Warner Method', 'Roskam Method','Average')

%cruise

figure(2)
subplot(2,2,1)
plot(alphadeg,CL11cruise,alphadeg,CL21,alphadeg,CLavg1cruise)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at -15 deg elevator deflection at Cruise Speed')
legend('Warner Method', 'Roskam Method','Average')

subplot(2,2,2)
plot(alphadeg,CL12cruise,alphadeg,CL22,alphadeg,CLavg2cruise)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at 0 deg elevator deflection at Cruise Speed')
legend('Warner Method', 'Roskam Method','Average')

subplot(2,2,3)
plot(alphadeg,CL13cruise,alphadeg,CL23,alphadeg,CLavg3cruise)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at 15 deg elevator deflection at Cruise Speed')
legend('Warner Method', 'Roskam Method','Average')

figure(3)
subplot(2,2,1)
plot(alphadeg,CL11max,alphadeg,CL21,alphadeg,CLavg1max)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at -15 deg elevator deflection at Max Speed')
legend('Warner Method', 'Roskam Method','Average')

subplot(2,2,2)
plot(alphadeg,CL12max,alphadeg,CL22,alphadeg,CLavg2max)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at 0 deg elevator deflection at Max Speed')
legend('Warner Method', 'Roskam Method','Average')

subplot(2,2,3)
plot(alphadeg,CL13max,alphadeg,CL23,alphadeg,CLavg3max)
xlabel('alpha (deg)')
ylabel('Coefficient of Lift')
title('Alpha vs Coefficient of Lift at 15 deg elevator deflection at Max Speed')
legend('Warner Method', 'Roskam Method','Average')

%stall

figure(4)
subplot(2,2,1)
plot(CLavg1stall,CDavg1stall)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at -15 deg elevator deflection at Stall Speed')

%cruise

subplot(2,2,2)
plot(CLavg1cruise,CDavg1cruise)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at -15 deg elevator deflection at Cruise Speed')

%max

subplot(2,2,3)
plot(CLavg1max,CDavg1max)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at -15 deg elevator deflection at Max Speed')

%stall

figure(5)
subplot(2,2,1)
plot(CLavg2stall,CDavg2stall)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at 0 deg elevator deflection at Stall Speed')

%cruise

subplot(2,2,2)
plot(CLavg2cruise,CDavg2cruise)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at 0 deg elevator deflection at Cruise Speed')

%max

subplot(2,2,3)
plot(CLavg2max,CDavg2max)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at 0 deg elevator deflection at Max Speed')

%stall

figure(6)
subplot(2,2,1)
plot(CLavg3stall,CDavg3stall)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at 15 deg elevator deflection at Stall Speed')

%cruise

subplot(2,2,2)
plot(CLavg3cruise,CDavg3cruise)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at 15 deg elevator deflection at Cruise Speed')

%max

subplot(2,2,3)
plot(CLavg3max,CDavg3max)
xlabel('CLavg')
ylabel('CDavg')
title('Drag Polar at 15 deg elevator deflection at Max Speed')

figure(7)
subplot(2,2,1)
plot(alphadeg,LtoD1max)
xlabel('alpha (deg)')
ylabel('Lavg/Davg')
title('Lift to Drag Ratio at -15 deg elevator deflection')

subplot(2,2,2)
LtoD2max=Lavg2max./Davg2max;
plot(alphadeg,LtoD2max)
xlabel('alpha (deg)')
ylabel('Lavg/Davg')
title('Lift to Drag Ratio at 0 deg elevator deflection')

subplot(2,2,3)
LtoDstall=Lavg3stall./Davg3stall;
plot(alphadeg,LtoD3stall)
xlabel('alpha (deg)')
ylabel('Lavg/Davg')
title('Lift to Drag Ratio at 15 deg elevator deflection')


figure(8)
subplot(2,2,1)
plot(alpha,endu1stall)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at -15 deg elevator deflection at Stall Speed')

subplot(2,2,2)
plot(alpha,endu1cruise)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at -15 deg elevator deflection at Cruise Speed')

subplot(2,2,3)
plot(alpha,endu1cruise)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at -15 deg elevator deflection at Cruise Speed')

figure(9)
subplot(2,2,1)
plot(alpha,endu2stall)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at 0 deg elevator deflection at Stall Speed')

subplot(2,2,2)
plot(alpha,endu2cruise)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at 0 deg elevator deflectionat Cruise Speed')

subplot(2,2,3)
plot(alpha,endu2max)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at 0 deg elevator deflection at Max Speed')

figure(10)
subplot(2,2,1)
plot(alpha,endu3stall)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at 15 deg elevator deflection at Stall Speed')

subplot(2,2,2)
plot(alpha,endu3cruise)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at 15 deg elevator deflection at Cruise Speed')

subplot(2,2,3)
plot(alpha,endu3max)
xlabel('alpha (deg)')
ylabel('Endurance Parameter')
title('Endurance at 15 deg elevator deflection at Max Speed')

figure(11)
plot(alphadeg,Cm1,alphadeg,Cm2,alphadeg,Cm3)
xlabel('alpha (deg)')
ylabel('Coefficient of Moment')
title('Moment Coefficient vs angle of attack')
legend('flap deflection=-15','flap deflection=0','flap deflection=15')

figure(12)
plot(alphadeg,Cmcg1stall,alphadeg,Cmcg2stall,alphadeg,Cmcg3stall)
xlabel('alpha (deg)')
ylabel('Cmcg')
title('Coefficient of Moment about the Cg at Stall Speed')
