function T=Thrust(Vtrueknots,Hpress,RPM,TempDegR,BHPmax,Din,C,m,b)
% function T=Thrust(Vtrueknots,Hpress,RPM,TempDegR,BHPmax,Din,C,m,b)
% Compute maximum thrust for a given pressure altitude
% and outside air temperature. m and b are from the propeller polar.
% INPUTS
% Vtrueknots, Hpress, RPM, TempDegR can be vectors of the same size
% Vtrueknots, knots 
% Hpress, ft 
% RPM is engine revolutions per minute of the propeller
% TempDegR is the Outside air temperature (degR) of the engine. 
% BHPmax=max BHP per engine
% Typically BHPmax=300 hp
% Din in propeller diameter in inches
% C in an input (typically C=.12)
% OUTPUTS
% T is thrust in pounds
% Notes
% This version Gagg-Farrar Rule for engine horsepower

	
% Get atmospheric pressure for the given pressure altitude.
[tempDegR,press,rho,Hgeopvector]=atmosphere4(Hpress,0);
% Compute density for given temperature (Temp) and pressure (press).
% rho=press./(1716.55*Temp);     % <<<<<<<<<<<<< OLD METHOD
sigma=rho/0.00237691267925741;
%C=.12; % C in an input
phi=(sigma-C)/(1-C);  % altitude drop-off factor
phiT=sqrt(tempDegR./TempDegR);
BHP=BHPmax*phi.*phiT; % Compute hp from altitude drop-off factor (Gagg-Farrar Rule)
% The above method uses the sigma from the standard atmosphere (as a
% function of Hp)
% and does a separate temperature correction
% This was not what I understood to be
% the correct interpretation of the Gagg-Farrar Rule, but it certainly
% appears to work better then my original interpretation which
% involved computing rho as a funtion of actual temperature and pressure.
Dft=Din/12; 
Vtruefps=Vtrueknots*1.687808;  % Convert knots to ft/sec
n=RPM/60; % Convert RPM to rev per sec
Cp=550*BHP./(rho.*n.^3*Dft^5);  % Find power coef.
J=Vtruefps./(n*Dft); % Find advance ratio
CpoverJ2=Cp./(J.*J);  % Cp/J^2
CtoverJ2=m*CpoverJ2+b;
Ct=CtoverJ2.*J.*J;
T=rho.*n.^2*Dft^4.*Ct; % Thrust, lbf