% Script to explore Glide analysis
% Data is from Piper Warrior III Information Manual manual
% CD0 e and k are computed using bootstrap equations as per experiment writeup.
% The drag polar will be based on dropping 12000 ft in 22.75 nm at Vi=73
% KIAS


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

echo on; format compact
% Bootstrap Data Plate: Warrior 3
S=170;         % wing area (ft^2)
B=35;          % Wing span (ft)
BHPmax=160;    % Rated engine horsepower (hp), See Owner's Manual p. 4-2
RPM=2700;      % Rated engine speed, RPM
Din=74;        % 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=2440;       % Standard weight, lbf: Max TO wt normal category
Wbg=2440;      % Weight for which Best glide os known, lbf, page 5-27
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; 


disp('The drag polar will be based on dropping 12000 ft in 22.75 sm at Vi=73 KIAS.')
DropDist=22.75; %nm drop 12000 feet in this distance
Vibg= 73 % indicated best glide speed knots
flap=0;
% function VcalibratedKnots=AS1Warrior3(VindicatedKnots,flap)
Vcbg=AS1Warrior3(Vibg,flap)   % Calibrated best glide speed, knots
Hpcalibrated=0;
[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=11000:100:13000;
n=length(Hdrop);
CD0vect=zeros(size(Hdrop));
evect=zeros(size(Hdrop));
kvect=zeros(size(Hdrop));

for i=1:n   
    tanGbg=Hdrop(i)/(DropDist*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=12000;
tanGbg=HdropNominal/(DropDist*5280); % tangent of best glide angle
Gbg=atan(tanGbg)
GbgDEG=Gbg*180/pi
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. 22.75 nm of gliding distance for each 12000 ft of altitude drop at Vi=73 KIAS'];
figure(1)
plot(Hdrop,CD0vect,'-',HdropNominal,CD0Nominal,'rx')
xlabel('Altitude dropped in 22.75 nm (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 22.75 nm (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 22.75 nm (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)
