% Patrick LeGresley
% Drag Calculations
% February 19, 1999


% Jenn too
% Modified for AIAA Design Project
% Last modified: May 1, 1999

function [CD,CL]=drag_opt(tiltwing,b_wing,S_ref,taper_wing,sweep_wing_25c,...
   tc_avg_wing,x_wing,flap_to_chord,eta_o_flap,eta_i_flap,flap_angle,slat_angle,b_ht,Sht_ref,...
   taper_ht,sweep_ht_25c,tc_avg_ht,x_ht,b_vt,Svt_ref,taper_vt,sweep_vt_25c,tc_avg_vt,...
   h,V,L_eng,D_eng,N_e,x_nosegear,gear,x_cg_aft,x_cg_forward,CL_data,Cmu_data,...
   wing_angle,plot_number,output)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs
% [CL,CD] = vector that makes up the drag polar
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fixed Wing Input Parameters

% mission profile
%h = 35000;				% altitude (ft)
%V = 350;				% speed (knots)
%W_to = 48635;				% takeoff weight (lbs)

% fuselage
Swet_fuselage = 1995;			% fuselage wetted area (sq ft)
L = 60;					% fuselage length (ft)
h_max = 9.5;				% maximum fuselage height (ft)
w_max = 13;				% maximum fuselage width (ft)
hl = 0.033;				% fuselage upsweep parameter (see AA241 course notes)
x_cg = (x_cg_aft+x_cg_forward)/2;	% Average cg location

% wing
%S_ref = 1000;				% wing reference area (sq ft)
%b_wing = 47.05;				% wing span (ft)
%taper_wing = 0.8;			% wing taper ratio
sweep_wing = sweep_wing_25c;				% wing quarter-chord sweep angle (deg)
%tc_avg_wing = 0.14;			% wing average thickness to chord ratio

% h. tail
%Sht_ref = 325;				% h. tail reference area (sq ft)
%b_ht = 29;				% h. tail span (ft)
%taper_ht = 0.8;				% h. tail taper ratio
sweep_ht = sweep_ht_25c;				% h. tail quarter-chord sweep angle (deg)
%tc_avg_ht = 0.12;			% h. tail average thickness to chord ratio

% v. tail
N_vt = 2;				% number of vertical tails
%Svt_ref = 275;				% TOTAL v. tail reference area (sq ft)
%b_vt = 12.32;				% v. tail span (ft)
%taper_vt = 0.8;				% v. tail taper ratio
sweep_vt = sweep_vt_25c;				% v. tail quarter-chord sweep angle (deg)
%tc_avg_vt = 0.1;			% v. tail taper ratio

% nacelle
%N_e = 2;                      		% Number of engines
D_nacelle_max = 1.1*D_eng;
D_nacelle_turbine = D_eng;
L_cowl_inlet = 0.6*D_nacelle_max;
L_cowl_nozzle = 0.65*D_nacelle_turbine;
L_n = L_eng + L_cowl_inlet + L_cowl_nozzle;	% Nacelle Length (ft)
D_inlet = 0.8*D_eng;                   		% Inlet Diameter (ft)
D_outlet = 0.6*D_eng;
D_exit = D_outlet;
A_max_n = (D_nacelle_max/2)^2*3.14159; % Maximum nacelle cross-sectional area (ft^2)
A_exit_n = (D_outlet/2)^2*3.14159;       % Total exit area (ft^2)
A_inlet_n = 0.7*A_max_n;      		% Inlet area (ft^2)
A_inflow_n = 0.8*A_inlet_n;   		% Effective inlet area based on mass flow (ft^2)

% pylon
tc_avg_pyl = 0.10;			% pylon average thickness to chord ratio
L_pyl = 0.8*L_n;				% pylon length (ft)
h_pyl = 1;				% pylon height (ft)
sweep_pyl = 45;				% pylon quarter-chord sweep angle (deg)

% flaps
%flap_to_chord = 0.3;               	% Flap to chord ratio  
%flap_angle = 0;                  	% Flap deflection angle (in degrees)
%eta_o_flap = 0.7;                       % Outboard flap station
%eta_i_flap = 0.25;                       % Inboard flap station
delta_CL_flap = 0;						% Change in airplane CL due to flaps (tiltwing)

% slats
slat_ratio = 1.17/25*slat_angle; % mac_wing/mac_wing_slat

% nose landing gear
%gear = 0;								% 1 if landing gear down, 0 if landing gear up
b_t_nose = 0.5;							% Nose gear tire width (ft)
d_t_nose = 1.2;							% Nose gear tire diameter (ft)
a_nose = x_nosegear;								% Distance from nose of the aircraft to nose gear (ft)
e_nose = 2.5;								% Distance from ground to fuselage (at nose gear) (ft)
p_nose = -0.3;								% Induced drag factor (0 to -0.5) (Airplane Design part VI, pg 97)

% main landing gear
b_t_main = 0.75;							% Main gear tire width (ft)
d_t_main = 1.5;							% Main gear tire diameter (ft)
a_main = 2.0;								% Width of the main gear (ft)
b_main = 2.5;								% Height of the main gear (ft)
p_main = -0.4;								% Induced drag factor (0 to -0.5) (Airplane Design part VI, pg 97)
L_main = 2.0;								% Length of the main strut (ft)
D_main = 0.4;								% Diameter of main strut (ft)

% transition locations (range from 0-1, fully turbulent to fully laminar)
xtc_wing = 0;
xtc_fus = 0;
xtc_ht = 0;
xtc_vt = 0;
xtc_nac = 0;
xtc_pyl = 0;

CDpmisc = 3;				% markup in CDp due to misc items
Cf_markup = 10;				% markup in Cf due to roughness
f_markdown = 3;				% markdown in Cf due to Mach number effects
Cf_inter = 3;				% markup in Cf due to interference effects


% end of inputs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Mission Profile Parameters
T_o = 518.7;
rho_o = 0.002377;
R = 1716;
gamma = 1.4;

% Calculate temperature and density at altitude
% Equations from Airplane Aerodynamics and Performance by Roskam and Lan

if h < 36089
   theta = 1- 6.875E-6*h;
   Temp = T_o*theta;
   rho = rho_o*theta^4.2561;
else
   Temp = -69.7 + 459.67;
   rho = rho_o*0.2971*exp(-((h-36089)/20806.7));
end

% Calculate speed of sound
a = sqrt(gamma*R*Temp);	

% Calculate viscosity coefficient
mu = (1E-10*0.3170)*Temp^1.5*(734.7/(Temp+216));

% Convert speed from knots to feet/sec
V = 1.6878*V;

% Calculate Mach number
M = V/a;

% Find weight for cruise segment (values computed inside performance)
W_beg_cr = 44962;
W_end_cr = 34910;

% Calculate operating lift coefficient at start/end of cruise segment
CL_beg_cr = W_beg_cr/(0.5*rho*V^2*S_ref);
CL_end_cr = W_end_cr/(0.5*rho*V^2*S_ref);

% end of Mission Profile Parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Lookup table data

% Lookup table to determine Cf from Re and xt/c
% from AA 241 course notes:
% http://aero.stanford.edu/aa241/drag/skinfriction.html
Reyn=[1E5 1E6 1E7 1E8 1E9;1E5 1E6 1E7 1E8 1E9;1E5 1E6 1E7 1E8 1E9;1E5 1E6 1E7 1E8 1E9;...
      1E5 1E6 1E7 1E8 1E9;1E5 1E6 1E7 1E8 1E9;1E5 1E6 1E7 1E8 1E9;1E5 1E6 1E7 1E8 1E9;...
      1E5 1E6 1E7 1E8 1E9];
trans=[0 0 0 0 0;0.1 0.1 0.1 0.1 0.1;0.2 0.2 0.2 0.2 0.2;0.3 0.3 0.3 0.3 0.3;...
      0.4 0.4 0.4 0.4 0.4;0.5 0.5 0.5 0.5 0.5;0.6 0.6 0.6 0.6 0.6;0.8 0.8 0.8 0.8 0.8;...
      1 1 1 1 1];
Cfinterp = [0.0074 0.0045 0.0031 0.0022 0.0016; 0.0072 0.0044 0.0028 0.002 0.0015;...
      0.0073 0.0042 0.0027 0.0018 0.0014;0.0070 0.0038 0.0024 0.0017 0.0013;...
      0.0067 0.0036 0.0022 0.0015 0.0012; 0.0064 0.0033 0.0019 0.0012 0.0010;...
      0.0060 0.0029 0.0017 0.0011 0.0008; 0.0052 0.0022 0.0011 0.0007 0.0005;...
      0.0042 0.0013 0.0004 0.0002 0.0001];

% Lookup table to determine Minf_perp from t/c_perp and CL/cos^2(lambda)
% from AA 241 course notes:% http://aero.stanford.edu/aa241/drag/mdiv.html
tcperp = [0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;...
      0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;...
      0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;...
      0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20;...
      0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20];
CLover = [0 0 0 0 0 0 0 0;0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1;0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2;...
      0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3;0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4;...
      0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5;0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6;...
      0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7;0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8];
Minfperp = [0.845 0.815 0.785 0.76 0.74 0.72 0.70 0.68;...
      0.83 0.80 0.77 0.75 0.725 0.705 0.685 0.67;...
      0.81 0.78 0.755 0.73 0.71 0.69 0.67 0.655;...
      0.79 0.765 0.74 0.72 0.70 0.68 0.66 0.645;...
      0.775 0.75 0.725 0.70 0.68 0.665 0.65 0.63;...
      0.76 0.73 0.71 0.69 0.67 0.65 0.64 0.62;...
      0.74 0.715 0.695 0.675 0.66 0.64 0.625 0.61;...
      0.725 0.70 0.68 0.66 0.645 0.63 0.615 0.6;...
      0.71 0.735 0.665 0.65 0.63 0.615 0.60 0.59];

% Data for interpolation of
% Figure 12: Incremental Drag Coefficient Due to Compressibility
% http://aero.stanford.edu/aa241/drag/computingcdc.html
Mratio = [0.75 0.85 0.95 1 1.025 1.05];
deltaCDc = [0 0.00025 0.00075 0.00125 0.002125 0.005];

% Data for profile drag increment of double-slotted flaps
% Figure 4.47
% Airplane Design part VI, page 84
flap_deflection = [0 30 60;0 30 60;0 30 60];
fc_ratio = [0.2 0.2 0.2;0.3 0.3 0.3;0.4 0.4 0.4];
delta_CDp_flap = [0 0.042 0.083;0 0.07 0.14;0 0.12 0.24];

% Data for flap k-factor
% Figure 4.53
% Airplane Design part VI, page 89
AR = [2 4 6 8;2 4 6 8;2 4 6 8];
out = [0.4 0.4 0.4 0.4;0.6 0.6 0.6 0.6;0.8 0.8 0.8 0.8];
kflap = [0.37 0.35 0.34 0.335;0.265 0.25 0.24 0.235;0.2 0.18 0.175 0.17];

% Data for nose gear drag
% Figure 4.58
% Airplane Design part VI, page 96
a_dt = [0.75 2.15 3.6;0.75 2.15 3.6;0.75 2.15 3.6];
e_dt = [1 1 1;2 2 2;3 3 3];
cd_nose = [0.27 0.45 0.38;0.47 0.66 0.63;0.7 0.92 0.9];

% Data for main gear drag
% Figure 4.59
% Airplane Design part VI, page 96
cd_main = polyfit([0.3 0.4 0.5 0.6 0.7 0.8],[0.23 0.45 0.68 1.1 1.8 2.7],3);

% Data for interpolation of Lifting Fan effects in horizontal tail
% (tiltwing case only)
% Figures 9-27
% Aerodynamics of V/STOL Flight, page 256
fan_CF_data = [0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0];
fan_CL_data = [0 0.5 1.0 1.5 1.9 2.3 2.9 3.3 3.8 4.2 4.75];
fan_CD_data = [0 0.05 0.2 0.32 0.42 0.52 0.6 0.68 0.73 0.85 0.9];

% Data for interpolation of 30 deg  rearward deflection of exit vanes
% for lifting fan, tiltwing case only
% Figures 9-29
% Aerodynamics of V/STOL Flight, page 258
vane_CL_data = [0 0.5 1.0 1.3 1.7 2.2 2.5 2.9 3.15 3.5 4.05];
vane_CD_data = [0 0.12 0.13 0.13 0.12 0.1 0 -0.1 -0.2 -0.3 -0.4];

% End of constants / lookup table data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fuselage Calculations

% Calculate fuselage maximum cross-sectional area
A_max = 3.14159*(h_max/2)*(w_max/2);

% Calculate fuselage Reynolds number
Re_fus = rho*V*L/mu;

% Calculate fineness ratio from geometrical data
Deff = ((4*A_max)/3.14159)^(1/2);
fineness = L/Deff;

% Find skin friction coefficient from lookup table
% including markdown for Mach number effects and
% markup for roughness
Cf_fus = ((100-f_markdown)/100)*((100+Cf_markup)/100)*...
   interp2(Reyn,trans,Cfinterp,Re_fus,xtc_fus);

% Calculate form factor from
% Curve fit to data on course notes page:
% http://aero.stanford.edu/aa241/drag/bodyformfactor.html
K_fus = 0.0092*fineness^2 - 0.1765*fineness + 1.9579;

% Calculate upsweep drag (based on fuselage cross sectional area)
CD_pi_upsweep = 0.75*hl;

% Calculate equivalent parasite area for fuselage
% including markup for interference effects
f_fus = ((100-Cf_inter)/100)*K_fus*Cf_fus*Swet_fuselage+...
   CD_pi_upsweep*A_max;

% Calculate fuselage parasite drag coefficient
CD_p_fus = f_fus/S_ref;

% End fuselage calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin wing calculations

% Calculate geometry for wing
AR_wing = b_wing^2/S_ref;
c_r_wing = 2*S_ref/(b_wing*(1+taper_wing));
c_t_wing = taper_wing*c_r_wing;
mac_wing = (2/3)*c_r_wing*((1+taper_wing+taper_wing^2)/(1+taper_wing));

% Calculate form factor for wing
% from  http://aero.stanford.edu/aa241/drag/lsformfactor.html
K_wing = 1 + (2.2*(tc_avg_wing)*cos(sweep_wing/57.3)^2)/...
   sqrt(1-0.5^2*cos(sweep_wing/57.3)^2)+...
   ((1.1)^2*cos(sweep_wing/57.3)^2*(tc_avg_wing)^2*(1+5*cos(sweep_wing/57.3)^2))/...
   (2*(1-0.5^2*cos(sweep_wing/57.3)^2));

% Calculate exposed area of wing based on fuselage diameter
r_fus = 0.5*w_max;
S_wing_exposed = S_ref*(1-(r_fus/(b_wing/2)))*(2-(1-taper_wing)*...
   ((r_fus/(b_wing/2))+1))/(1+taper_wing);

% Calculate the wetted area of the exposed wing
S_wet_wing = 2*S_wing_exposed*(1+0.2*tc_avg_wing);

% Calculate Reynolds number for wing
Re_wing = rho*V*mac_wing/mu;

% Find skin friction coefficient from lookup table
% including markdown for Mach number effects and
% markup for roughness
Cf_wing = ((100-f_markdown)/100)*((100+Cf_markup)/100)*...
   interp2(Reyn,trans,Cfinterp,Re_wing,xtc_wing);

% Calculate parasite area for wing
f_wing = K_wing*Cf_wing*S_wet_wing;

% Calculate f for control surface gap drag assuming the affected area is
% 30% of the wing reference area
f_gaps_wing = 0.0002*cos(sweep_wing/57.3)^2*(0.3*S_ref);

% Calculate parasite drag coefficient for wing
CD_p_wing = (f_wing + f_gaps_wing)/S_ref;

% end of wing calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Begin empennage calculations

% Calculate geometry for horizontal tail
AR_ht = b_ht^2/Sht_ref;
c_r_ht = 2*Sht_ref/(b_ht*(1+taper_ht));
c_t_ht = taper_ht*c_r_ht;
mac_ht = (2/3)*c_r_ht*((1+taper_ht+taper_ht^2)/(1+taper_ht));

% Calculate geometry for vertical tail(s)
AR_vt = b_vt^2/((1/N_vt)*Svt_ref);
c_r_vt = 2*(1/N_vt)*Svt_ref/(b_vt*(1+taper_vt));
c_t_vt = taper_vt*c_r_vt;
mac_vt = (2/3)*c_r_vt*((1+taper_vt+taper_vt^2)/(1+taper_vt));

% Calculate form factor for horizontal tail
% from  http://aero.stanford.edu/aa241/drag/lsformfactor.html
K_ht = 1 + (2.2*(tc_avg_ht)*cos(sweep_ht/57.3)^2)/...
   sqrt(1-0.5^2*cos(sweep_ht/57.3)^2)+...
   ((1.1)^2*cos(sweep_ht/57.3)^2*(tc_avg_ht)^2*(1+5*cos(sweep_ht/57.3)^2))/...
   (2*(1-0.5^2*cos(sweep_ht/57.3)^2));

% Calculate form factor for vertical tail
% from  http://aero.stanford.edu/aa241/drag/lsformfactor.html
K_vt = 1 + (2.2*(tc_avg_vt)*cos(sweep_vt/57.3)^2)/...
   sqrt(1-0.5^2*cos(sweep_vt/57.3)^2)+...
   ((1.1)^2*cos(sweep_vt/57.3)^2*(tc_avg_vt)^2*(1+5*cos(sweep_vt/57.3)^2))/...
   (2*(1-0.5^2*cos(sweep_vt/57.3)^2));

% Calculate the wetted area of the horizontal tail
% assuming exposed horizontal tail area equals the ref area (T-tail)
S_wet_ht = 2*Sht_ref*(1+0.2*tc_avg_ht);

% Calculate the wetted area of the vertical tail
S_wet_vt = 2*Svt_ref*(1+0.2*tc_avg_vt);

% Calculate Reynolds number for horizontal tail
Re_ht = rho*V*mac_ht/mu;

% Calculate Reynolds number for vertical tail
Re_vt = rho*V*mac_vt/mu;

% Find skin friction coefficient (for horizontal tail) from lookup table
% including markdown for Mach number effects and markup for roughness
Cf_ht = ((100-f_markdown)/100)*((100+Cf_markup)/100)*...
   interp2(Reyn,trans,Cfinterp,Re_ht,xtc_ht);

% Find skin friction coefficient (for vertical tail) from lookup table
% including markdown for Mach number effects and markup for roughness
Cf_vt = ((100-f_markdown)/100)*((100+Cf_markup)/100)*...
   interp2(Reyn,trans,Cfinterp,Re_vt,xtc_vt);

% Calculate parasite area for horizontal tail
f_ht = K_ht*Cf_ht*S_wet_ht;

% Calculate parasite area for vertical tail
f_vt = K_vt*Cf_vt*S_wet_vt;

% Calculate f for control surface gap drag assuming the affected area is
% 100% & 90% of the horizontal & vertical tail reference areas, respectively
f_gaps_ht = 0.0002*cos(sweep_ht/57.3)^2*(1.0*Sht_ref);
f_gaps_vt = 0.0002*cos(sweep_vt/57.3)^2*(0.9*Svt_ref);

% Calculate parasite drag coefficient for horizontal tail
CD_p_ht = (f_ht + f_gaps_ht)/S_ref;

% Calculate parasite drag coefficient for vertical tail
CD_p_vt = (f_vt + f_gaps_vt)/S_ref;

% end of empennage calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Nacelle and Pylon Drag (part of the parasite drag calculation)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Nozzle base drag
CD_nozzle_base = (0.5/12)*pi*D_exit*(0.2/S_ref);

% Calculate fineness ratio for the nacelles
fineness_n = (L_n + D_inlet)/(sqrt((4/3.14159)*(A_max_n - ...
   ((A_exit_n + A_inflow_n)/2))));

% Calculate nacelle form factor from
% Curve fit to data on course notes page:
% http://aero.stanford.edu/aa241/drag/bodyformfactor.html
K_nac = 0.0092*fineness_n^2 - 0.1765*fineness_n + 1.9579;

% Calculate the Reynolds Number for the nacelle
Re_nac = rho*V*L_n/mu;

% Find skin friction coefficient from lookup table
Cf_nac = ((100+Cf_markup)/100)*interp2(Reyn,trans,Cfinterp,Re_nac,xtc_nac);

% Nacelle wetted area calculation
D_eff_nac = (1/fineness_n)*L_n;
Swet_nac = 3.14159*D_eff_nac*L_n;

% Calculate equivalent parasite area for nacelle
f_nac = ((100-f_markdown)/100)*K_nac*Cf_nac*Swet_nac;

% Calculate nacelle parasited drag coefficient
CD_p_nac = f_nac/S_ref;

% Calculate pylon reference area from length and height
Spyl_ref = L_pyl*h_pyl;

% Calculate form factor for pylon
% from  http://aero.stanford.edu/aa241/drag/lsformfactor.html
K_pyl = 1 + (2.2*(tc_avg_pyl)*cos(sweep_pyl/57.3)^2)/...
   sqrt(1-0.5^2*cos(sweep_pyl/57.3)^2)+...
   ((1.1)^2*cos(sweep_pyl/57.3)^2*(tc_avg_pyl)^2*(1+5*cos(sweep_pyl/57.3)^2))/...
   (2*(1-0.5^2*cos(sweep_pyl/57.3)^2));

% Calculate the wetted area of the pylon
S_wet_pyl = 2*Spyl_ref*(1+0.2*tc_avg_pyl);

% Calculate the Reynolds Number for the pylon
Re_pyl = rho*V*L_pyl/mu;

% Find skin friction coefficient (for pylon) from lookup table
if Re_pyl == 0
   Cf_pyl = 0;
else
   Cf_pyl = ((100+Cf_markup)/100)*interp2(Reyn,trans,Cfinterp,Re_pyl,xtc_pyl);
end

% Calculate equivalent parasite area for pylon
f_pyl = ((100-f_markdown)/100)*K_pyl*Cf_pyl*S_wet_pyl;

% Calculate pylon parasite drag coefficient
CD_p_pyl = f_pyl/S_ref;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flap Drag
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Flap reference area
S_ref_flap = S_ref*(eta_o_flap - eta_i_flap)*(2-(1-taper_wing)*...
   (eta_i_flap + eta_o_flap))/(1+taper_wing);

% Calculate profile drag increment due to flaps
% Including interference effects
CD_p_flap = 0.4*cos(sweep_wing/57.3)*(S_ref_flap/S_ref)*...
   interp2(flap_deflection,fc_ratio,delta_CDp_flap,flap_angle,flap_to_chord);

if tiltwing == 1
   if flap_angle ~= 0
      % Calculate K (a factor) for the flaps (from lookup table)
      K_flap = interp2(AR,out,kflap,AR_wing,eta_o_flap);
      
      % Induced drag due to flaps (used for tiltwing only)
      CD_i_flap = K_flap^2*delta_CL_flap^2*cos(sweep_wing/57.3);
   else
      CD_i_flap = 0;
   end
else
   CD_i_flap = 0;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Slat Drag
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Area affected by slats
S_ref_slat = S_wing_exposed;   

CD_p_slat = 0.1*cos(sweep_wing/57.3)*CD_p_wing*(slat_ratio)*...
   (S_ref_slat/S_ref);

if gear == 1
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% Landing Gear Drag
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% nose calculations
	% Area of nose gear
   S_nose = b_t_nose*d_t_nose;
   
   a_dt_ratio = a_nose/d_t_nose;
      
   if a_dt_ratio > 3.6
      a_dt_ratio = 3.6;
   end
    
	% Calculate Cd for nose gear using lookup table
	CDp_nose = interp2(a_dt,e_dt,cd_nose,a_dt_ratio,e_nose/d_t_nose)*S_nose/S_ref;

	% main gear calculations
	% Area of main gear
	S_main = b_t_main*d_t_main;

	% Frontal area-to-total area ratio
	m_main = (S_main + (D_main*L_main))/(a_main*b_main);

	% Calculate Cd for main gear using curve fit
	CDp_main = polyval(cd_main,m_main)*S_main/S_ref;

	% Total landing gear drag
   CDp_gear = CDp_nose + CDp_main;
else
   CDp_gear = 0;
end
% End landing gear drag
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Caculate total parasite drag coefficient
CD_p = ((100-CDpmisc)/100)*(CD_p_fus + CD_p_wing + CD_p_ht + CD_p_vt ...
   + N_e*(CD_nozzle_base + CD_p_nac + CD_p_pyl) + CD_p_flap + CD_p_slat...
   + CDp_gear);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Drag Polar Calculations    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Assume taper ratio chosen such that u = 0.99 (per notes)
u = 0.99;

% Calculation lift-dependent drag factor for fuselage interference
s = 1 - 0.0407*(2*r_fus/b_wing) - 1.792*(2*r_fus/b_wing)^2;
sht = 1 - 0.0407*(2*r_fus/b_ht) - 1.792*(2*r_fus/b_ht)^2;

% Calculate e_inviscid assuming that the wing is twisted so that it achieves
% a nearly elliptical distribution of lift
e_inviscid = u*s;
eht_inviscid = u*sht;

% Assume 'K' factor for viscous induced drag is 0.38
% as discussed in AA241 course notes
K = 0.38;

% Calculate e
e = 1/(pi*AR_wing*(1/(pi*AR_wing*e_inviscid) + K*CD_p));
eht = 1/(pi*AR_ht*(1/(pi*AR_ht*eht_inviscid) + K*CD_p));

% Calculate total wetted and parasite area
Swet_total = Swet_fuselage+S_wet_vt+S_wet_ht+S_wet_wing...
   + Swet_nac + S_wet_pyl;
f_total = CD_p*S_ref;

% Calculate airplane drag polar
n=0;

if flap_angle == 0
   for i = 0:0.05:0.7
      
      n=n+1;
      CL(n) = i;
      
      Mperp = interp2(tcperp,CLover,Minfperp,...
         tc_avg_wing/cos(sweep_wing/57.3),CL(n)*cos(sweep_wing/57.3)^2);
      Mcc = Mperp/cos(sweep_wing/57.3);
      deltaCD_c = cos(sweep_wing)^3*interp1(Mratio,deltaCDc,M/Mcc,'spline');
      
      % Find CL for the horizontal tail and the wing/fuselage using trim code
      [CL_wf,CL_ht] = trim(b_wing,S_ref,taper_wing,sweep_wing_25c,...
         tc_avg_wing,x_wing,b_ht,Sht_ref,taper_ht,sweep_ht_25c,x_ht,x_cg,CL(n),M);
      
      % Induced Drag Terms
      CD_i = (CL_wf)^2/(pi*AR_wing*e) + (CL_ht)^2/(pi*AR_ht*eht)*(Sht_ref/S_ref);
      
      if gear == 1
         CD_i = CD_i + (p_nose*CL(n)*S_nose/S_ref) + (p_main*CL(n)*S_main/S_ref);
      end
      
      % Total Drag Coefficient for the wing/body/tail Cd = Cdo + Cdi + deltaCdc
      CD(n) = CD_p + CD_i + deltaCD_c;
   end
else
   if tiltwing == 1
      for i = 0:0.05:3.0
         n=n+1;
         CL(n) = i;
         
         % Find CL for the horizontal tail and the wing/fuselage using trim_tilt code
         [CL_wf,CL_ht] = trim_tilt(b_wing,S_ref,taper_wing,sweep_wing_25c,...
            tc_avg_wing,x_wing,b_ht,Sht_ref,taper_ht,sweep_ht_25c,x_ht,x_cg,CL(n),M,wing_angle);
         
         % Find the drag increment due to the lifting fan in the horizontal tail
         CF_fan = interp1(fan_CL_data,fan_CF_data,CL_ht,'spline');
         delta_CD_fan = interp1(fan_CF_data,fan_CD_data,CF_fan,'spline');
         CF_vane = interp1(vane_CL_data,fan_CF_data,CL_ht,'spline');
         delta_CD_vane = interp1(fan_CF_data,vane_CD_data,CF_vane,'spline');
         
         % Induced Drag Terms
         CD_i = (CL_wf)^2/(pi*AR_wing*e)+(CL_ht)^2/(pi*AR_ht*eht)*(Sht_ref/S_ref)+...
            CD_i_flap + delta_CD_fan*(Sht_ref/S_ref) + delta_CD_vane*(Sht_ref/S_ref);
         
         if gear == 1
            CD_i = CD_i + (p_nose*CL(n)*S_nose/S_ref) + (p_main*CL(n)*S_main/S_ref);
         end
         
         % Total Drag Coefficient for the wing/body/tail Cd = Cdo + Cdi + deltaCdc
         CD(n) = CD_p + CD_i;
      end
   else
      dummy = max(CL_data);
      for i = 0:0.1:dummy
         n=n+1;
         CL(n) = i;
         
         % Data for blown flaps
         blown_flaps = polyfit(CL_data,Cmu_data,2);
         
         % Calculate Cmu value for each CL
         Cmu = polyval(blown_flaps,CL(n));
         
         % Find CL for the horizontal tail and the wing/fuselage using trim code
         [CL_wf,CL_ht] = trim(b_wing,S_ref,taper_wing,sweep_wing_25c,...
            tc_avg_wing,x_wing,b_ht,Sht_ref,taper_ht,sweep_ht_25c,x_ht,x_cg,CL(n),M);
              
         % Correction factor from Induced Drag, Figure 7-13 in McCormick (VSTOL)
         k = 0.95 + (0.6/2)*CL_wf^2/(pi*AR_wing*e + 2*Cmu);
         
         % Induced Drag Terms
         CD_i = k*(CL_wf)^2/(pi*AR_wing*e + 2*Cmu) + (CL_ht)^2/(pi*AR_ht*eht)*(Sht_ref/S_ref);
         
         if gear == 1
            CD_i = CD_i + (p_nose*CL(n)*S_nose/S_ref) + (p_main*CL(n)*S_main/S_ref);
         end
         
         % Total Drag Coefficient for the wing/body/tail Cd = Cdo + Cdi
         CD(n) = CD_p + CD_i;
      end
   end
   
end
   
   
   if output == 1
      
      % Calculate L/D at each lift coefficient
      for i = 1:length(CL)
         LD(i) = CL(i)/CD(i);
      end
      
      if plot_number == 1
         fprintf('	TAKEOFF CONFIGURATION\n\n');
      elseif plot_number == 2
         fprintf('	CRUISE CONFIGURATION\n\n');
      elseif plot_number == 3
         fprintf('	LANDING CONFIGURATION\n\n');
      else end
      
      % Parasite Drag Component Breakdown
      fprintf('-------------------------------------------------------\n');
      fprintf('		PARASITE DRAG BREAKDOWN\n');
      fprintf('-------------------------------------------------------\n');
      fprintf('\n');
      fprintf('Fuselage:				%1.6f\n',CD_p_fus);
      fprintf('Wing:					%1.6f\n',CD_p_wing);
      fprintf('Horizontal Tail:			%1.6f\n',CD_p_ht);
      fprintf('Vertical Tail:				%1.6f\n',CD_p_vt);
      fprintf('Nacelle:				%1.6f\n',CD_p_nac);
      fprintf('Pylon:					%1.6f\n',CD_p_pyl);
      fprintf('Flaps:					%1.6f\n',CD_p_flap);
      fprintf('Slats:					%1.6f\n',CD_p_slat);
      fprintf('Landing Gear:				%1.6f\n',CDp_gear);
      fprintf('Miscellaneous Items:			%1.6f\n',((100-CDpmisc)/100)*CD_p + CD_nozzle_base);
      fprintf('TOTAL PARASITE DRAG COEFFICIENT:	%1.6f\n',CD_p);
      fprintf('\n')
      fprintf('-------------------------------------------------------\n');
      fprintf('\n');
      
      fprintf('-------------------------------------------------------\n');
      fprintf('\n');
      fprintf('Airplane Wetted Area:			%6.2f\n',Swet_total);
      fprintf('Airplane Equivalent Parasite Area:	%6.2f\n',f_total);
      fprintf('Airplane Span Efficiency Factor:	%1.6f\n',e);
      fprintf('\n')
      fprintf('-------------------------------------------------------\n');
      fprintf('\n');
      
      % Data for drag polar
      dragpolar = polyfit(CL,CD,3);
      
      if plot_number == 2
         % Calculate operating point at start of cruise
         CD_beg_cr = polyval(dragpolar,CL_beg_cr);
         
         % Calculate operating point at end of cruise
         CD_end_cr = polyval(dragpolar,CL_end_cr);
      end
   
      
      % Data for L/D vs. CL
      ldratio = polyfit(CL,LD,3);
      
      % Scale the L/D data
      CD_val_max = max(CD);
      LD_val_max = max(LD);
      scale = CD_val_max/LD_val_max;
      LD = scale*LD;
      
      if plot_number == 2
         % Calculate value for start of cruise
         LD_beg_cr = polyval(ldratio,CL_beg_cr);
         LD_beg_cr = scale*LD_beg_cr;
         
         % Calculate value for end of cruise
         LD_end_cr = polyval(ldratio,CL_end_cr);
         LD_end_cr = scale*LD_end_cr;
      end
      
      % Plotting the drag polar AND L/D ratio
      figure(plot_number)
      hold on
      plot(CL,CD,'k-',CL,LD,'k:');
      xlabel('Lift Coefficient');ylabel('Drag Coefficient');
      t1 = sprintf('Drag Polar and L/D vs. CL for M = %1.2f',M);
      title(t1);
      legend('Drag Polar','L/D vs. CL')
      
      if plot_number == 2
         plot(CL_beg_cr,CD_beg_cr,'X');
         plot(CL_end_cr,CD_end_cr,'O');
         plot(CL_beg_cr,LD_beg_cr,'X');
         plot(CL_end_cr,LD_end_cr,'O');
      end
      hold off
      
      figure(plot_number + 4)
      plot(CL,LD/scale,'k:');
      title('Lift to Drag Ratio vs. Lift Coefficient')
      xlabel('Lift Coefficient');ylabel('Lift-to-Drag Ratio')
      
   end
   