%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Messerschmitt Me-262 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc
format short g

% Aircraft data
b = 12.50              % wingspan (m)
S = 21.70              % wing area (m^2)
hw = 1.5               % wing height from ground (m)
Wg = 6400*9.8          % gross weight (N)
W = 5000*9.8           % average weight (N)
We = 3800*9.8          % empty weight (N)
cL_max = 2.0           % max lift coefficient
cDo = 0.020            % zero-lift drag coefficient
e = 0.7                % Oswald efficiency factor
AR = b^2/S             % aspect ratio
Ta_max = 2*900*9.8     % max trhust available (N)

%% Mission briefing:
%% 1) Takeoff
%% 2) Climb at 6,000 meters
%% 3) Cruise
%% 4) Gliding
%% 5) Landing

%% 1) Calculation of the takeoff distance on runway
rho = 1.225                         % air density at sea level
g = 9.8                             % gravity acceleration
mr = 0.02                           % runway friction coefficient
Vmin = sqrt(2*Wg/(rho*S*cL_max))    % stall speed
Vto = 1.2 * Vmin                    % takeoff speed

phi = (16*hw/b)^2/(1+(16*hw/b)^2)   % coefficient due to ground effect
D = 0.5*rho*(0.7*Vto)^2*S*(cDo+phi*cL_max^2/(pi*AR*e))  % average drag
L = 0.5*rho*(0.7*Vto)^2*S*cL_max                        % average lift

dto = 1.44*Wg^2/(rho*S*cL_max*g*(Ta_max-D-mr*(Wg-L)))   % takeoff distance

%% 2) Calculation of the max rate of climb
h = 6000                                 % cruise altitude
rho_h = 0.66011                          % air density at 6000 meters

for i = 1:300
    V(i) = i;
    q(i) = rho*V(i)^2/2;                 % dynamic pressure
    cL(i) = Wg/(q(i)*S) ;                % lift coefficient
    cD(i) = cDo + cL(i)^2/(pi*e*AR);     % drag coefficient
    Pr_zl(i) = q(i)*S*V(i)*cDo;          % power due to zero-lift drag
    Pr_id(i) = q(i)*S*V(i)*cL(i)^2/(pi*AR*e); % power due to induced drag
    Pr(i) = Pr_zl(i) + Pr_id(i);         % power required at sea level
    Pa(i) = Ta_max * V(i);               % MAX power available at sea level
    EoP(i) = (Pa(i) - Pr(i));            % excess of power
    RC(i) = (EoP(i))/Wg;                 % R/C = excess of power / weight
    RC_mpm(i) = RC(i) * 60;              % meters per minute
    RC_mpm_h(i) = RC_mpm(i) * rho_h/rho; % rate of climb at 6000 meters
end

RCmax = max(EoP)/Wg                      % max rate of climb at sea level
RCmax_mpm = RCmax * 60
RCmax_h = RCmax * rho_h/rho              % max rate of climb at 6000 meters
RCmax_h_mpm = RCmax_h * 60

plot (V,Pa, 'r.')
hold on
plot (V,Pr_zl, 'g.')
hold on
plot (V,Pr_id, 'b.')
hold on
plot (V,Pr, 'k.')
grid on
title ('Me-262 - Power required and power available at sea level')
xlabel ('Velocity (m/s)')
ylabel ('Power (W)')
axis ([0 300 0 5e6])
legend ('MAX power available', 'power due to zero-lift drag',...
    'power due to induced drag', 'power required')

figure
plot (V,RC_mpm, 'g.')
hold on
plot (V,RC_mpm_h, 'b.')
grid on
title ('Me-262 - Rate of Climb')
xlabel ('Velocity (m/s)')
ylabel ('Rate of Climb (m/min)')
axis ([0 300 0 2000])
legend ('sea level', '6000 meters')

av_RC = (RCmax + RCmax_h)/2       % average rate to climb
av_RC_mpm = av_RC * 60
time = h/av_RC                    % approx. time to climb

%% 3) Cruise flight
Ta = 0.8 * Ta_max * rho_h/rho     % cruise thrust available at 6000 meters

for i = 1:300
    V(i) = (i+34)*sqrt(rho/rho_h);
    q(i) = rho_h*V(i)^2/2;
    Tr_zl(i) = q(i)*S*cDo;             % thrust due to zero-lift drag
    Tr_id(i) = W^2/(q(i)*S*pi*e*AR);   % thrust due to induced drag
    Tr(i) = Tr_zl(i) + Tr_id(i);       % thrust required
    cL(i) = W/(q(i)*S);
    cD(i) = cDo + cL(i)^2/(pi*e*AR);
    E(i) = cL(i)/cD(i);                % aerodynamic efficiency
end

figure
plot (V,Tr_zl, '.g')
hold on
plot (V,Tr_id, '.r')
hold on
plot (V,Tr, '.k')
hold on
plot (V,Ta, '.b')
grid on
axis ([0 300 0 1.5e4])
title ('Me-262 - Thrust required curves at 6000 meters')
xlabel ('Velocity (m/s)')
ylabel ('Thrust (N)')
legend ('thrust due to zero-lift drag','thrust due to induced drag',...
    'thrust required', 'thrust available')

figure
plot (V,E, '.')
grid on
axis ([0 300 0 20])
title ('Me-262 - Aerodynamic efficiency')
xlabel (' Velocity (m/s)')
ylabel ('E = cL/cD')

Tmin = min(Tr)                    % min thrust
Emax = max(E)                     % max aerodynamic efficiency
VTmin = sqrt(Tmin/(rho_h*S*cDo))  % min thrust (max efficiency) airspeed
VTmin_KPH = VTmin * 3.6

%% 4) Calculation of the max range in gliding flight
theta = atand(1/Emax)             % glide angle at max efficiency
range = h * Emax                  % max range in gliding at 6000 meters

%% 5) Calculation of the landing distance on the grass
Vmin = sqrt(2*We/(rho*S*cL_max))  % stall speed
Vl = 1.3 * Vmin                   % landing speed
mg = 0.60                         % grass + brake friction coefficient

D = 0.5*rho*(0.7*Vl)^2*S*(cDo+phi*cL_max^2/(pi*e*AR))      % average drag
L = 0.5*rho*(0.7*Vl)^2*S*cL_max                            % average lift

dl = 1.69*We^2/(g*rho*S*cL_max*(D+mg*(W-L)))           % landing distance