% Script to explore Glide analysis
% Data is from Cessna 182P Owner's manual
% CD0 e and k are computed using bootstrap equations as per experiment writeup.
% The drag polar will be based on dropping 16000 ft in 30 sm at Vi=80 mph


clear all; close all
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<')

echo on; format compact
% Bootstrap Data Plate
S=174;         % wing area (ft^2)
B=36;          % Wing span (ft)
BHPmax=230;    % Rated engine horsepower (hp), See Owner's Manual p. 4-2
RPM=2600;      % Rated engine speed, RPM
Din=82;        % Propeller diameter, in
%CD0=Unknown;  % Paracite drag coefficient
%e=Unknown;    % Airplane efficiency factor
%m=Unknown;    % Slope of propeller polar (from Step 2)
%b=Unknown;    % Intercept of propeller polar (from Step 2)

% Auxilary data
W0=2950;       % Standard weight, lbf
Wbg=2950;      % Weight for which Best glide os known, lbf
AR=B*B/S       % Aspect ratio
%k=unknown=1/(pi*AR*e)  % Parameter of drag polar
nblades=2;     % Number of blades of propeller
C=.12;         % Altitude drop-off parameter for Gagg-Farrar Rule
echo off; 

mph2k=0.868976;
k2mph=1/0.868976;


disp('The drag polar will be based on dropping 16000 ft in 30 sm at Vi=80 mph.')
Vibgmph=80 % indicated best glide speed, mph
Vibg= Vibgmph*0.868976; % indicated best glide speed knots
flap=0;
% function VcalibratedKnots=AS1Cessna182(VindicatedKnots,flap)
Vcbg=AS1Cessna182(Vibg,flap);   % Calibrated best glide speed, knots
Hpcalibrated=5000;
[temp,press,rho,Hgeopvector]=atmosphere4(Hpcalibrated,0);
[Vebg,qcOverp0,qcOverp]=AS2(Vcbg,press);
[Vtbg,M,aknots,rho2,sigma]=AS3(Vebg,temp,Hpcalibrated);
VtbgFPS=Vtbg*1.687808; % Convert knots to ft/sec       ***


Hdrop=15000:100:17000;
n=length(Hdrop);
CD0vect=zeros(size(Hdrop));
evect=zeros(size(Hdrop));
kvect=zeros(size(Hdrop));

for i=1:n   
    tanGbg=Hdrop(i)/(30*5280); % tangent of best glide angle
    Gbg=atan(tanGbg);
    CD0=Wbg*sin(Gbg)/(rho*S*VtbgFPS*VtbgFPS); % original equation
    CD0=Wbg*tan(Gbg)/(rho*S*VtbgFPS*VtbgFPS); % modified   ***
    CD0vect(i)=CD0;
    e=4*CD0/(pi*AR*(tan(Gbg))^2) ;                        % ***
    evect(i)=e;
    k=1/(pi*AR*e);
    kvect(i)=k;
end

% Compute nominal values for CD0 and e and k using drop altitude of 1000 ft
HdropNominal=16000;
tanGbg=HdropNominal/(30*5280); % tangent of best glide angle
Gbg=atan(tanGbg);
CD0=Wbg*sin(Gbg)/(rho*S*VtbgFPS*VtbgFPS); % original equation
CD0=Wbg*tan(Gbg)/(rho*S*VtbgFPS*VtbgFPS); % modified   ***
CD0Nominal=CD0
eNominal=4*CD0Nominal/(pi*AR*(tan(Gbg))^2)                        % ***
kNominal=1/(pi*AR*eNominal)


strT=['Glide ratio is approx. 30 sm of gliding distance for each 16000 ft of altitude drop at Vi=80mph'];
figure(1)
plot(Hdrop,CD0vect,'-',HdropNominal,CD0Nominal,'rx')
xlabel('Altitude dropped in 30 sm (ft)')
ylabel('CD0')
srt1=['Nominal CD0 is ',num2str(CD0Nominal),' at HdropNominal of ',num2str(HdropNominal),' feet.'];
text2(.45,.52,srt1)
title('Glide analysis using data from POH')
text2(.02,.95,strT)

figure(2)
plot(Hdrop,evect,'-',HdropNominal,eNominal,'rx')
xlabel('Altitude dropped in 30 sm (ft)')
ylabel('e')
srt2=['Nominal e is ',num2str(eNominal),' at HdropNominal of ',num2str(HdropNominal),' feet.'];
text2(.45,.52,srt2)
title('Glide analysis using data from POH')
text2(.02,.95,strT)

figure(3)
plot(Hdrop,kvect,'-',HdropNominal,kNominal,'rx')
xlabel('Altitude dropped in 30 sm (ft)')
ylabel('k')
srt3=['Nominal k is ',num2str(kNominal),' at HdropNominal of ',num2str(HdropNominal),' feet.'];
text2(.45,.52,srt3)
title('Glide analysis using data from POH')
text2(.02,.95,strT)
