%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Patrick LeGresley
% Blown Flaps / High Lift
% Last modified: May 1, 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
%  flap_to_chord = flap to chord ratio
%  eta_i_flap = flap inboard station
%  eta_o_flap = flap outboard station
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs:
%  CL_max = clean wing maximum lift coefficient at zero Mach number
%  flap_setting, C_mu_data, CL_data = data for interpolation of lift coefficient
%     for blown flaps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [CL_max, C_mu_data, CL_data] = highlift(tiltwing,b_wing,S_ref,taper_wing,...
	sweep_wing_25c,tc_avg_wing,flap_to_chord,eta_i_flap,eta_o_flap,slat_angle,flap_angle)


if tiltwing ~=1
   
   % Misc. stuff
   flap_angle = flap_angle/57.3;
   kappa = 0.97;
   beta = 1;
   alpha = 5/57.3;
   alpha_max_Cmu0 = 17;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Find clean wing maximum lift coefficient -- CL_max
   
   [rho, Temp, a, mu] = atm(0);
   sigma = 1;
   
   % calculate 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;
   
   % find chord & t/c for outer panel (at 75% semi-span)
   c_outerpanel = c_r_wing - 0.75*(c_r_wing-c_t_wing);
   tc_outerpanel = 0.9*tc_avg_wing;
   
   % find local Reynolds number for the outer panel at landing stall speed
   V_S_L = 70;
   Re_outerpanel = sigma*rho*V_S_L*c_outerpanel/mu;
   
   % airfoil data from Theory of Wing Section by Abbott and von Doenhoff
   clmax_ref = 1.65;
   Re_ref = 9*10^6;
   
   % section clmax corrected for Reynolds number effects
   clmax = clmax_ref*(Re_outerpanel/Re_ref)^0.1;
   
   % lookup table for Figure 3: Effect of taper and sweep on wing/outer panel Clmax
   % for positive sweep angles only
   taper = [0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;];
   sweep = [0 0 0 0 0;10 10 10 10 10;20 20 20 20 20;30 30 30 30 30;40 40 40 40 40];
   C1 = [0.885 0.93 0.925 0.90 0.865; 0.865 0.92 0.925 0.915 0.88;0.83 0.905 0.92 0.915 0.90;...
         0.815 0.885 0.915 0.915 0.91;0.79 0.875 0.915 0.92 0.915];
   
   % find wing maximum lift coefficient (3-D) from figure 3
   CLmax_wing = interp2(taper,sweep,C1,taper_wing,sweep_wing_25c)*clmax;
   
   % wing CLmax based on FAA definition
   CL_max = 1.11*CLmax_wing;
   
   % end of calculations for clean wing
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Calculations for blown flaps
   
   
   
   % find flap affected area from inboard and outboard stations
   S_flap = S_ref*(eta_o_flap-eta_i_flap)*(2-(1-taper_wing)*...
      (eta_i_flap+1))/(1+taper_wing);
   
   
   % 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));
   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);
   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));
   
   dCl_da_inf = C_L_alpha_wing;
   dCl_da_inf_Cmu0 = C_L_alpha_wing;
   
   
   
   j = 0;
   
   for C_mu = 0:0.1:16
      
      j = j + 1;
      
      % Eqn. 7-3 (for c_f/c = 0)
      dCl_ddelta = sqrt(4*pi*C_mu*(1+0.151*C_mu^0.5+0.139*C_mu));
      
      % Eqn. 7-4 (for c_f/c = 0)
      dCl_da = 2*pi*(1+0.151*C_mu^0.5+0.219*C_mu);
      
      % Figure 7-2, p. 195
      cf_c_data = polyfit([0 0.2 0.3 0.5 1],[0 3.5 4 5.3 6.5],3);
      
      dCl_ddelta = dCl_ddelta*polyval(cf_c_data,flap_to_chord);
      dCl_da = dCl_da*polyval(cf_c_data,flap_to_chord);
      
      % Eqn. 7-6, p. 196
      F = (AR_wing + 2*C_mu/pi)/(AR_wing + 2 + 0.604*C_mu^0.5 + 0.876*C_mu);
      
      % Eqn. 7-22
      lambda_flap = S_flap/S_ref;
      b_flap = (eta_o_flap - eta_i_flap)*b_wing;
      
      % Eqn. 7-23
      v = (S_flap*dCl_da_inf + (S_ref - S_flap)*dCl_da_inf_Cmu0)/(S_ref*dCl_da_inf);
      
      % Calculate C_mu based on S
      C_mu = lambda_flap*C_mu;
      
      % Eqn. 7-24
      C_L = F*(1 + tc_avg_wing)*(lambda_flap*dCl_ddelta*flap_angle + v*dCl_da*alpha) - tc_avg_wing*C_mu*(flap_angle+alpha);
      
      C_mu_data(j) = C_mu*b_flap;
      
      C_L_aero = F*(1 + tc_avg_wing)*(v*dCl_da*alpha);
      
      if C_L_aero > 1*AR_wing
         CL_data(j) = -1;
      else
         CL_data(j) = C_L;
      end
      
      
      % Figure 7-8, p. 202
      %alphamax_data = polyfit([0 1 2 3 4 5 6 7],[0 -3 0 2 3.3 4 5 6],3);
      
      %alpha_max = alpha_max_Cmu0 + polyval(alphamax_data,C_mu);
      
      
   end
   
   
   
   
else
   
   % Misc. stuff
   flap_angle = flap_angle/57.3;
   kappa = 0.97;
   beta = 1;
   alpha = 5/57.3;
   alpha_max_Cmu0 = 17;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Find clean wing maximum lift coefficient -- CL_max
   
   [rho, Temp, a, mu] = atm(0);
   sigma = 1;
   
   % calculate 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;
   
   % find chord & t/c for outer panel (at 75% semi-span)
   c_outerpanel = c_r_wing - 0.75*(c_r_wing-c_t_wing);
   tc_outerpanel = 0.9*tc_avg_wing;
   
   % find local Reynolds number for the outer panel at landing stall speed
   V_S_L = 70;
   Re_outerpanel = sigma*rho*V_S_L*c_outerpanel/mu;
   
   % airfoil data from Theory of Wing Section by Abbott and von Doenhoff
   clmax_ref = 1.65;
   Re_ref = 9*10^6;
   
   % section clmax corrected for Reynolds number effects
   clmax = clmax_ref*(Re_outerpanel/Re_ref)^0.1;
   
   % lookup table for Figure 3: Effect of taper and sweep on wing/outer panel Clmax
   % for positive sweep angles only
   taper = [0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;0.2 0.4 0.6 0.8 1;];
   sweep = [0 0 0 0 0;10 10 10 10 10;20 20 20 20 20;30 30 30 30 30;40 40 40 40 40];
   C1 = [0.885 0.93 0.925 0.90 0.865; 0.865 0.92 0.925 0.915 0.88;0.83 0.905 0.92 0.915 0.90;...
         0.815 0.885 0.915 0.915 0.91;0.79 0.875 0.915 0.92 0.915];
   
   % find wing maximum lift coefficient (3-D) from figure 3
   CLmax_wing = interp2(taper,sweep,C1,taper_wing,sweep_wing_25c)*clmax;
   
   % end of calculations for clean wing
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Calculations for Slats
   
   % curve fit for Figure 6: clmax increment due to slats
   C2 = polyfit([0 5 10 15 20 25],[0 0.15 0.4 0.55 0.75 1.25],3);
   
   % curve fit for Figure 7: effect of sweep on slat CLmax increment
   C3 = polyfit([0 10 20 30 40 50 60],[1 0.95 0.85 0.725 0.55 0.35 0.1],3);
   
   % calculate increment in wing lift coefficient (3-D) due to flaps
   dCLmax_slats = polyval(C2,slat_angle)*polyval(C3,sweep_wing_25c);
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Calculations for Flaps
   
   % curve fit for figure 8: section clmax increment due to flaps for
   % double slotted flaps, 25% flap, 50 deg defl
   %
   % multiply by 0.93 for single slotted flaps
   % multiply by 1.08 for triple slotted flaps
   C4 = polyfit([0 5 10 15 20],[1 1.12 1.4 1.97 2.27],3);
   % assume single slotted flaps
   dClmax_ref = 1*polyval(C4,tc_outerpanel*100);
   
   % calculate flap chord correction factor from Figure 9
   K1 = (1/25)*flap_to_chord;
   
   % curve fit for Figure 10: flap motion correction factor
   C5 = polyfit([0 10 20 30 40 50],[0 0.26 0.5 0.71 0.87 1],3);
   K2 = polyval(C5,flap_angle*57.3);
   
   % correction factor for sweep
   K_sweep = (1-0.08*cos(sweep_wing_25c/57.3)^2)*cos(sweep_wing_25c/57.3)^0.75;
   
   % increment in section lift coefficient (2-D) due to flaps
   dClmax_flaps = K1*K2*dClmax_ref;
   
   % find flap affected area from inboard and outboard stations
   S_flaps_eff = S_ref*(eta_o_flap-eta_i_flap)*(2-(1-taper_wing)*...
      (eta_i_flap+1))/(1+taper_wing);
   
   % increment in wing lift coefficient (3-D) due to flaps
   dCLmax_flaps = (S_flaps_eff)*dClmax_flaps*K_sweep;
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Calculate total CLmax with flaps and slats
   CL_max = CLmax_wing + dCLmax_flaps + dCLmax_slats;
   
   % don't calculate these for the tiltwing
   C_mu_data = 0;
   CL_data = 0;
   
end

