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

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,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)

% 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

% % Compute thrust from tabulated propeller data and gngine graphical data
% for i=1:length(Hp2)
%     VtrueKnots(i);
%     HcVector(i);
%     RPM;
%     tVector(i);
%     nblades;
%     BHPmax;
%     Din;
% 	[BHP(i),thrust(i),Mtip(i)]=Gthrust2(VtrueKnots(i),HcVector(i),RPM,tVector(i),nblades,BHPmax,Din);
% end

% Compute thrust from propeller polar
% C=.12; m=1.40; b=-0.032 % good agreement with Gthrust2
% C=.12; m=1.38; b=-0.032 % best agreement with drag polar
% C=.12; m=1.38; b=-0.042 % changing b changes CD0
% C=.12; m=1.38; b=-0.032;  % Changing m spreads out the drag data that should be coincident
% C=.12; m=1.76210205056895; b=-0.0884525793118286; % 3 blades, 76 in diameter
% C=.12; m=1.1106; b=-0.025846; % Step 1: This is from a fit of the tabulated prop data
% C=.12; m=1.38; b=-0.025846;   % Step 2: This change gets the data to form one line
%                               and gets the slope k correct. Now we have
%                               m and k
% C=.12; m=1.38; b=-0.032;
% C=.12; m=1.38; b=0.02923;
% C=.12; m=1.38; b=0.01011;
% C=.12; m=1.38; b=0.02445575;
% C=.12; m=1.38; b=-0.042;
% C=.12; m=1.38; b=-0.0433289779767523;
% C=.12; m=1.38; b=-0.0443259813783418;
% C=.12; m=1.38488275847452; b=-0.04577167413532;
% C=.12; m=1.38488275847452; b=-0.0468088203213117;
% C=.12; m=1.37318901443435; b=-0.044775570759613;
C=.12; m=1.3666; b=-0.040327; 


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

% % Plot thrust comparison
% ifig=ifig+1; figure(ifig)
% plot(thrust,'rd'); hold on; plot(thrust2,'gx'); hold off
% xlabel('Data point number');ylabel('One Engine Thrust (lbf)')
% legend('Engine and Prop Tabulated data','Propeller polar')
% str2=['Propeller Polar: Cp/J*J= ',num2str(m),'*Cp/J*J + ',num2str(b),', C= ',num2str(C)];
% text2(.1,.1,str2)
% 
% ifig=ifig+1; figure(ifig)
% plot(BHP,'rd'); hold on; plot(BHP2,'gx'); hold off
% legend('Engine Tabulated data','Gagg-Ferrar rule (function PHI)')
% str2=['Propeller Polar: Cp/J*J= ',num2str(m),'*Cp/J*J + ',num2str(b),', C= ',num2str(C)];
% text2(.1,.1,str2)
% xlabel('Data point number');ylabel('One Engine BHP (hp)')

% 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)-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);
    Drag=Thrust2engine.*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 with CD= ',num2str(CD0temp),' + ',num2str(ktemp),'*CL^2'];
text2(.15,.5,strx)



CD0=.02538   % Paracite drag coefficient
e=.66556      % Airplane efficiency factor

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 from 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)
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')



