% Script to determine drag polar from rate of climb data.
% For the Warrior 3 we can use rate of climb data from the 
% Pilot's Information Manual, p. 5-19
close all
clear
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'); disp(' ');
load  Warrior3RCdata.txt       % Pilot's Information Manual, page 5-19
Hp2=  Warrior3RCdata(:,1);     % Indicated Pressure altitude (ft)
OAT2= Warrior3RCdata(:,2);     % degF        <<<<< Note units
hdot2=Warrior3RCdata(:,3);     % rate of climb (ft/min)
Wt=   Warrior3RCdata(:,4);     % Weight (lbf)
ViKTS=Warrior3RCdata(:,5);     % Indicated airspeed, knots



echo on; format compact
% Bootstrap Data Plate: Warrior 3
S=170;         % wing area (ft^2)
B=35;          % Wing span (ft)
BHPmax=160;    % Rated engine horsepower (hp), See Owner's Manual p. 4-2
RPM=2700;      % Rated engine speed, RPM
Din=74;        % Propeller diameter, in
CD0glide=.03727;    % Paracite drag coefficient including drag of windmilling prop
%DeltaCD0windmillingprop=.01152; % Computed by method in Roskam
DeltaCD0windmillingprop=0;
CD0climb=CD0glide-DeltaCD0windmillingprop; % Paracite drag for climb is less
e=.6598;       % Airplane efficiency factor
%m=Unknown;    % Slope of propeller polar (from Step 2)
%b=Unknown;    % Intercept of propeller polar (from Step 2)

% Auxilary data
W0=2440;       % Standard weight, lbf: Max TO wt normal category
Wbg=2440;      % Weight for which Best glide os known, lbf, page 5-27
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; 


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,ViKTS,'x')
ylabel('Vi (kts)')


%pause

flapVector=0*ones(size(Hp2)); % No Flaps
ViVector=ViKTS; % knots
VcVector=AS1Warrior3(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,TempDegR,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_Warrior3(ViVector,flapVector,Hp2);
    % See equation 12.11 page 122 of Kimberlin.
    multiplier=(1+VtFtperSec.*dVdHpersec/32.17);
    CD=CD0climb+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
RPMmat=RPM*ones(size(Hp2));
%[Ct,CtoverJ2,Cp,CpoverJ2]=ThrustCoef(VtrueKnots,HcVector,RPMmat,TempDegR,Thrust1engine,BHPmax,Din,C);
[Ct,CtoverJ2,Cp,CpoverJ2]=ThrustCoef2(VtrueKnots,HcVector,RPMmat,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),', CD0climb= ',num2str(CD0climb), ' 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)
%   CD0glide=.037    Paracite drag coefficient with prop drag
%   CD0glide=.026    Paracite drag coefficient without drag from
%   windmilling prop
%   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(['CD0glide=',num2str(CD0glide),'  e=',num2str(e)])
disp(['CD0climb=',num2str(CD0climb)])
disp(['m=',num2str(m),'  b=',num2str(b)])
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]=...
   boot3(S,AR,BHPmax,C,Din/12,CD0glide,CD0climb,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=22.75*5280*tan(gammabgDEG/(180/pi));
str2=['This aircraft will drop ',num2str(-deltaH),' ft in 22.75 nm. at VbgeKNOTS of ',num2str(VbgeKNOTS),' knots.'];
str3=['This aircraft is supposed to drop approximately 12000 ft in 22.75m. at Vbgikts of 73 knots.'];
disp(str2); disp(str3)
format loose

% This propeller efficiency data is from Jeppesen page 82 for PA28-180
data=[0.23559	0.43468
0.24841	0.45430
0.26065	0.47391
0.27348	0.49582
0.28806	0.51429
0.30381	0.53507
0.32364	0.55817
0.34231	0.57782
0.36215	0.59633
0.38316	0.61484
0.40359	0.63335
0.42460	0.65071
0.44853	0.67038
0.47421	0.68892
0.50048	0.70632
0.52968	0.72487
0.56121	0.74230
0.59507	0.75973
0.61960	0.77137
0.64938	0.78533
0.67742	0.79584
0.70896	0.80406
0.73817	0.80998
0.76680	0.81359
0.79017	0.81373
0.81823	0.80814
0.83694	0.80250
0.85507	0.79341
0.87144	0.78316
0.88900	0.76832
0.90772	0.75003
0.92470	0.72829
0.94052	0.70310
0.95283	0.67903
0.96397	0.65381
0.97453	0.62283
0.98276	0.59415
0.98923	0.56545
0.99394	0.53789]; % This data is from Jeppesen page 82 for PA28
J=data(:,1); etaP=data(:,2);
ifig=ifig+1; figure(ifig)
plot(J,etaP)
xlabel('Advance Ratio J (nd)'); ylabel('Propeller Efficiency (nd)')
J2=VtFtperSec/((RPM/60)*(Din/12));
% [Ct,CtoverJ2,Cp,CpoverJ2]
etaP2=J2.*Ct./Cp;
hold on; plot(J2,etaP2,'rx')
axis([.2 1 0 1])
text2(.1,.15,'Blue line is from Jeppesen page 83 for PA-28-180 fixed pitch propeller')
text2(.1,.1,'Red x is form the climb data')
text2(.1,.05,'etaP=J*Ct/Cp')

