%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Patrick LeGresley
% Trim
% Last modified: April 28, 1999
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Known problems with this code:
%
%  none :)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Inputs:
%  b_wing = wing span (ft)
%  S_ref = wing reference area (ft^2)
%  taper_wing = wing taper ratio
%  sweep_wing_25c = wing quarter chord sweep angle (deg)
%  tc_avg_wing = wing average thickness to chord ratio
%  x_wing = x-position of wing root leading edge (ft)
%  b_ht = horizontal tail span (ft)
%  Sht_ref = horizontal tail reference area (ft^2)
%  taper_ht = horizontal tail taper ratio
%  sweep_ht_25c = horizontal tail quarter chord sweep angle (deg)
%  x_ht = x-position of horizontal tail root leading edge (ft)
%  x_cg = x-position of center of gravity (ft)
%  CL_1 = steady state airplane lift coefficient
%  M = Mach number
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs:
%  CL_wf = lift coefficient of the wing-fuselage
%  CL_ht = horizontal tail lift coefficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [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_1,M,wing_angle)


% misc. stuff
alpha_o_l = -2;
epsilon_t = -5;
i_w = 3;
c_m_o = -0.04;
n_ht = 0.95;
kappa = 0.97;
beta = sqrt(1-M^2);
L_fus = 60;
w_fus = 13;
h_fus = 9.5;
A_fus = 3.14159*(w_fus/2)*(h_fus/2);
D_eff_fus = 2*sqrt(A_fus/3.14159);

% AA241 Course Notes
Kf_data = polyfit([.1 .2 .3 .4 .5 .6 .7],[.115 .172 .344 .487 .688 .888 1.146],3);

% Figure 8.99, A/P Design Part VI, page 306
CmoM = polyfit([0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8],[1 1 1 1 1.04 1.07 1.13 1.2 1.28],3);

% Figure 8.42, A/P Design Part VI, page 247
Mcos=[0 0.3 0.4 0.5 0.6 0.7;0 0.3 0.4 0.5 0.6 0.7;0 0.3 0.4 0.5 0.6 0.7;...
   0 0.3 0.4 0.5 0.6 0.7;0 0.3 0.4 0.5 0.6 0.7];
tc=[0.08 0.08 0.08 0.08 0.08 0.08;0.1 0.1 0.1 0.1 0.1 0.1;0.12 0.12 0.12 0.12 0.12 0.12;...
      0.14 0.14 0.14 0.14 0.14 0.14;0.16 0.16 0.16 0.16 0.16 0.16];
alpha_ratio = [1 1 1 1 1 1;1 1 1 1 1 1;1 1 1 1 1 0.6;1 1 1 1 0.8 -0.05;1 1 1 0.85 0.4 -1];

% Figure 8.41, A/P Design Part VI, page 247
% curve fit to data for taper=1 and 0 sweep 
delta_alpha = polyfit([1.5 3.5 6 10],[-.408 -.418 -.4275 -.4375],3);

% Figure 8.98, A/P Design Part VI, page 304
% data for taper=1
sweep=[0 20 40 60;0 20 40 60;0 20 40 60;0 20 40 60;0 20 40 60];
AR=[1.5 1.5 1.5 1.5;3.5 3.5 3.5 3.5;6 6 6 6;8 8 8 8;10 10 10 10];
dCmo = [0 -0.0005 -0.0001 -0.0015; 0 -0.0015 -0.003 -0.006;0 -0.0035 -0.0075 -0.0135;...
      0 -0.006 -0.012 -0.2;0 -0.008 -0.017 -0.026];

% Figure 8.111, A/P Design Part VI, page 319
massfactor = polyfit([4 8 12 16 20],[0.77 0.915 0.96 0.97 0.98],3);


% find wing geometry
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));
y_mac_wing = b_wing*(1+2*taper_wing)/6/(1+taper_wing);
sweep_wing_0c = atan(tan(sweep_wing_25c/57.3) + (1-taper_wing)/AR_wing/(1+taper_wing));
sweep_wing_50c = tan(sweep_wing_0c) - 4*(0.5)*(1-taper_wing)/AR_wing/(1+taper_wing);
n_mac_wing = tan(sweep_wing_0c)*y_mac_wing;
x_ac_wing = x_wing + n_mac_wing + 0.25*mac_wing;
C_L_alpha_wing = 2*pi*AR_wing/(2+sqrt(AR_wing^2*beta^2/kappa^2*(1+tan(sweep_wing_50c/57.3)^2/beta^2)+4));

% find fuselage induced shift in aerodynamic center
Kf = polyval(Kf_data,(x_wing+0.25*c_r_wing)/L_fus);
dCm_CL_fus = Kf*w_fus^2*L_fus/S_ref/mac_wing/C_L_alpha_wing;
x_ac_wf = x_ac_wing - mac_wing*dCm_CL_fus;

% find horizontal tail geometry
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));
y_mac_ht = b_ht*(1+2*taper_ht)/6/(1+taper_ht);
sweep_ht_0c = atan(tan(sweep_ht_25c/57.3) + (1-taper_ht)/AR_ht/(1+taper_ht));
sweep_ht_50c = tan(sweep_ht_0c) - 4*(0.5)*(1-taper_ht)/AR_ht/(1+taper_ht);
n_mac_ht = tan(sweep_ht_0c)*y_mac_ht;
x_ac_ht = x_ht + n_mac_ht + 0.25*mac_ht;
C_L_alpha_ht_o = 2*pi*AR_ht/(2+sqrt(AR_ht^2*beta^2/kappa^2*(1+tan(sweep_ht_50c/57.3)^2/beta^2)+4));
deps_dalpha = 2*C_L_alpha_wing/(3.14159*AR_wing);
C_L_alpha_ht = C_L_alpha_ht_o*(1-deps_dalpha)*n_ht;

% find Cmobar
alpha_o_L_w = (alpha_o_l + polyval(delta_alpha,AR_wing)*epsilon_t)*...
   interp2(Mcos,tc,alpha_ratio,M*cos(sweep_wing_25c/57.3),tc_avg_wing);
C_L_o_wf = ((i_w - alpha_o_L_w)/57.3)*C_L_alpha_wing;
C_L_o = C_L_o_wf + 0;
C_m_o_w = ((AR_wing*cos(sweep_wing_25c)^2)/(AR_wing+2*cos(sweep_wing_25c)))*c_m_o +...
   interp2(sweep,AR,dCmo,sweep_wing_25c,AR_wing)*epsilon_t;
fus_var = 0;
for i=1:13
   fus_var = fus_var + w_fus^2*(i_w + alpha_o_L_w + 0)/57.3*(L_fus/13);
end
C_m_o_f = polyval(massfactor,L_fus/D_eff_fus)/(36.5*S_ref*mac_wing)*fus_var;
C_m_o_wf = (C_m_o_w + C_m_o_f)*polyval(CmoM,M);
C_m_o = C_m_o_wf;
SM = -(x_cg-x_ac_wing)/mac_wing + (x_ac_ht-x_cg)*Sht_ref*C_L_alpha_ht/(mac_wing*S_ref*C_L_alpha_wing) - dCm_CL_fus;
Cmobar = C_m_o + SM*C_L_o;

% find CL_ht & CL_wf
C_1 = (x_cg - x_ac_wf)/mac_wing;
C_2 = (x_cg - x_ac_ht)/mac_wing*n_ht*Sht_ref/S_ref;
CL_wf = (CL_1*C_2 + Cmobar*(Sht_ref/S_ref)*n_ht)/(C_2 - C_1*(Sht_ref/S_ref)*n_ht);
CL_ht = -(Cmobar + C_1*CL_wf)/C_2;
