%% Calculating Best Glide Quantities
% This example shows how to perform  glide calculations for a Cessna 172 following 
% Example 9.1 in Reference 1 using the Aerospace Toolbox(TM) software. 
%
% Best glide calculations provide values (velocity and glide angle) that 
% minimize drag and maximize lift-drag ratio (also called the glide ratio).

% Copyright 2006-2019 The MathWorks, Inc.

%% Aircraft Specifications
% The aircraft parameters are declared as follows. 
W = 2400; % weight, lbf
S = 174;  % wing reference area, ft^2;
A = 7.38; % wing aspect ratio
C_D0 = 0.037; % flaps up parasite drag coefficient
e = 0.72; % airplane efficiency factor
%% Conditions
% Set the current aircraft conditions. The bank angle (phi) is zero for this case.
h = 4000; % altitude, ft
phi = 0; % bank angle, deg
%%
% Convert altitude to meters using |convlength|. The atmospheric calculations in the next step require values in metric units.
h_m = convlength(h,'ft','m');
%%
% Calculate atmospheric parameters based on altitude using |atmoscoesa|:
[T, a, P, rho] = atmoscoesa(h_m, 'Warning');
%%
% Convert density from metric to English units using |convdensity|:
rho = convdensity(rho,'kg/m^3','slug/ft^3');
%% Best Glide Data
% Best glide velocity is calculated using the following equation. TAS (true airspeed in feet per second) is the velocity of the aircraft relative to the surrounding air mass.
%
% $$TAS_{bg}=\sqrt{\frac{2W}{\rho S}}\times \lbrack\frac{1}{4{C_{D_0}}^2+C_{D_0}\pi e A\cos^2{\phi}}\rbrack^\frac{1}{4}$$
%
TAS_bg = sqrt((2*W) / (rho*S))...
         *(1./(4*C_D0.^2 + C_D0.*pi*e*A*cos(phi)^2)).^(1/4); % TAS, fps
%%
% Convert velocity from fps to kts using |convvel|. KTAS is true airspeed in knots.
KTAS_bg = convvel(TAS_bg,'ft/s','kts')';
%%
% Convert KTAS to KCAS using |correctairspeed|. KCAS (calibrated airspeed in knots) is the velocity corrected for instrument error and position error. This position error comes from inaccuracies in static pressure measurements at different points in the flight envelope.
KCAS_bg = correctairspeed(KTAS_bg,a,P,'TAS','CAS')';
%%
% Best glide angle is calculated using:
%
% $$\sin{\gamma_{bg}}=-\sqrt{\frac{4C_{D_0}}{\pi e A\cos^2{\phi}+4C_{D_0}}}$$
%
% This is the angle between the flight path and the ground that provides the highest L/D ratio.
gamma_bg_rad = asin( -sqrt((4.*C_D0')./(pi*e*A*cos(phi)^2 + 4.*C_D0')) );
%%
% Convert glide angle from radians to degrees using |convang|:
gamma_bg = convang(gamma_bg_rad,'rad','deg');
%%
% Best glide drag is calculated using:
%
% $$D_{min}=D_{bg}=\frac{1}{2}\rho (TAS^2_{bg})S(2C_{D_0})=-W\sin{\gamma_{bg}}$$
%
D_bg = -W*sin(gamma_bg_rad);
%%
% Best glide lift is calculated using:
%
% $$L_{bg}=L_{max}=W\cos{\gamma_{bg}}=\sqrt{W^2-D^2_{bg}}$$
%
L_bg =  W*cos(gamma_bg_rad);
%%
% Calculate dynamic pressure using |dpressure|:
qbar = dpressure([TAS_bg' zeros(size(TAS_bg,2),2)], rho);
%%
% Calculate drag and lift coefficients using:
%
% $$C_{D_{bg}}=\frac{D_{bg}}{\bar{q}S}$$
%
% $$C_{L_{bg}}=\frac{L_{bg}}{\bar{q}S}$$
%
C_D_bg = D_bg./(qbar*S);
C_L_bg = L_bg./(qbar*S);
%% Summary of Best Glide Values
% Here are the best glide values:
% 
% $$KCAS_{bg}=71.9 \hspace{1pt} KCAS$$
%
% $$\gamma_{bg}=-5.38\deg$$
%
% $$C_{D_{bg}}=0.074$$
%
% $$C_{L_{bg}}=0.7859$$
%
% $$D_{bg}=224.9 \hspace{1pt} lbf$$
%
% $$L_{bg}=2389.4 \hspace{1pt} lbf$$
%
%% Verification
% These plots show drag and lift-drag ratio plots for the aircraft as a function of KCAS. The plots are used to verify the best glide calculations.
%%
% Set range of airspeeds and convert to KCAS using |convvel| and |correctairspeed|:
TAS = (70:200)'; % true airspeed, fps
KTAS = convvel(TAS,'ft/s','kts')'; % true airspeed, kts
KCAS = correctairspeed(KTAS,a,P,'TAS','CAS')'; % corrected airspeed, kts
%%
% Calculate dynamic pressure for new airspeeds using |dpressure|:
qbar = dpressure([TAS zeros(size(TAS,1),2)], rho);
%%
% Calculate parasite drag using:
%
% $$D_p=\frac{1}{2}\rho S C_{D_0}(TAS^2)$$
%
Dp = qbar*S.*C_D0;
%% 
% Calculate induced drag using:
%
% $$D_i=\frac{2W^2}{\rho S\pi e A}\frac{1}{(TAS^2)}$$
%
Di = (2*W^2)/(rho*S*pi*e*A).*(TAS.^-2);
%%
% Calculate total drag using:
%
% $$D=D_p+D_i$$
%
D = Dp + Di;
%%
% Approximate lift as weight (assuming small glide angle and small angle of attack). At this speed, assuming
%
% $$C_L=2\pi\alpha$$
%
% and using
%
% $$C_{L_{bg}}$$
%
% from above, the angle of attack is about 7 degrees. Adding the flight path angle (i.e. best glide angle) from above shows the fuselage pitch (attitude angle theta) to be about 2 degrees.
L = W;
%%
% *Plot L/D versus KCAS*
%
% As expected, the maximum L/D occurs at approximately the best glide velocity calculated above.
h1 = figure;
plot(KCAS,L./D);
title('L/D vs. KCAS');
xlabel('KCAS'); ylabel('L/D');
hold on
plot(KCAS_bg,L_bg/D_bg,'Marker','o','MarkerFaceColor','black',...
    'MarkerEdgeColor','black','Color','white');
hold off
legend('L/D','L_{bg}/D_{bg}','Location','Best');
annotation('textarrow',[0.49 0.49],[0.23 0.12],'String','KCAS_{bg}');
%%
% *Plot parasite, induced, and total drag curves*
%
% Notice the minimum total drag (i.e. D_bg) occurs at approximately the same best glide velocity calculated above.
h2 = figure;
plot(KCAS,Dp,KCAS,Di,KCAS,D); 
title('Parasite, induced, and total drag curves');
xlabel('KCAS'); ylabel('Drag, lbf'); 
hold on
plot(KCAS_bg,D_bg,'Marker','o','MarkerFaceColor','black',...
    'MarkerEdgeColor','black','Color','white');
hold off
legend('Parasite, D_p','Induced, D_i','Total, D','D_{bg}','Location','Best');
annotation('textarrow',[0.49 0.49],[0.23 0.12],'String','KCAS_{bg}');
%%
%
close(h1,h2);

%% Reference
%
%  [1] Lowry, J. T., "Performance of Light Aircraft", AIAA(R) Education Series,
%      Washington, DC, 1999.

%%
displayEndOfDemoMessage(mfilename)
