% MATLAB Script to study the Propeller Polar for Baron 58
% This script uses Gthrust.m. Gthrust.m uses Gengine.m and
% the generic propeller data.

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

% Measured quantities
Vindicated1=60:1:89; %knots
Vindicated2=90:5:180; %knots
Vindicated=[Vindicated1,Vindicated2];

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.
nblades=2;
Din=78
Dft=Din/12;
RPM=2750
n=RPM/60
for i=1:length(Vtrueknots)
	[BHP(i),thrust(i),Mtip(i)]=Gthrust(Vtrueknots(i),Hvector(i),RPM,temp(i),nblades);
end
J=Vtruefps./(n*Dft);
%Ct=
Cpi=550*BHP./(rho.*n.^3*Dft^5);
%Ct/J*J
%Cp/J*J
%plot(Cp/J*J,Ct/J*J)
%thrust=Cti.*rho.*n.^2*Dft^4
Cti=thrust./(rho.*n.^2*Dft^4);

figure(1)
x=Cpi./(J.*J);
y=Cti./(J.*J);
subplot(211)
plot(Vindicated,y)
xlabel('Vindicated (knots)')
ylabel('Ct / J^2')
title('Comparison of Hartzel Chart-based Propeller model with Lowery''s Propeller Polar')
subplot(212)
plot(x,y)
string1=['Propeller Polar: ',num2str(nblades),' bladed prop, ',num2str(Din),' inches in diameter'];
title(string1)
xlabel('Cp / J^2')
ylabel('Ct / J^2')
%function [m,b]=straightline(x,y)
[m,b]=straightline(Cpi./(J.*J),Cti./(J.*J))
hold on
plot(Cpi./(J.*J),m*(Cpi./(J.*J))+b,'x')
hold off
string2='From Charts';
string3='From Prop Polar';
string4='Equation of Lowery''s Propeller Polar Ct/J^2 = m*(Cpi/J^2) + b';
string5=['where for this case m= ',num2str(m),' b= ',num2str(b)];
text(.15,.1,string4)
text(.15,.05,string5)
legend(string2,string3)
Cti2=(m*(Cpi./(J.*J))+b).*(J.*J);
thrust2=Cti2.*rho.*n.^2*Dft^4;

figure(2)
subplot(311)
plot(Vindicated,thrust)
%xlabel('Vindicated (knots)')
ylabel('thrust(lbf)')
title('Thrust using prop data')
text(70,350,'thrust=Ct*rho*n^2*Dft^4')
subplot(312)
plot(Vindicated,thrust2)
%xlabel('Vindicated (knots)')
ylabel('thrust(lbf)')
title('Thrust using propeller polar')
subplot(313)
deltat=thrust2-thrust;
plot(Vindicated,(deltat./thrust)*100)
xlabel('Vindicated (knots)')
ylabel('% error in thrust')
title('Thrust modeling error')

figure(3)
subplot(211)
x=J./Cpi.^(1/3);
etai=Cti.*J./Cpi;
plot(x,etai)
xlabel('J/Cp^(1/3)')
ylabel('eta')
title('Propeller efficiency')
hold on
eta=b*x.^3+m*J;
plot(x,eta,'x')
hold off
legend(string2,string3)
text(1.5,.5,'eta= b*(J/Cp^(1/3))^3  + m*J')
text(1.5,.45,string5)
subplot(212)
plot(Vindicated,etai)
ylabel('eta')
xlabel('Vindicated (knots)')

figure(4)
plot(Vindicated,Mtip)
xlabel('Vindicated')
ylabel('Mtip')
title('Tip Mach number')
string7='Vtip1=n*Dft*pi';
string8='Vtip=sqrt(Vtip1*Vtip1+Vtruefps*Vtruefps)';
string9='Mtip=Vtip/afps';
text(100,0.655,string7)
text(100,0.650,string8)
text(100,0.645,string9)

