% 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
clear
disp(' '); disp('Start Here'); disp(' ');
load Baron_Glide_Data.txt  % Pilot's Operating Handbook page 5-30
Hp2=Baron_Glide_Data(:,1)   % Pressure altitude (ft)
OAT2=Baron_Glide_Data(:,2)  % degC
hdot2=Baron_Glide_Data(:,3) % rate of climb (ft/min)
Wt=Baron_Glide_Data(:,4)    % Weight (lbf)
ViVector=Baron_Glide_Data(:,5)    % Vindicated knots

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
VcVector=AS1Baron58(ViVector,flapVector); % knots
HcVector=ALTBaron58(Hp2,ViVector,flapVector); % 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=9*OAT2/5+32+459.67;
rho2=pVector./(R*TempDegR) % Non-standard atmosphere density. slugs/ft^3
rho                        % Use standatd atmosphere density. slugs/ft^3
%rho=rho2                   % Alternatively, use rho2. 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)

S=199.2;    % wing area (ft^2)
B=37+10/12; % Wing span (ft)
AR=B*B/S    % Aspect ratio


% 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(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)

%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;
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'];
text2(.15,.2,str1)
str2=['and comes from a least squares fit to the green data.'];
text2(.15,.15,str2)





