% MATLAB Script to study the climb performance for Baron 58
% Using a drag polar determined by Experiment #4
% This script uses Gthrust.m. Gthrust.m uses Gengine.m and
% the generic propeller data.

clear all
disp(' ')
disp('Start here')
disp(' ')
disp('Generic 6 cylinder 300 HP aircraft engine')
disp('Generic propeller')

% Measured quantities
Vindicated=60:5:180; %knots
Hvector=1000*ones(size(Vindicated));  %Hvector is Geopotential altitude (ft)
flap=0;

% Computed quanities
GeometricFlag=0; %Hvector is Geopotential altitude
[temp,p,rho,Hgeopvector]=      atmosphere4(Hvector,GeometricFlag);
Vcalibrated=                   AS1Baron58(Vindicated,flap);
[Vequivalent,qcOverp0,qcOverp]=AS2(Vcalibrated,p); %Knots
[Vtrueknots,M,aknots,rho,sigma]=AS3(Vequivalent,temp,Hvector);
Vtruefps=Vtrueknots*1.687808; % Convert knots to ft/sec

% Get engine performance
% [BHP,thrust,Mtip]=Gthrust(Vtrueknots,Hpress,RPM,Temp,Nblades)
BHP=zeros(size(Vtrueknots)); % 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.
for i=1:length(Vtrueknots)
	[BHP(i),thrust(i),Mtip(i)]=Gthrust(Vtrueknots(i),Hvector(i),2450,temp(i),2);
end
W=5500  %lbf
S=199.2 %ft^2

qbar=.5*rho.*Vtruefps.*Vtruefps; % Dynamic pressure
CL=W./(qbar*S); % Lift coefficient

% Drag Polar
%k= 0.052548;   CD0= 0.019626;  %(Piper Twin, PA-30, Tc='44, flap=0, over 9 points)
%k= 0.049698;   CD0= 0.020861   %(Piper Twin, PA-30, Tc='44, flap=0, over 8 points)
%k= 0.044651;   CD0= 0.022674   %(Piper Twin, PA-30, Tc='44, flap=0, over 7 points)
%k= 0.08;   CD0= 0.03;  %Guess

% CD=CD0+CL.*CL*k
% Cd= 0.026232 + 0.13927 *CL^2
%k=0.13927
%k=.1
k=.08
CD0=0.025
CD=CD0+CL.*CL*k;

% Power required
tstring=['Altitude= ',num2str(Hvector(1))];
PR=S*qbar.*CD.*Vtruefps; % ft-lbf/sec

% Power Available
% PA=thrust.*Vtruefps % ft-lbf/sec one engine
PA=2*thrust.*Vtruefps; % ft-lbf/sec two engines


% Plot power required and Power Available
figure(1)
clf
plot(Vindicated,PR,'x',Vindicated,PA,'o')
text2(.1,0.02,tstring)
xlabel('Indicated airspeed, knots')
ylabel('Power, ft-lbf/sec')
legend('Power Required','Power Available (2 engine)')
title(['Baron Power Required and Power Available for ',tstring])
z=axis;
axis([z(1)         z(2)       0      250000])

% Climb performance equation
hdot=(PA-PR)/W;

% Plot climb performance
figure(2); 
clf
plot(Vindicated,hdot*60)
xlabel('Indicated airspeed, knots')
ylabel('dh/dt, ft/min')
axis([z(1)         z(2)       0      2400])
title(['Baron Climb Performance for ',tstring])
text2(0.1,0.1,tstring)
grid on

