%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Patrick LeGresley
% Performance
% Last modified: May 1, 1999
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clc;
%clear all;
%close all;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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)
%  h = cruise altitude (ft)
%  V = cruise speed (KTAS)
%  T_to_uninst = takeoff uninstalled thrust (lbs)
%  N_e = number of engines
%  CL_max = clean wing maximum lift coefficient at zero Mach number
%  Cmu_data_to = momentum coefficient data for takeoff configuration
%  CL_data_to = lift coefficient data for takeoff configuration
%  Cmu_data_land = momentum coefficient data for landing configuration
%  CL_data_land = lift coefficient data for landing configuration
%  CD_to = drag coefficient data for takeoff configuration
%  CL_to = lift coefficient data for takeoff configuration
%  CD_cr = drag coefficient data for cruise configuration
%  CL_cr = lift coefficient data for cruise configuration
%  CD_land = drag coefficient data for landing configuration
%  CL_land = lift coefficient data for landing configuration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Outputs:
%  W_f_used = weight of fuel used during mission (lbs)
%  R_bl = block range (nm)
%  R_cruise = maximum cruise range (nm)
%  t_cl = time to climb to cruise altitude (hours)
%  T_to = takeoff installed thrust (lbs)
%  S_g_to = takeoff ground roll (ft)
%  S_g_land = landing ground roll (ft)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%function [W_f_used,R_bl,R_cruise,t_cl,T_to,S_g_to,S_g_land] = performance(b_wing,S_ref,taper_wing,...
%   sweep_wing_25c,h,V,T_to_uninst,N_e,W_payload,CL_max,Cmu_data_to,CL_data_to,Cmu_data_land,...
%   CL_data_land,W_to,W_e,CD_to,CL_to,CD_cr,CL_cr,CD_land,CL_land)


function [W_f_used,R_bl,R_cruise,t_cl,T_to,S_g_to,S_g_land] = performance(tiltwing,b_wing,S_ref,...
   taper_wing,sweep_wing_25c,flap_angle_to,flap_angle_land,wing_angle_to,wing_angle_land,h,V,...
   T_to_uninst,N_e,W_payload,CL_max,CL_max_to,Cmu_data_to,CL_data_to,CL_max_land,Cmu_data_land,...
   CL_data_land,W_to,W_e,CD_to,CL_to,CD_cr,CL_cr,CD_land,CL_land)



% Misc. calculations
W_tfo = 0.005*W_to;
W_crew = 400;
W_fuel = W_to - W_e - W_tfo - W_crew - W_payload;
DP_to = polyfit(CL_to,CD_to,3);
DP_cr = polyfit(CL_cr,CD_cr,3);
DP_land = polyfit(CL_land,CD_land,3);


if tiltwing ~= 1
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Climb & Cruise Range
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   % create installed data
   % Note: T_to_un is the TOTAL airplane uninstalled thrust
   [T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
         wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
         wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix]=installed(W_to,T_to_uninst);
   
   W = W_to;
   
   % start, warmup, taxi (total of 10 minutes)
   
   T = T_idle_in(1,1);
   
   fuelrate=interpolation(4,0,0.3,T,...
      T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
      wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
      wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
   
   W = W - N_e*fuelrate*(10/60);
   
   
   
   % takeoff (assume 3 minutes at TO power)
   
   T = T_to_in(1,1);
   T_to = N_e*T;
   
   fuelrate=interpolation(1,0,0.3,T,...
      T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
      wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
      wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
   
   W = W - N_e*fuelrate*(3/60);
   
   
   R = 0;
   t_cl = 0;
   
   % climb
   
   CD_p = polyval(DP_cr,0);
   B = 0.5^2/(polyval(DP_cr,0.5)-CD_p);
   
   T_climb = T_to_in(1,:);
   h_climb = h_matrix(1,:);
   Tav_climb = polyfit(h_climb,T_climb,3);
   
   for alt=500:1000:h-500
      
      [rho, Temp, a, mu] = atm(alt);
      
      T_av = polyval(Tav_climb,alt);
      
      C_L = -(N_e*T_av/W)*B/2 + sqrt(((N_e*T_av/W)*B/2)^2 + 3*CD_p*B);
      V_climb = sqrt(2*W/(rho*S_ref*C_L));
      M = V_climb/a;
      
      qbar = 0.5*rho*V_climb^2;
      C_D = polyval(DP_cr,C_L);
      D = C_D*qbar*S_ref;
      
      % assume climb at constant EAS
      % RC in ft/sec
      RC = (N_e*T_av - D)*V_climb/W/(1+0.567*M^2);
      
      % dt in hours
      dt = (1000/RC)/3600;
      
      % dR in nm
      dR = 0.5925*V_climb*dt;
      
      % dW in lbs
      [fuelrate,T_av]=interpolation(1,alt,0.3,T_av,...
         T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
         wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
         wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
      
      dW = N_e*fuelrate*dt;
      
      W = W - dW;
      R = R + dR;
      t_cl = t_cl + dt;
   end
   
   R_climb = R;
   W_begcr = W;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   Wfuel_dlt = 0;
   
   % descent
   
   RD = -2500;	% fpm
   dt = (1000/(-RD))/60;	% hours
   
   for alt = h-500:-1000:500
      fuelrate=interpolation(1,0,0.3,T,...
         T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
         wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
         wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
      
      dW = fuelrate*dt;	
      Wfuel_dlt = Wfuel_dlt + dW;
   end
   
   
   
   % land (4 minutes)
   
   T = T_to_in(1,1);
   
   fuelrate=interpolation(1,0,0.3,T,...
      T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
      wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
      wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
   
   Wfuel_dlt = Wfuel_dlt + fuelrate*(4/60);
   
   
   
   % taxi (5 minutes)
   
   T = T_idle_in(1,1);
   
   fuelrate=interpolation(4,0,0.3,T,...
      T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
      wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
      wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
   
   Wfuel_dlt = Wfuel_dlt + fuelrate*(5/60);
   
   
   
   W_endcr = W_e + W_tfo + W_crew + W_payload + 0.1*W_fuel + Wfuel_dlt;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   % cruise
   
   [rho, Temp, a, mu] = atm(h);
   M = 1.6878*V/a;
   qbar = 0.5*rho*(1.6878*V)^2;
   
   for w=W_begcr:-10:W_endcr
      C_L = w/(qbar*S_ref);
      C_D = polyval(DP_cr,C_L);
      D = C_D*qbar*S_ref;
      
      fuelrate=interpolation(2,h,M,D/N_e,...
         T_to_in,T_max_in,T_1_in,T_2_in,T_3_in,T_4_in,T_5_in,T_6_in,T_idle_in,...
         wf_matrix_to,wf_matrix_max,wf_matrix_1,wf_matrix_2,wf_matrix_3,wf_matrix_4,wf_matrix_5,...
         wf_matrix_6,wf_matrix_idle,h_matrix,M_matrix);
      
      
      SR = V/(N_e*fuelrate);
      dR = SR*10;
      R = R + dR;
      
   end
   
   
   R_bl = R;
   R_cruise = R - R_climb;
   W_f_used = W_fuel/1.1;
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Takeoff
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   CL_to = .6;                % Lift Coefficient During TO Ground Roll
   CD_to = CL_to/8;           % Drag Coefficient During TO Ground Roll
   
   g = 32.2;                  % Gravitational Acceleration
   rho = 0.00226;             % Air Density
   mu_fric = .025;            % Ground Friction Coefficient
   mu_brake = 0.3;            % Braking Friction Coefficient
   stepsize = 1;              % Set stepsize for numerical integration
   
   a = sqrt(1.4*1716*(59+459.7+27));
   CL_max_to = max(CL_data_to);
   
   
   Vs_to = sqrt(2*W_to/(rho*S_ref*CL_max_to)); % Stall speed in TO configuration
   V_lo=1.10*Vs_to;                        % Calculate liftoff speed
   
   
   v = [0:stepsize:V_lo];
   sg_to(1)=0;
   
   for i=2:length(v)
      
      if (v(i)>1.1*Vs_to)                  % Values for after plane has been rotated
         CL = CL_max_to/1.21;              % Lift coefficient after plane has rotated
         CD = CL/6;                        % Drag Coefficient after plane has rotated
      end
      
      drag = .5*rho*v(i)^2*S_ref*CD_to;    % Calculate drag force
      lift = .5*rho*v(i)^2*S_ref*CL_to;    % Calculate lift force
      
      T = T_to_uninst*(-(v(i)/a)+0.97);           % Variation of T w/ M from curve fit w/ AA241 course notes
      
      % Calculate incremental takeoff distances and sum to numerically integrate
      sg_to(i)=sg_to(i-1)+stepsize*((v(i))/(g/W_to*(T-drag-mu_fric*(W_to-lift-T))));
      
   end
   
   S_g_to = sg_to(length(v));
   
   
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Landing
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   stepsize = 1;
   
   W_L = 0.85*W_to;           % Landing Weight
   CL_max_L = 8;              % Maximum Lift Coefficient (Landing Configuration)
   
   LD_app = 6;                % Lift to Drag Ratio in Approach Configuration
   CL_L = 0;                  % Lift Coefficient During Landing Ground Roll
   CD_L = 0.5;               % Drag Coefficient During Landing Ground Roll
   
   
   
   
   Vs_L = sqrt(2*W_L/(rho*S_ref*CL_max_L));    % Stall speed in TO configuration
   
   V_td=1.15*Vs_L;                         % Calculate touchdown speed
   
   v=[V_td:-stepsize:0];                   % Numerically integrate from touchdown speed to zero
   
   sg_l(1)=V_td;                           % Assume 1 second for rotation during landing
   
   for i=2:length(v)
      T=0;                                 % Assume there is no thrust during landing ground roll
      drag=.5*rho*v(i)^2*S_ref*CD_L;       % Calculate drag during landing
      % lift=.5*rho*v(i)^2*S_ref*CL_L;     % Calculate lift during landing
      lift = 0;                            % Assume there is no lift during landing ground roll
      
      % Calculate incremental landing distances and sum to numerically integrate
      sg_l(i)=sg_l(i-1)-stepsize*(v(i)/(g*((T/W_L)+(lift*mu_brake-drag)/W_L+0.08*(mu_brake-mu_fric)-mu_brake)));
   end
   
   S_g_land = sg_l(length(v));           % Output landing ground roll
   
   
   
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
else
   
   P_to_uninst = T_to_uninst;
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Climb & Cruise Range
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   % create installed data
   % Note: %  P_to_un is the uninstalled power for one engine

   [wf_to_matrix,wf_c_matrix,wf_cr_matrix_A,wf_cr_matrix_B,wf_cr_matrix_C,wf_cr_matrix_D,...

         h_to_matrix,h_c_matrix,h_cr_matrix,v_to_matrix,v_c_matrix,v_cr_matrix,...

         P_to_in,P_c_in,P_cr_in_A,P_cr_in_B,P_cr_in_C,P_cr_in_D]=prop_inst_funct(P_to_uninst/N_e,W_to);

   
   W = W_to;
   
   % start, warmup, taxi (total of 10 minutes)
   

   fuelrate = 646*(P_to_uninst/N_e)/2500;

   
   W = W - N_e*fuelrate*(10/60);
   
   
   
   % takeoff (assume 5 minutes at TO power)
   

   fuelrate = 1437*(P_to_uninst/N_e)/2500;

   
   W = W - N_e*fuelrate*(5/60);

   

   
   R = 0;
   t_cl = 0;
   
   % climb
   
   %CD_p = polyval(DP_cr,0);
   %B = 0.5^2/(polyval(DP_cr,0.5)-CD_p);
   
   P_climb = P_c_in(3,:);
   h_climb = h_c_matrix(1,:);
   Pav_climb = polyfit(h_climb,P_climb,3);
   

   
   for alt=500:1000:h-500
      
      [rho, Temp, a, mu] = atm(alt);
      
      P_av = polyval(Pav_climb,alt);
      V_climb = 250;

      T_av = P_av*550/(1.6878*V_climb);

      
      qbar = 0.5*rho*(1.6878*V_climb)^2;

      C_L = W/qbar/S_ref;
      C_D = polyval(DP_cr,C_L);
      D = C_D*qbar*S_ref;
      
      % RC in ft/sec
      RC = (N_e*T_av - D)*(1.6878*V_climb)/W; 
      
      % dt in hours
      dt = (1000/RC)/3600;
      
      % dR in nm
      dR = V_climb*dt;
      
      % dW in lbs
      fuelrate=prop_fuel_flow(2,alt,V_climb,P_av,wf_to_matrix,wf_c_matrix,wf_cr_matrix_A,...

         wf_cr_matrix_B,wf_cr_matrix_C,wf_cr_matrix_D,h_to_matrix,h_c_matrix,h_cr_matrix,...

         v_to_matrix,v_c_matrix,v_cr_matrix,P_to_in,P_c_in,P_cr_in_A,P_cr_in_B,P_cr_in_C,P_cr_in_D);
      
      dW = N_e*fuelrate*dt;
      
      W = W - dW;
      R = R + dR;
      t_cl = t_cl + dt;
   end
   
   R_climb = R;
   W_begcr = W;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   Wfuel_dlt = 0;
   
   % descent
   
   RD = -2500;	% fpm
   dt = (1000/(-RD))/60;	% hours
   
   for alt = h-500:-1000:500

      

      P = 0.65*W/8*(250*0.5925)/550;

      

      if alt > 20000
         fuelrate=prop_fuel_flow(3,alt,250,P,wf_to_matrix,wf_c_matrix,wf_cr_matrix_A,...

            wf_cr_matrix_B,wf_cr_matrix_C,wf_cr_matrix_D,h_to_matrix,h_c_matrix,h_cr_matrix,...

            v_to_matrix,v_c_matrix,v_cr_matrix,P_to_in,P_c_in,P_cr_in_A,P_cr_in_B,P_cr_in_C,P_cr_in_D);

      end
      
      dW = fuelrate*dt;	
      Wfuel_dlt = Wfuel_dlt + dW;
   end
   
   
   
   % land (8 minutes)
   
   fuelrate = 1437*(P_to_uninst/N_e)/2500;
   
   Wfuel_dlt = Wfuel_dlt + fuelrate*(8/60);
   
   
   
   % taxi (5 minutes)
   
   fuelrate = 646*(P_to_uninst/N_e)/2500;
   
   Wfuel_dlt = Wfuel_dlt + fuelrate*(5/60);
   
   W_endcr = W_e + W_tfo + W_crew + W_payload + 0.1*W_fuel + Wfuel_dlt;
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   % cruise
   
   [rho, Temp, a, mu] = atm(h);
   %M = 1.6878*V/a;
   qbar = 0.5*rho*(1.6878*V)^2;
   
   for w=W_begcr:-25:W_endcr
      C_L = w/(qbar*S_ref);
      C_D = polyval(DP_cr,C_L);
      D = C_D*qbar*S_ref;

      

      P = D*1.6878*V/550;
      
      fuelrate=prop_fuel_flow(3,h,V,P/4,wf_to_matrix,wf_c_matrix,wf_cr_matrix_A,...

         wf_cr_matrix_B,wf_cr_matrix_C,wf_cr_matrix_D,h_to_matrix,h_c_matrix,h_cr_matrix,...

         v_to_matrix,v_c_matrix,v_cr_matrix,P_to_in,P_c_in,P_cr_in_A,P_cr_in_B,P_cr_in_C,P_cr_in_D);
      
      
      SR = V/(N_e*fuelrate);
      dR = SR*25;
      R = R + dR;
      
   end
   
   
   R_bl = R;
   R_cruise = R - R_climb;
   W_f_used = W_fuel/1.1;
   

   

   

   %

   % Lookup table for variation of thrust with forward speed and angle of rotor

   %

   % Based on data from Aerodynamics of V/STOL Flight by McCormick

   % Figure 8-5: Variation of Available Thrust with Speed for Constant Ideal Power, page 217

   %

   % x = V/w_o,  y = theta = rotor angle (alpha in book),  z = T/T_o

   %

   x=[0 1 2 3 4 5; 0 1 2 3 4 5; 0 1 2 3 4 5; 0 1 2 3 4 5;0 1 2 3 4 5; 0 1 2 3 4 5; 0 1 2 3 4 5;...

         0 1 2 3 4 5];

   y=[0 0 0 0 0 0; 30 30 30 30 30 30; 40 40 40 40 40 40; 50 50 50 50 50 50; 60 60 60 60 60 60;...

         70 70 70 70 70 70; 80 80 80 80 80 80; 90 90 90 90 90 90];

   z=[0.995 0.65 0.45 0.32 0.28 0.20;0.996 0.75 0.50 0.40 0.30 0.25;0.997 0.78 0.57 0.43 0.33 0.26;...

         .998 0.84 0.65 0.5 0.38 0.33; 0.999 0.9 0.77 0.6 0.45 0.4; 1.0 0.97 0.92 .78 .64 .56;...

         1.001 1.002 1.0 1.0 1.0 0.95; 1.002 1.17 1.45 1.75 2.02 2.23];

   

   %

   % Lookup table for variation of induced speed with forward speed and angle of rotor

   %

   % Based on data from Aerodynamics of V/STOL Flight by McCormick

   % Figure 8-6: Variation of Induced Speed with Speed for Constant Ideal Power, page 217

   %

   % x1 = V/w_o,  y1 = theta = rotor angle (alpha in book),  z1 = w/w_o

   %

   x1=[0 1 2 3 4 5;0 1 2 3 4 5; 0 1 2 3 4 5; 0 1 2 3 4 5; 0 1 2 3 4 5];

   y1=[0 0 0 0 0 0; 60 60 60 60 60 60; 70 70 70 70 70 70; 80 80 80 80 80 80; 90 90 90 90 90 90];

   z1=[0.996 0.50 0.25 0.13 0.05 0.03; 0.997 0.63 0.34 0.20 0.12 0.09; 0.998 0.70 0.42 0.25 0.18 0.12;...

         0.999 0.78 0.54 0.37 0.26 0.20; 1.00 0.86 0.68 0.58 0.50 0.46];

   

   

   
   

   %

   % Data from NASA Technical Report TN D-4030: Aerodynamic Data on Large Semispan Tiltwing Wing

   % with 0.5 Diameter Chord, Single-Slotted Flap, and Single Propeller 0.08 Chord Below Wing

   %

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   

   A = 25.16;

   S = 19.6;

   

   alpha_lu = [5 10 15 20 25 30 35 40 45 50 55;5 10 15 20 25 30 35 40 45 50 55;...

         5 10 15 20 25 30 35 40 45 50 55;5 10 15 20 25 30 35 40 45 50 55];

   

   CTs_lu = [0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3;0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6;...

         0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8;0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9 0.9];

   

   % Table 29: Aerodyanmic data for up at tip rotation,

   % full span slat with fences on, delta_flap = 20 deg

   CLs_20_lu = [0.592 1.008 1.445 1.779 2.090 2.258 2.381 2.386 2.301 2.195 1.772;...

         0.541 0.828 1.175 1.470 1.681 1.882 1.994 1.988 1.970 1.938 1.882;...

         0.497 0.697 0.924 1.144 1.327 1.470 1.581 1.639 1.677 1.665 1.631;...

         0.417 0.589 0.756 0.913 1.057 1.185 1.297 1.379 1.425 1.462 1.457];

   

   CDs_20_lu = [-0.252 -.150 -.022 .149 .330 .562 .778 .976 1.118 1.206 1.054;...

         -0.620 -.558 -.435 -.277 -.107 .080 .271 .461 .607 .730 .836;...

         -0.886 -.828 -.741 -.615 -.466 -.307 -.131 .045 .198 .325 .422;...

         -1.013 -.966 -.898 -.797 -.681 -.549 -.408 -.254 -.107 .039 .161];

   

   Cms_20_lu = [-0.095 -.142 -.145 -.144 -.133 -.149 -.140 -.127 -.106 -.067 -.013;...

         -0.038 -.056 -.069 -.074 -.072 -.073 -.064 -.063 -.045 -.017 .021;...

         0.008 .012 .004 -.005 -.023 -.024 -.030 -.031 -.025 -.003 .029;...

         0.038 .041 .039 .029 .032 .018 0 -.004 -.008 -.002 .007];

   % end Table 29

   

   % Table 30: Aerodyanmic data for up at tip rotation,

   % full span slat with fences on, delta_flap = 40 deg

   CLs_40_lu = [0.893 1.348 1.777 2.078 2.285 2.389 2.460 2.382 2.238 2.084 1.49;...

         0.756 1.115 1.452 1.719 1.901 2.040 2.066 2.036 1.967 1.904 1.802;...

         0.631 0.874 1.121 1.339 1.496 1.623 1.676 1.712 1.697 1.653 1.591;...

         0.527 0.706 0.897 1.043 1.198 1.307 1.397 1.470 1.490 1.500 1.465];

   

   CDs_40_lu = [-0.122 0.021 0.188 0.381 0.581 0.815 1.009 1.172 1.249 1.293 1.002;...

         -0.514 -0.399 -0.242 -0.055 0.131 0.326 0.508 0.664 0.769 0.858 0.925;...

         -0.793 -0.704 -0.577 -0.417 -0.245 -0.073 0.101 0.264 0.392 0.481 0.550;...

         -0.947 -0.877 -0.779 -0.653 -0.522 -0.358 -0.193 -0.036 0.110 0.240 0.349];

   

   Cms_40_lu = [-0.181 -0.205 -0.228 -0.218 -0.203 -0.228 -0.193 -0.154 -0.120 -0.059 -0.017;...

         -0.102 -0.187 -0.161 -0.160 -0.157 -0.146 -0.134 -0.109 -0.075 -0.037 0.005;...

         -0.040 -0.058 -0.079 -0.096 -0.106 -0.112 -0.117 -0.105 -0.083 -0.046 -0.010;...

         -0.003 -0.003 -0.017 -0.029 -0.051 -0.065 -0.071 -0.070 -0.066 -0.054 -0.037];

   % end Table 30

   

   % Table 31: Aerodyanmic data for up at tip rotation,

   % full span slat with fences on, delta_flap = 60 deg

   CLs_60_lu = [1.190 1.602 2.009 2.287 2.445 2.481 2.428 2.233 2.070 1.8 1.6;...

         1.054 1.356 1.650 1.869 2.005 2.076 2.055 1.998 1.909 1.810 1.666;...

         0.836 1.100 1.308 1.485 1.603 1.671 1.698 1.697 1.644 1.576 1.491;...

         0.667 0.863 1.032 1.166 1.294 1.382 1.452 1.471 1.466 1.446 1.395];

   

   CDs_60_lu = [0.044 0.205 0.427 0.670 0.881 1.080 1.210 1.249 1.269 1.29 1.32;...

         -0.348 -0.200 -0.012 0.019 0.389 0.570 0.714 0.853 0.928 0.963 0.957;...

         -0.661 -0.517 -0.350 -0.178 -0.013 0.147 0.303 0.452 0.535 0.586 0.636;...

         -0.840 -0.740 -0.592 -0.450 -0.285 -0.118 0.045 0.190 0.320 0.422 0.468];

   

   Cms_60_lu = [-0.241 -0.267 -0.277 -0.280 -0.275 -0.258 -0.195 -0.136 -0.091 -0.05 -0.01;...

         -0.178 -0.203 -0.208 -0.221 -0.217 -0.198 -0.162 -0.136 -0.095 -0.034 0.013;...

         -0.100 -0.133 -0.167 -0.169 -0.170 -0.172 -0.149 -0.142 -0.099 -0.055 -0.029;...

         -0.056 -0.063 -0.093 -0.110 -0.134 -0.140 -0.138 -0.137 -0.124 -0.110 -0.068];

   % end Table 31

   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   

   

   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Takeoff
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   g = 32.2;                  % Gravitational Acceleration
   rho = 0.00226;             % Air Density
   mu_fric = .025;             % Ground Friction Coefficient
   mu_brake = 0.3;            % Braking Friction Coefficient
   
   D_rotor = 2*(0.5*b_wing - 6.5 - 2 - 3);
   A_disk = 2*pi*D_rotor^2/4;
   %A_disk = 673;
   theta_to = wing_angle_to/57.3;
   i_w = 0/57.3;				   % Wing Incidence Angle (degs / 57.3)

   stepsize = 0.5;              % Stepsize for numerical integration (ft/sec)
   %flap_angle_to = 1;
   
   v_to = [0:stepsize:200];

   sg_to(1) = 0;

   
   P_to = 0.95*550*P_to_uninst;
   
   % Calculate static induced velocity & thrust
   w_o = (P_to/(2*rho*A_disk))^(1/3);
   T_o = P_to^(2/3)*(2*rho*A_disk)^(1/3);
   alpha_to = theta_to + i_w; 
   

   for i=2:length(v_to)

      
      % Calculate induced velocity & thrust
      w = w_o*interp2(x1,y1,z1,v_to(i)/w_o,theta_to*180/pi);
      T = T_o*interp2(x,y,z,v_to(i)/w_o,theta_to*180/pi);
      thrust(i) = T;
      
      % Calculation dynamic pressures
      q = 0.5*rho*v_to(i)^2;
      qs = q + (T/A_disk);
      
      % Calculate thrust coefficient
      CTs = T/(qs*A_disk);
      if CTs > 0.8
         CTs = 0.8;
      end
      if CTs < 0.3
         CTs = 0.3;
      end
      
      
      % Find lift and drag coefficients
      A = 25.16;
      S = 19.6;
      CLs_20 = interp2(alpha_lu,CTs_lu,CLs_20_lu,alpha_to*57.3,CTs);
      CDs_20 = interp2(alpha_lu,CTs_lu,CDs_20_lu,alpha_to*57.3,CTs);
      CLs_40 = interp2(alpha_lu,CTs_lu,CLs_40_lu,alpha_to*57.3,CTs);
      CDs_40 = interp2(alpha_lu,CTs_lu,CDs_40_lu,alpha_to*57.3,CTs);
      CLs_60 = interp2(alpha_lu,CTs_lu,CLs_60_lu,alpha_to*57.3,CTs);
      CDs_60 = interp2(alpha_lu,CTs_lu,CDs_60_lu,alpha_to*57.3,CTs);
      CLs_fit = polyfit([20 40 60],[CLs_20 CLs_40 CLs_60],3);
      CDs_fit = polyfit([20 40 60],[CDs_20 CDs_40 CLs_60],3);
      
      if flap_angle_to >= 20
         CLs = polyval(CLs_fit,flap_angle_to);
         CDs = polyval(CDs_fit,flap_angle_to);
      else
         CLs = CLs_fit(1)*flap_angle_to^3 + CLs_fit(2)*flap_angle_to^2 + CLs_fit(3)*flap_angle_to + CLs_fit(4);
         CDs = CDs_fit(1)*flap_angle_to^3 + CDs_fit(2)*flap_angle_to^2 + CDs_fit(3)*flap_angle_to + CDs_fit(4);
      end
      
      CN = CLs/(1-CTs);
      CA = (CDs + CTs*(A/S))/(1-CTs);
      CL_to(i) = CN*cos(alpha_to)-CA*sin(alpha_to);
      CD_to(i) = CN*sin(alpha_to)+CA*cos(alpha_to);
      
      % Calculate drag & lift forces

      drag = q*S_ref*CD_to(i);

      lift = qs*S_ref*CL_to(i);

      

      % Calculate incremental takeoff distances and sum to numerically integrate

      sg_to(i)=sg_to(i-1)+stepsize*((v_to(i))/(g/W_to*(T*cos(theta_to)-drag...
         -mu_fric*(W_to-lift-T*sin(theta_to)))));
      
      % Forces in horizontal & vertical direction
      F_horiz(i) = T*cos(theta_to)-drag-mu_fric*(W_to-lift-T*sin(theta_to));
      F_vert(i) = T*sin(theta_to) + lift - W_to;
      
      % Calculate lift coefficient necessary to takeoff

      CL_reqd(i) = (W_to-T*sin(theta_to))/(qs*S_ref);
      
   end

   
   sg_to(1) = NaN;
   v_to(1) = NaN;
   
   % Calculate distance at which CL_reqd is less than CL_to
   j=length(v_to);
   while CL_reqd(j) > CL_to(j)
      j=j-1;   
   end
   
   liftoffspeed = v_to(j);
   takeoffdist = sg_to(j);
   S_g_to = takeoffdist;
   
   
   %figure(1)
   %plot(v_to,CL_to,v_to,CL_reqd)
   %legend('Operating','Required')
   
   %figure(2)
   %plot(v_to,F_horiz,v_to,F_vert)
   %legend('Horizontal','Vertical')
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % Landing
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   stepsize = 5;
   v = [0:stepsize:100];
   
   
   W_land = 0.85*W_to;           			% Landing Weight
   
   CL_L = 0;                  			% Lift Coefficient During Landing Ground Roll
   CD_L = 0.45;                			% Drag Coefficient During Landing Ground Roll
   
   theta_land = wing_angle_land/57.3;		   		% Rotor/Wing Angle for Landing (degs / 57.3)
   alpha_land = theta_land + i_w; 
   
   P_land = 0.01*P_to;
   w_o = (P_land/(2*rho*A_disk))^(1/3);
   T_o = P_land^(2/3)*(2*rho*A_disk)^(1/3);
   %flap_angle_land = 20;
   
   alpha_land = theta_land + i_w; 
   
   
   for i = 2:length(v)
      
      w = w_o*interp2(x1,y1,z1,v(i)/w_o,theta_land*180/pi);
      T = T_o*interp2(x,y,z,v(i)/w_o,theta_land*180/pi);
      
      % Calculation dynamic pressures
      q = 0.5*rho*v(i)^2;
      qs = q + (T/A_disk);
      
      % Calculate thrust coefficient
      CTs = T/(qs*A_disk);
      if CTs > 0.8
         CTs = 0.8;
      end
      if CTs < 0.3
         CTs = 0.3;
      end
      
      
      % Find lift and drag coefficients
      A = 25.16;
      S = 19.6;
      CLs_20 = interp2(alpha_lu,CTs_lu,CLs_20_lu,theta_land*57.3,CTs);
      CDs_20 = interp2(alpha_lu,CTs_lu,CDs_20_lu,theta_land*57.3,CTs);
      CLs_40 = interp2(alpha_lu,CTs_lu,CLs_40_lu,theta_land*57.3,CTs);
      CDs_40 = interp2(alpha_lu,CTs_lu,CDs_40_lu,theta_land*57.3,CTs);
      CLs_60 = interp2(alpha_lu,CTs_lu,CLs_60_lu,theta_land*57.3,CTs);
      CDs_60 = interp2(alpha_lu,CTs_lu,CDs_60_lu,theta_land*57.3,CTs);
      CLs_fit = polyfit([20 40 60],[CLs_20 CLs_40 CLs_60],3);
      CDs_fit = polyfit([20 40 60],[CDs_20 CDs_40 CLs_60],3);
      if flap_angle_land > 20
         CLs = polyval(CLs_fit,flap_angle_land);
         CDs = polyval(CDs_fit,flap_angle_land);
      else
         CLs = CLs_fit(1)*flap_angle_land^3 + CLs_fit(2)*flap_angle_land^2 + CLs_fit(3)*flap_angle_land + CLs_fit(4);
         CDs = CDs_fit(1)*flap_angle_land^3 + CDs_fit(2)*flap_angle_land^2 + CDs_fit(3)*flap_angle_land + CDs_fit(4);
      end
      CN = CLs/(1-CTs);
      CA = (CDs + CTs*(A/S))/(1-CTs);
      CL(i) = CN*cos(theta_land)-CA*sin(theta_land);
      CD(i) = CN*sin(theta_land)+CA*cos(theta_land);
      
      % Calculate drag & lift forces

      drag = q*S_ref*CD(i);

      lift = qs*S_ref*CL(i);
      
      vert(i) = lift + T*sin(theta_land) - W_land;
      gamma(i) = 57.3*asin((T*cos(theta_land) - drag)/W_land);
      
   end
   
   v_data = polyfit(vert,v,3);
   gamma_data = polyfit(v,gamma,3);
   touchdown = polyval(v_data,0);

   flightpath = polyval(gamma_data,touchdown);
   
   V_td=touchdown;                         % Calculate touchdown speed
   v=[V_td:-stepsize:0];                   % Numerically integrate from touchdown speed to zero
   sg_l(1)=V_td;                           % Assume 1 second for rotation during landing
   
   for i=2:length(v)
      T=0;                                 % Assume there is no thrust during landing ground roll
      drag=.5*rho*v(i)^2*S_ref*CD_L;       % Calculate drag during landing
      % lift=.5*rho*v(i)^2*S_ref*CL_L;     % Calculate lift during landing
      lift = 0;                            % Assume there is no lift during landing ground roll
      % Calculate incremental landing distances and sum to numerically integrate
      sg_l(i)=sg_l(i-1)-stepsize*(v(i)/(g*((T/W_land)+(lift*mu_brake-drag)/W_land+0.08*(mu_brake-mu_fric)-mu_brake)));
   end
   
   touchdown;
   flightpath;
   
   if length(v) > 0
      landingdist = sg_l(length(v));
   else
      landingdist = 0;
   end
   
   S_g_land = landingdist;
   
   P_to = P_to/550;
   T_to = P_to;
   
end