% Script to generate simulated glide data at different airspeeds
% In Step_1_Glide_Cessna182.m we used the glide data from the Cessna 18a Owner's Manual to
% compute a drag polar. That information was that the aircraft when gliding would
% drop 16000 ft in 30 sm at Vi=80 mph. 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; 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=0.048649;  % Paracite drag coefficient based on dropping 16000 ft in 30 sm at Vi=80 mph
e=0.81507;     % Airplane efficiency factor based on dropping 16000 ft in 30 sm at Vi=80 mph
%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
W=2950;        % Test day weight lbf
AR=B*B/S       % Aspect ratio
k=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;


Vimph=70:5:95;   % mph
Viknots=Vimph*0.86897;
Hipress=5000*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
% function VcalibratedKnots=AS1Cessna182(VindicatedKnots,flap)
VcVector=AS1Cessna182(Viknots,flapVector);
HcVector=Hipress;
[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*AR*e))*CL.*CL;
Drag=0.5*rho.*VtFtperSec.*VtFtperSec*S.*CD;
dVdHpersec=dVdH_Cessna(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*AR*e))*CL.*CL;
Drag=0.5*rho.*VtFtperSec.*VtFtperSec*S.*CD;
dVdHpersec=dVdH_Cessna(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*AR*e))*CL.*CL;
Drag=0.5*rho.*VtFtperSec.*VtFtperSec*S.*CD;
dVdHpersec=dVdH_Cessna(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*k2mph,hdot*60)
title('Cessna 182 Unpowered glide analysis')
xlabel('Vi mph')
ylabel('hdot ft/min')

tanGbg=16000/(30*5280); % tangent of best glide angle
Gbg=-atan(tanGbg);


Vibg=(80)*0.86897; % Indicated knots
Vc=AS1Cessna182(Vibg,0); % Calibrated knots
Hc=Hipress(1); % 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*k2mph,hdotPOH*60,'rx'); 
text2(.02,.15,str1)
text2(.02,.1,'The red x is from dropping 16000 ft in 30 sm at Vi=80 mph.')
text2(.02,.05,'The red x does not consider the dV/dH term from gliding at constant Vi (not constant Vt).')
grid on; hold off



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

tanGbg=16000/(30*5280); % tangent of best glide angle
Gbg=-atan(tanGbg);

Vibg=(80)*0.86897; % Indicated knots
Vc=AS1Cessna182(Vibg,0); % Calibrated knots

Hc=Hipress(1); % 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,.85,str1); 
text2(.02,.1,'The red x is from dropping 16000 ft in 30 sm at Vi=80 mph.')
text2(.02,.05,'The red x does not consider the dV/dH term from gliding at constant Vi (not constant Vt).')
grid on; hold off

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

tanGbg=16000/(30*5280); % tangent of best glide angle
Gbg=-atan(tanGbg);

Vibg=(80)*0.86897; % Indicated knots
Vc=AS1Cessna182(Vibg,0); % Calibrated knots
Hc=Hipress(1); % 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*k2mph,Gbg*57.3,'rx'); text2(.05,.85,str1); 
text2(.02,.1,'The red x is from dropping 16000 ft in 30 sm at Vi=80 mph.')
text2(.02,.05,'The red x does not consider the dV/dH term from gliding at constant Vi (not constant Vt).')
grid on; hold off

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

disp(' '); disp('Copy the above data into a text file called Cessna_Glide_Data.txt.')
disp('This file is the input file for Step_1Alternate_Glide_Cessna.m.')
format loose