% 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
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=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=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)')


%pause

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.


% 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
Thrust1engine=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);
    CD=CD0+k*CL.*CL;
    Drag=CD.*(.5*rho.*VtFtperSec.*VtFtperSec*S); 
    Thrust1engine=(1./cos(alpharad)).*(Drag+Wt.*singamma.*multiplier);
    alpharad=(CL-Cl0)/Clalpharad;
    alphadeg=alpharad*57.3;
end
% OK up to here

[Ct,CtoverJ2,Cp,CpoverJ2]=ThrustCoef2(VtrueKnots,HcVector,RPM,TempDegR,Thrust1engine,BHPmax,Din,C);

% function [m,b]=straightline(x,y)
disp(' '); disp('Propeller polar determined from Sawtooth Climb/Glide data')
[m,b]=straightline(CpoverJ2',CtoverJ2')


ifig=ifig+1; figure(ifig)
plot(CpoverJ2,CtoverJ2,'x')
xlabel('Cp/J^2')
ylabel('Ct/J^2')
title('Propeller polar determined from Sawtooth Climb/Glide data')

x=.1:.01:.24;
y=m*x+b;
hold on; plot(x,y,'r-'); 
str1=['Red line is the propeller polar with m= ',num2str(m),' and b= ',num2str(b),'.'];
text2(.1,.1,str1)
str2=['This is based on a drag polar with k= ',num2str(k),', CD0= ',num2str(CD0), ' and e= ',num2str(e),'.'];
text2(.1,.05,str2)
hold off


%[VMeKNOTS,VmeKNOTS,VyeKNOTS,hdotmaxFPM,VxeKNOTS,gammaxDEG,VbgeKNOTS,gammabgDEG,VmdeKNOTS,hdotmdFPM]=...
% boot2(S,A,P0HP,C,d,CD0,e,m,b,RPM0,W0,sigma,W,nengine)
% INPUTS AND TYPICAL VALUES FOR CESSNA 172 (The Bootstrap Data Plate)
%   S=174       Wing reference area (ft^2)
%   A=7.38      Wing aspect ratio
%   P0HP=160    Rated sea level horsepower per engine, hp
%   C=.12       Altitude dropoff parameter
%   d=6.25      Propeller diameter (ft)
%   CD0=.037    Paracite drag coefficient
%   e=.72       Airplane efficiency factor
%   m=1.7       Propeller polar slope
%   b=-.0564    Propeller polar intercept
%   RPM0=2700   Prop revolutions per min
%   W0=2400     Standard weight lbf
%   sigma=.786  Density ratio on given test day
%   W=1800      Test day weight, lbf
%   nengine     Number of engines

disp(' '); disp(' ')
sigma=1 % density ratio for sea level
disp(' '); disp('Bootstrap Outputs for sea level performance'); format compact
disp('%   VMeKNOTS    Equivalent velocity for max level flight speed, knots')
disp('%   VmeKNOTS    Equivalent velocity for min level flight speed, knots')
disp('%   VyeKNOTS    Equivalent velocity for best rate of climb, knots')
disp('%   hdotmaxFPM  Best climb rate, ft/min')
disp('%   VxeKNOTS    Equivalent velocity for best climb angle, knots')
disp('%   gammaxDEG   Best climb angle, deg')
disp('%   VbgeKNOTS   Equivalent velocity for best (longest) glide, knots')
disp('%   gammabgDEG  Best (longest) glide angle, deg')
disp('%   VmdeKNOTS   Equivalent velocity for minimum descent rate, knots')
disp('%   hdotmdFPM   Minimum rate of descent, ft/min')
[VMeKNOTS,VmeKNOTS,VyeKNOTS,hdotmaxFPM,VxeKNOTS,gammaxDEG,VbgeKNOTS,gammabgDEG,VmdeKNOTS,hdotmdFPM]=...
   boot2(S,AR,BHPmax,C,Din/12,CD0,e,m,b,RPM,W0,sigma,W0,1)
disp(' '); disp('Change units from knots to mph')
VMemph=VMeKNOTS*k2mph
Vmemph=VmeKNOTS*k2mph
Vyemph=VyeKNOTS*k2mph
Vxemph=VxeKNOTS*k2mph
Vbgemph=VbgeKNOTS*k2mph
Vmdemph=VmdeKNOTS*k2mph


deltaH=30*5280*tan(gammabgDEG/(180/pi));
str2=['This aircraft will drop ',num2str(-deltaH),' ft in 30 sm. at Vbgemph of ',num2str(VbgeKNOTS*k2mph),' mph.'];
str3=['This aircraft is supposed to drop approximately 16000 ft in 30 sm. at Vbgimph of 80 mph.'];
disp(str2); disp(str3)
format loose
