% Script to generate simulated glide data at different airspeeds
% In Step_1_Glide_Baron58.m we used the glide data from the Baron 58 POH to
% compute a drag polar. That information was that the aircraft when gliding would
% drop 1000 ft in 2 nm at Vi=115 kts. Using the Bootstrap approach we
% computed a drag polar. In this script we use that drag polar to 
% generate rate of descent data at other airspeeds. This simulates data
% we would like to have from sawtooth climb/glide flight testing.


clear all
A=7.17 % Aspect Ratio
S=199.2 % wing area, ft^2
W=5500

CD0=0.02538  %drag polar based on dropping 1000 ft in 2 nm at Vi=115 kts
e=0.66556    %drag polar based on dropping 1000 ft in 2 nm at Vi=115 kts

k=1/(pi*A*e)
Viknots=100:5:130;
Hipress=2500*ones(size(Viknots));
flapVector=0*ones(size(Viknots));
str1=['Drag Polar: CD0=',num2str(CD0),', e=',num2str(e),', k=',num2str(k),','];

VcVector=AS1Baron58(Viknots,flapVector);
HcVector=ALTBaron58(Hipress,Viknots,flapVector); % a vector same length as Hp2
[tVector,pVector,rho,Hgeopvector]=atmosphere4(HcVector,0); % altitude is geopotential, p is a vector
[VeVector,qcOverp0,qcOverp]=AS2(VcVector,pVector);
[VtrueKnots,M,aknots,rho,sigma]=AS3(VeVector,tVector,HcVector);
VtFtperSec=VtrueKnots*1.687808; % Convert knots to ft/sec

cosgamma=ones(size(Viknots));
Lift=W./cosgamma;
CL=2*Lift./(rho.*VtFtperSec.*VtFtperSec*S);
CD=CD0+(1/(pi*A*e))*CL.*CL;
Drag=0.5*rho.*VtFtperSec.*VtFtperSec*S.*CD;
dVdHpersec=dVdH(Viknots,flapVector,Hipress);
% See equation 12.11 page 122 of Kimberlin.
multiplier=(1+VtFtperSec.*dVdHpersec/32.17);
%Drag=-Wt.*singamma.*multiplier;    
hdot=(-Drag.*VtFtperSec)./(W*multiplier);

% iterate second time on gamma
singamma=hdot./VtFtperSec;
gamma=asin(singamma);
cosgamma=cos(gamma);
Lift=W./cosgamma;
CL=2*Lift./(rho.*VtFtperSec.*VtFtperSec*S);
CD=CD0+(1/(pi*A*e))*CL.*CL;
Drag=0.5*rho.*VtFtperSec.*VtFtperSec*S.*CD;
dVdHpersec=dVdH(Viknots,flapVector,Hipress);
% See equation 12.11 page 122 of Kimberlin.
multiplier=(1+VtFtperSec.*dVdHpersec/32.17);
%Drag=-Wt.*singamma.*multiplier;    
hdot=(-Drag.*VtFtperSec)./(W*multiplier);

% iterate a third time on gamma
singamma=hdot./VtFtperSec;
gamma=asin(singamma);
cosgamma=cos(gamma);
Lift=W./cosgamma;
CL=2*Lift./(rho.*VtFtperSec.*VtFtperSec*S);
CD=CD0+(1/(pi*A*e))*CL.*CL;
Drag=0.5*rho.*VtFtperSec.*VtFtperSec*S.*CD;
dVdHpersec=dVdH(Viknots,flapVector,Hipress);
% See equation 12.11 page 122 of Kimberlin.
multiplier=(1+VtFtperSec.*dVdHpersec/32.17);
%Drag=-Wt.*singamma.*multiplier;    
hdot=(-Drag.*VtFtperSec)./(W*multiplier);


GammaGlide=57.3*asin(hdot./VtFtperSec);


figure(2); clf
plot(Viknots,hdot*60)
title('Unpowered glide analysis')
xlabel('Vi knots')
ylabel('hdot ft/min')

tanGbg=1000/(2*6076.1) % tangent of best glide angle
Gbg=-atan(tanGbg)
Vibg=115 %knots
Vc=AS1Baron58(Vibg,0);
Hc=ALTBaron58(Hipress(1),Vibg,0); % a vector same length as Hp2
[t,p,rho,Hgeopvector]=atmosphere4(Hc,0); % altitude is geopotential, p is a vector
[Ve,qcOverp0,qcOverp]=AS2(Vc,p);
[VtrueKnots,M,aknots,rho,sigma]=AS3(Ve,t,Hc);
Vt=VtrueKnots*1.687808; % Convert knots to ft/sec
hdotPOH=Vt*sin(Gbg);
hold on; plot(Vibg,hdotPOH*60,'rx'); 
text2(.02,.15,str1)
text2(.02,.1,'The red x is from dropping 1000 ft in 2 nm at Vi=115kts.')
text2(.02,.05,'The red x does not consider the dV/dH term from gliding at constant Vi (not constant Vt).')
hold off



% ****
figure(3); clf
plot(VtFtperSec,GammaGlide)
title('Unpowered glide analysis')
xlabel('Vtrue ft/sec')
ylabel('gamma deg')

tanGbg=1000/(2*6076.1) % tangent of best glide angle
Gbg=-atan(tanGbg)
Vibg=115 %knots
Vc=AS1Baron58(Vibg,0);
Hc=ALTBaron58(Hipress(1),Vibg,0); % a vector same length as Hp2
[t,p,rho,Hgeopvector]=atmosphere4(Hc,0); % altitude is geopotential, p is a vector
[Ve,qcOverp0,qcOverp]=AS2(Vc,p);
[VtrueKnots,M,aknots,rho,sigma]=AS3(Ve,t,Hc);
Vt=VtrueKnots*1.687808; % Convert knots to ft/sec
hdotPOH=Vt*sin(Gbg);
hold on; plot(Vt,Gbg*57.3,'rx'); text2(.05,.95,str1); 
text2(.02,.1,'The red x is from dropping 1000 ft in 2 nm at Vi=115kts.')
text2(.02,.05,'The red x does not consider the dV/dH term from gliding at constant Vi (not constant Vt).')
hold off

% ****
figure(1); clf
plot(Viknots,GammaGlide)
title('Unpowered glide analysis')
xlabel('Vi knots')
ylabel('gamma deg')

tanGbg=1000/(2*6076.1) % tangent of best glide angle
Gbg=-atan(tanGbg)
Vibg=115 %knots
Vc=AS1Baron58(Vibg,0);
Hc=ALTBaron58(Hipress(1),Vibg,0); % a vector same length as Hp2
[t,p,rho,Hgeopvector]=atmosphere4(Hc,0); % altitude is geopotential, p is a vector
[Ve,qcOverp0,qcOverp]=AS2(Vc,p);
[VtrueKnots,M,aknots,rho,sigma]=AS3(Ve,t,Hc);
Vt=VtrueKnots*1.687808; % Convert knots to ft/sec
hdotPOH=Vt*sin(Gbg);
hold on; plot(Vibg,Gbg*57.3,'rx'); text2(.05,.95,str1); 
text2(.02,.1,'The red x is from dropping 1000 ft in 2 nm at Vi=115kts.')
text2(.02,.05,'The red x does not consider the dV/dH term from gliding at constant Vi (not constant Vt).')
hold off

% Write data to command window
n=length(Hipress)
TdegC=(5/9)*(tVector-32-459.67);
hdot=hdot*60; % Convert from ft/sec to ft/min
str2=['% Hp(ft)   OAT (degC)   Hdot(ft/min) Wt(lbf) Vi(knots)'];
disp(str2)
for i=1:n
    strx=[num2str(Hipress(i)),'       ',num2str(TdegC(i)),'       ',num2str(hdot(i)),'       ',num2str(W),'       ',num2str(Viknots(i))];
    disp(strx)
end

disp(' '); disp('Copy the above data into a text file called Baron_Glide_Data.txt.')
disp('This file is the input file for Step_1Alternate_Glide.m.')