% Script to determine drag polar from rate of climb data.
% For the Baron 58 we can use rate of climb data from the 
% Pilot's Operating Handbook (POH) page 5-30
% This code assumes you know a propeller polar is given.
% Using this it backs out the drag polar.
% It will overplot another input drag polar.
% This code can be used iteratively to find the propeller polar
% by repeated guessing of m and b of the propeller polar.
% It can also be used to check a drag polar and a propeller polar
% found antother way such as by using
% Step_1Glide to compute the drag polar
% and then using Step_2_Climb to compute the propeller polar.


close all
clear

disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<'); disp(' ');
load  Cessna182RCdata.txt       % Pilot's Operating Handbook page 5-30
Hp2=  Cessna182RCdata(:,1);     % Indicated Pressure altitude (ft)
OAT2= Cessna182RCdata(:,2);     % degF        <<<<< Note units
hdot2=Cessna182RCdata(:,3);     % rate of climb (ft/min)
Wt=   Cessna182RCdata(:,4);     % Weight (lbf)
Vimph=Cessna182RCdata(:,5);     % Indicated airspeed, mph  <<<<Note additional col.

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=.048662; % Paracite drag coefficient (from Step 1)
e=.7984;     % Airplane efficiency factor from Step 1
m=1.7172;    % Slope of propeller polar (from Step 2)
b=-0.040323; % 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=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('Cessna 182 Climb Data From Owner''s Manual page 6-3')
ylabel('pressure altitude (ft)')
subplot(322);plot(index,OAT2,'x')
ylabel('OAT (degF)')
subplot(323); plot(index,hdot2,'x')
ylabel('climb rate (ft/min)')
subplot(324); plot(index,Wt,'x')
ylabel('Weight (lbf)')
subplot(325); plot(index,Vimph,'x')
ylabel('Vi (mph)')

flapVector=0*ones(size(Hp2)); % No Flaps
% Vcalibrated=AS1Baron58(Vindicated,flap)
% function [Vequivalent,qcOverp0,qcOverp]=AS2(Vcalibrated,p)
% function [temp,press,rho,Hgeopvector]=atmosphere4(Hvector,GeometricFlag)
%  function Hcalibrated=ALTBaron58(Hindicated,Vindicated,flap)
% function [VtrueKnots,M,aknots,rho,sigma]=AS3(Vequivalent,TdegR,Hpcalibrated)

ViVector=Vimph*mph2k; % knots
VcVector=AS1Cessna182(ViVector,flapVector); % knots
HcVector=Hp2; % a vector same length as Hp2
[tVector,pVector,rhoX,Hgeopvector]=atmosphere4(HcVector,0); % altitude is geopotential, p is a vector


R=1716.55;	%ft^2/(sec^2degR)
TempDegR=OAT2+459.67;   %degR computed from degF
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)

% Get per engine performance
BHP=zeros(size(Hp2)); % Set up matrices of the correct size.
thrust=BHP;                  % Set up matrices of the correct size.
Mtip=BHP;                    % Set up matrices of the correct size.




[BHP2,thrust2]=Sthrust4(VtrueKnots,HcVector,RPM,tVector,nblades,BHPmax,Din,C,m,b);
Thrust1engine=thrust2; % Use thrust from propeller polar


% 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)-Thrust1engine.*sin(alpharad);
    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=Thrust1engine.*cos(alpharad)-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,'rx')
xlabel('CL^2')
ylabel('CD')
title('Drag Polar, climb')
%  function [m,b]=straightline(x,y)
[ktemp,CD0temp]=straightline(CL'.*CL',CD');
CL2temp=0:.1:1;
CD2temp=CD0temp+ktemp*CL2temp.*CL2temp;
hold on; plot(CL2temp.*CL2temp,CD2temp,'r-'); hold off
strx=['Red Line is from drag polar from red x where with CD= ',num2str(CD0temp),' + ',num2str(ktemp),'*CL^2'];
text2(.15,.5,strx)



CL2=0:.1:1;
k=1/(pi*AR*e);
CD2=CD0+k*CL2.*CL2;
hold on; plot(CL2.*CL2,CD2); hold off
str1=['Blue Line is the assumed drag polar with CD= ',num2str(CD0),' + ',num2str(k),'*CL^2'];
str2=['Propeller Polar: Cp/J*J= ',num2str(m),'*Cp/J*J + ',num2str(b)];
text2(.15,.1,str2)
text2(.15,.2,str1)
text2(.05,.9,'If the propeller polar is known, we can compute the drag polar using the red x data.')
text2(.05,.85,'Unfortunately, from this method we cannot separate CD0 from b of the propeller polar.')
text2(.05,.8,'We must use glide data to determine CD0 and perhaps k.')
text2(.05,.75,'This plot can be used to determinme m of the propeller polar and k of the drag polar.')
text2(.05,.70,'Adjust m until the CD data lies on a single straight line.')


% Adjustment rule for making the climb data fall on top of the drag oplar data
f=(S/(4*Din/12)^2);
disp('You might consider changing b to the following value')
bnew=f*(CD0-CD0temp)+b 

qbar1=.5*0.00237691267925741*VeVector.*VeVector*1.687808*1.687808;
qbar2=.5*rho.*VtFtperSec.*VtFtperSec; % same as qbar1

ifig=ifig+1; figure(ifig)
subplot(411)
plot(index,rho,'x')
ylabel('density slug/ft^3')

subplot(412)
plot(index,qbar1,'x')
ylabel('qbar1')

subplot(413)
plot(index,CD,'rx')
hold on; plot(CD0+k*CL.*CL,'gx'); hold off
legend('POH','Polar')
ylabel('CD')

subplot(414)
plot(index,CL,'x')
ylabel('CL')

format loose

