% Performance Sizing for SSTOL Fixed Wing Aircraft
%
% Patrick LeGresley
% February 2, 1999
%
% Last modified March 3, 1999

close all
clear all


% Define/Specify Constants
W_to = 48650;              % Takeoff Weight
W_L = 0.85*W_to;           % Landing Weight
g = 32.2;                  % Gravitational Acceleration
CL_max_L = 7;              % Maximum Lift Coefficient (Landing Configuration)
CL_max_to=0.7*CL_max_L;    % Maximum Lift Coefficient for TO configuration
rho = 0.00226;             % Air Density
CL_to = .6;                % Lift Coefficient During TO Ground Roll
CD_to = CL_to/8;           % Drag Coefficient During TO Ground Roll
LD_app = 6;                % Lift to Drag Ratio in Approach Configuration
CL_L = 0;                  % Lift Coefficient During Landing Ground Roll
CD_L = 0.25;               % Drag Coefficient During Landing Ground Roll
mu_fric = .025;            % Ground Friction Coefficient
mu_brake = 0.3;            % Braking Friction Coefficient
stepsize = 0.1;            % Set stepsize for numerical integration

a = sqrt(1.4*1716*(59+459.7+27));

% Make calculations for wing loadings from 5 lb/ft^2 to 60 lb/ft^2

n = 0;

for WS=5:5:60
   
n = n+1;
wingload(n) = WS;   
S_ref = W_to/WS;                        % Calculate wing area from takeoff wing loading

Vs_to = sqrt(2*W_to/(rho*S_ref*CL_max_to)); % Stall speed in TO configuration
V_lo=1.15*Vs_to;                        % Calculate liftoff speed

takeoffdist = 1000;                     % Initialize takeoff distance to start loop
TW = 0.05;                              % Initialize T/W to start loop

% Iterate on T/W to find minimum value to satisfy takeoff distance requirement

while takeoffdist > 300
   
   if takeoffdist > 300
      TW = TW + 0.01;
   else
      TW = TW - 0.01;
   end
   
T_to = W_to*TW;   
   
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*(-(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

takeoffdist = sg_to(length(v));

end
thrustweight(n) = TW;


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

landing(n) = sg_l(length(v));           % Output landing ground roll
end

figure(1)
plot(wingload,thrustweight);
xlabel('Wing Loading (lb/sq ft)');
ylabel('Thrust to Weight Ratio');
t1 = sprintf('%s%1.1f','Performance Matching Plot for CL_{max_{to}} = ',CL_max_to);
t2 = sprintf('Takeoff Distance = 300 ft');
title(t1);
text(5,0.8,t2);

figure(2)
plot(wingload,landing);
xlabel('Wing Loading (lb/sq ft)');
ylabel('Landing Distance (ft)');
t3 = sprintf('%s%1.1f','Performance Matching Plot for CL_{max_{L}} = ',CL_max_L);
title(t3);
