% Script to determine drag polar from rate of GLIDE data.
% Glide data is arranged in columns of a test file
% Hp(ft)   OAT (degC)   Hdot(ft/min) Wt(lbf) Vi(knots)

close all; echo off
clear
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'); disp(' ');
load Cessna_Glide_Data.txt         % Pilot's Operating Handbook page 5-30
Hp2=     Cessna_Glide_Data(:,1);   % Pressure altitude (ft)
OAT2=    Cessna_Glide_Data(:,2);   % degF      Note units <<<<<<<<<<<<<<
hdot2=   Cessna_Glide_Data(:,3);   % rate of climb (ft/min)
Wt=      Cessna_Glide_Data(:,4);   % Weight (lbf)
ViVector=Cessna_Glide_Data(:,5);   % Vindicated knots


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 (from Step 1)
%e=Unknown;     % Airplane efficiency factor from Step 1
%m=Unknown;    % Slope of propeller polar (from Step 2)
%b=Unknown; % Intercept of propeller polar (from Step 2)

% Auxilary data
W0=2950;       % Standard weight and test day weight, 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;


ifig=0;
index=1:length(Hp2);
ifig=ifig+1; figure(ifig)
subplot(321);plot(index,Hp2,'x')
title('Baron 58 2-Engine Glide Data')
ylabel('pressure altitude (ft)')
subplot(322);plot(index,OAT2,'x')
ylabel('OAT (degC)')
subplot(323); plot(index,hdot2,'x')
ylabel('climb rate (ft/min)')
subplot(324); plot(index,Wt,'x')
ylabel('Weight (lbf)')
subplot(325); plot(index,ViVector,'x')
ylabel('Vi (knots)')

flapVector=0*ones(size(Hp2)); % No Flaps
% function VcalibratedKnots=AS1Cessna182(VindicatedKnots,flap)
VcVector=AS1Cessna182(ViVector,flapVector); % knots
HcVector=Hp2; % a vector same length as Hp2
[tVector,pVector,rho,Hgeopvector]=atmosphere4(HcVector,0); % altitude is geopotential, p is a vector


R=1716.55;	%ft^2/(sec^2degR)
TempDegR=OAT2+459.67;  % Convert from degF to degR
rho=pVector./(R*TempDegR); % Non-standard atmosphere density. slugs/ft^3

[VeVector,qcOverp0,qcOverp]=AS2(VcVector,pVector); % knots
[VtrueKnots,M,aknots,rho3,sigma]=AS3(VeVector,tVector,HcVector);
VtFtperSec=VtrueKnots*1.687808; % Convert knots to ft/sec
HdotFtperSec=hdot2/60; % convet to ft/sec
singamma=HdotFtperSec./VtFtperSec; % sine(gamma)
gamma=asin(singamma); % gamma in radians
cosgamma=cos(gamma); % cos(gamma)



% Assume alpha is non-zero, and iterate to find it.
Cl0=.353; Clalpharad=5.62;  % for PA-28
% Iterate on alpha and gamma
alpharad=zeros(size(Hp2)); % initialize and get matrix sizes correct
singamma=(hdot2/60)./VtFtperSec;
gammarad=asin(singamma);       
for i=1:20
    Lift=Wt.*cos(gammarad);
    CL=2*Lift./(rho.*VtFtperSec.*VtFtperSec*S);
    dVdHpersec=dVdH_Cessna(ViVector,flapVector,Hp2);
    % See equation 12.11 page 122 of Kimberlin.
    multiplier=(1+VtFtperSec.*dVdHpersec/32.17);
    Drag=-Wt.*singamma.*multiplier;    
    CD=Drag./(.5*rho.*VtFtperSec.*VtFtperSec*S);
    alpharad=(CL-Cl0)/Clalpharad;
    alphadeg=alpharad*57.3;
end

ifig=ifig+1; figure(ifig)
plot(CL.*CL,CD,'go')
xlabel('CL^2')
ylabel('CD')
title('Drag Polar, GLIDE')
str2=['Green o is from climb data'];
text2(.45,.35,str2)

disp(' ');
disp('Drag polar determined from Sawtooth Glide data')
%function [m,b]=straightline(x,y)
[kGLIDE,CD0GLIDE]=straightline((CL.*CL)',CD')
eGLIDE=1/(pi*AR*kGLIDE)
CLGLIDE=CL;
CDGLIDE=CD;

CL2=0:.1:1.4;
CD2=CD0GLIDE+kGLIDE*CL2.*CL2;
hold on; plot(CL2.*CL2,CD2); hold off
str1=['This line represents a drag polar with CD= ',num2str(CD0GLIDE),' + ',num2str(kGLIDE),'*CL^2 (e=',num2str(eGLIDE),').'];
text2(.15,.2,str1)
str2=['and comes from a least squares fit to the green data.'];
text2(.15,.15,str2)





