function [BHP,thrust]=Gthrust2(Vtrueknots,Hpress,TempdegR,Nblades,BHP0,RPM0,BAF,Din,Z,istractor,C)
% function [BHP,thrust]=Gthrust2(Vtrueknots,Hpress,TempdegR,Nblades,BHP0,RPM0,BAF,Din,Z,istractor,C)
% Function to compute PER ENGINE maximum brake horsepower and maximum thrust
% FOR CONSTANT SPEED PROPELLER.
% Full throttle operation is assumed. 
%
% INPUTS:
% 	Vtrueknots, Hpress, and TempdegR can be vectors of the same length.
% 	Vtrueknots, knots 
% 	Hpress, geopotential ft
% 	TempdegR is the Outside air temperature (degR) of the engine.
%	Z is the ratio of fuselage diameter to propeller diameter.
%	istractor=1 for tractor aircraft (engine in the front).
%	istractor=0 for pusher aircraft (engine in the rear).
%
% 	BHP0 is the maximum horsepower at sea level for rated RMP given by RPM0.
% 	RPM0 is engine revolutions per minute at rated power (a scalar)  
% 	BAF is the blade activity factor.
% 	Din is the diameter of the propeller in inches.
%	C is the Power dropoff parameter (default=.12)
% OUTPUTS:
%   BHP and thrust will be vectors if Vtrueknots, Hpress, and TemdegRp are vectors.
% 	BHP is brake horsepower PER ENGINE.
% 	thrust in lbf. per engine. 
% 	Total thrust for two engines is twice the returned value.
%
% Notes:
% 	The Gagg and Farrar factor is used to predict maximum engine 
%   performance as a function of altitude (see Lowry, page 141).
% 	Full throttle operation of the engine is assumed. 
%	A constant speed propeller (i.e., with a governor) is assumed.
%
% Reference: John T. Lowry, Performance of Light Aircraft, AIAA, 1999, pp. 176-182.
% Secondary Reference: Perkins and Hage, Airplane Perfformance, Stability and Control,
%                      Wiley, 1949, pp. 149-153.

Dft=Din/12; 
Vtruefps=Vtrueknots*1.687808;  % Convert knots to ft/sec	  % a vector
if (isempty(C)); C=.12; end

TAF=Nblades*BAF;   % a scalar, eqn. 6.55, page 178.
X=0.001515*TAF-0.0880;  % a scalar, eqn 6.57 page 180.

% Get atmospheric pressure for the given pressure altitude.
[temp,press,rho,Hgeopvector]=atmosphere4(Hpress,0);
% Compute density for given temperature (Temp) and pressure (press).
rho=press./(1716.55*TempdegR);  % a vector
n=RPM0/60;  % a scalar
% C=Power dropoff parameter (typically 0.12)
% sigma is the density ratio
sigma=rho/0.00237691267925741;  % a vector
phi=(sigma-C)/(1-C); % a vector
BHP=BHP0*phi;	     % a vector
Cp=550*BHP./(rho*n^3*Dft^5);  % a vector
Cpx=Cp/X;            % a vector
J=Vtruefps/(n*Dft);  % a vector
x=J./Cp.^(1/3);      % a vector

load fig619DataFile  % load data from Figure 6.19 of Lowry
load xDataFile       % load data from Figure 6.19 of Lowry
load yDataFile       % load data from Figure 6.19 of Lowry
% Get propeller efficiency from Figure 6.19 of Lowry

kk=find(Cpx<.15); 
if isempty(kk)
else
	disp('Cpx is off the scale of figure 6.19')
	Cpx
	Cpx(kk)=.15;
	disp('Cpx has been modified as follows')
	Cpx
end	

kk=find(Cpx>1.4); 
if isempty(kk)
else
	disp('Cpx is off the scale of figure 6.19')
	Cpx
	Cpx(kk)=1.4;
	disp('Cpx has been modified as follows')
	Cpx
end	

eta= interp2(xDataFile,yDataFile,fig619DataFile,x,Cpx);   % a vector
Ct=eta.*Cp./J;
if (istractor)
	SDEF=1.05263-0.00722*Z-0.16462*Z^2-0.18341*Z^3; % Tractor: Eqn 6.58 page 181
else
	SDEF=1.05263-0.04185*Z-0.01481*Z^2-0.62001*Z^3;  % Pusher: Eqn 6.59 page 18
end
thrust=SDEF*Ct.*rho*n^2*Dft^4; % Eqn 6.60 page 182



