% Script to explore Glide analysis
% Data is from Baron 58 POH
% CD0 e and k are computed using bootstrap equations as per experiment writeup.

clear all
close all
A=7.17 % Aspect Ratio
S=199.2 % wing area, ft^2
Wbg=5500
Vibg=115 % best glide speed knots
flap=0
Vcbg=AS1Baron58(Vibg,flap)
Hpcalibrated=1000
[temp,press,rho,Hgeopvector]=atmosphere4(Hpcalibrated,0)
[Vebg,qcOverp0,qcOverp]=AS2(Vcbg,press)
[Vtbg,M,aknots,rho,sigma]=AS3(Vebg,temp,Hpcalibrated)
VtbgFPS=Vtbg*1.687808 % Convert knots to ft/sec       ***


Hdrop=950:10:1050;
n=length(Hdrop);
CD0vect=zeros(size(Hdrop));
evect=zeros(size(Hdrop));
kvect=zeros(size(Hdrop));

for i=1:n   
    tanGbg=Hdrop(i)/(2*6076.1); % 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*A*(tan(Gbg))^2) ;                        % ***
    evect(i)=e;
    k=1/(pi*A*e);
    kvect(i)=k;
end

% Compute nominal values for CD0 and e and k using drop altitude of 1000 ft
HdropNominal=1000;
tanGbg=HdropNominal/(2*6076.1); % 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*A*(tan(Gbg))^2)                        % ***
kNominal=1/(pi*A*eNominal)


strT=['Glide ratio is approx. 2 nm of gliding distance for each 1000 ft of altitude drop at Vi=115 knots'];
figure(1)
plot(Hdrop,CD0vect,'-',HdropNominal,CD0Nominal,'rx')
xlabel('Altitude dropped in 2 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 2 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 2 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)
