% This script assumes you have a drag polar and propeller polar and 
% compares the maximum speed from the Piper Warrior 3 Pilot's Info Manual 
% with the max speed computed from the drag and propeller polars.
% Note: The speed data with hdot=0 can also be incorporated into 
% Step 2 and processed with the climb data. In this case the
% speed data will be included in the least squares for for
% the propeller polar as part of Step 2.

close all
clear
disp(' '); disp('Start Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<'); disp(' ');
load  Warrior3SPEEDdata.txt       % Pilot's Information Manual, page 5-20b
Hp2=  Warrior3SPEEDdata(:,1);     % Indicated Pressure altitude (ft)
OAT2= Warrior3SPEEDdata(:,2);     % degF        <<<<< Note units
Wt=   Warrior3SPEEDdata(:,3);     % Weight (lbf)
VtrueKnots=Warrior3SPEEDdata(:,4);     % Maximum True 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=0;
CD0climb=CD0glide-DeltaCD0windmillingprop; % Paracite drag for climb is less
e=.6598;       % Airplane efficiency factor
m=1.6414;    % Slope of propeller polar (from Step 2)
b=-.049476;    % Intercept of propeller polar (from Step 2)


% CD0glide=.03727;    % Paracite drag coefficient including drag of windmilling prop
% DeltaCD0windmillingprop=.01152; % Computed by method in Roskam
% CD0climb=CD0glide-DeltaCD0windmillingprop; % Paracite drag for climb is less
% e=.6598;       % Airplane efficiency factor
% m=1.6087;    % Slope of propeller polar (from Step 2)
% b=-.069832;    % 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;
ifig=ifig+1; figure(ifig)
plot(VtrueKnots,Hp2)
xlabel('TRUE AIRSPEED-KTS'); ylabel('PRESSURE ALTITUDE-FEET')
title('Maximum True Airspeed in Level Flight for Piper Warrior 3')


%pause
[Temp,pVector,rho,Hgeopvector]=atmosphere4(Hp2,0); % altitude is geopotential, p is a vector
TempDegR=OAT2+459.67;
VtFtperSec=VtrueKnots*1.687808; % Convert knots to ft/sec

n=length(Hp2);
% set up for a local loop on velocity
Vkts=100:1:140;
Vfps=Vkts*1.687808; % Convert knots to ft/sec
Vmat=zeros(size(Hp2));

for i=1:n ; % For each altitude
    % Power required for a given altitude Hp2(i)
    CL=2*Wt(i)./(rho(i).*Vfps.*Vfps*S); % Compute for many velocities
    CD=CD0climb+k*CL.*CL;
    Drag=CD.*(.5*rho(i).*Vfps.*Vfps*S); 
    PowerRequired=Drag.*Vfps;
    % Power Available for a given altitude
    T=Thrust(Vkts,Hp2(i),RPM,TempDegR(i),BHPmax,Din,C,m,b);
    PowerAvailable=T.*Vfps;
    ExcessPower=PowerAvailable-PowerRequired;
    %YI = INTERP1(X,Y,XI)
    VmaxtrueKnots=interp1(ExcessPower,Vkts,0); % Interpolate to the velocity with no excess power
    Vmat(i)=VmaxtrueKnots; % Store result for each altitude
end
hold on; plot(Vmat,Hp2,'rx'); 
%axis([90 120 0 10000])
text2(.1,.1,'Blue line is from the Warrior 3 Pilot''s Information Manual, page 5-20b')
text2(.1,.05,'Red x is computed from the drag and propeller polars')
text2(.1,.9,'Full Throttle, Best Power Mixture')
hold off; % Add to existing plot for comparison

Vmat
disp(['CD0glide=',num2str(CD0glide),'  e=',num2str(e)])
disp(['CD0climb=',num2str(CD0climb)])
disp(['m=',num2str(m),'  b=',num2str(b)])