% 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 Baron58RCdata2a.txt  % Pilot's Operating Handbook page 5-30
Hp2=Baron58RCdata2a(:,1)   % Indicated Pressure altitude (ft)
OAT2=Baron58RCdata2a(:,2)  % degC
hdot2=Baron58RCdata2a(:,3) % rate of climb (ft/min)
Wt=Baron58RCdata2a(:,4)    % Weight (lbf)
C=.12; % altitude drop-off parameter for Gagg-Farrar Rule

ifig=0;
index=1:length(Hp2);
ifig=ifig+1; figure(ifig)
subplot(411);plot(index,Hp2,'x')
title('Baron 58 2-Engine Climb Data From POH page 5-30, 105 knots')
ylabel('pressure altitude (ft)')
subplot(412);plot(index,OAT2,'x')
ylabel('OAT (degC)')
subplot(413); plot(index,hdot2,'x')
ylabel('climb rate (ft/min)')
subplot(414); plot(index,Wt,'x')
ylabel('Weight (lbf)')

%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)
Vi=105; % knots
ViVector=Vi*ones(size(Hp2)); % knots
VcVector=AS1Baron58(ViVector,flapVector); % knots
HcVector=ALTBaron58(Hp2,ViVector,flapVector); % 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=9*OAT2/5+32+459.67;
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.
BHPmax=300; % See POH page 1-11
nblades=3;  % See POH page 1-11
Din=76;     % We will stick to the diameter for which we have prop data for.
RPM=2700;
S=199.2;    % wing area (ft^2)
B=37+10/12; % Wing span (ft)
AR=B*B/S    % Aspect ratio

% Assume we know the drag polar form glide analysis
CD0=.0242   % Paracite drag coefficient 
e=.70      % Airplane efficiency factor
CD0=.02538   % Paracite drag coefficient 
e=.66556      % Airplane efficiency factor

k=1/(pi*AR*e)

% 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
Thrust2engine=zeros(size(Hp2)); % initialize and get matrix sizes correct
singamma=(hdot2/60)./VtFtperSec;
gammarad=asin(singamma);       
for i=1:20
    Lift=Wt.*cos(gammarad)-Thrust2engine.*sin(alpharad);
    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);
    CD=CD0+k*CL.*CL;
    Drag=CD.*(.5*rho.*VtFtperSec.*VtFtperSec*S); 
    Thrust2engine=(1./cos(alpharad)).*(Drag+Wt.*singamma.*multiplier);
    alpharad=(CL-Cl0)/Clalpharad;
    alphadeg=alpharad*57.3;
end
% OK up to here

Thrust1engine=Thrust2engine/2;
[Ct,CtoverJ2,Cp,CpoverJ2]=ThrustCoef(VtrueKnots,HcVector,RPM,TempDegR,Thrust1engine,BHPmax,Din,C);

% function [m,b]=straightline(x,y)
[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 climb data, Baron 58 POH')

x=.1:.01:.2;
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
W0=Wt(1)  % Standard weight and test day weight
sigma=rho(1)/.002376913 % density ratio
disp(' '); disp('Bootstrap Outputs'); format compact

[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,2)
deltaH=2*6076.1*tan(gammabgDEG/(180/pi));
str2=['This aircraft will drop ',num2str(-deltaH),' ft in 2 nm. at VbgeKNOTS of ',num2str(VbgeKNOTS),' knots.'];
str3=['This aircraft is supposed to drop approximately 1000 ft in 2 nm. at VbgiKNOTS of 115 knots.'];
disp(str2); disp(str3)
format loose
