function [BHP,thrust]=Sthrust4(Vtrueknots,Hpress,RPM,Temp,Nblades,BHPmax,Din,C,m,b)
% function [BHP,thrust]=Sthrust4(Vtrueknots,Hpress,RPM,Temp,Nblades,BHPmax,Din,C,m,b)
% Function to compute PER ENGINE brake horsepower and thrust.
%
% Vtrueknots, Hpress, RPM, Temp, BHP, thrust can be vectors of the same size
% Vtrueknots, knots 
% Hpress, ft 
% RPM is engine revolutions per minute of the propeller
% Temp is the Outside air temperature (degR) of the engine.
% BHP is brake horsepower PER ENGINE. <<<<<<<<<<
% thrust in lbf. per engine. 
% Total thrust for two engines is twice the returned value.
% This version uses the altitude drop-off factor with C=.12
% and either 2-bladed or 3-bladed generic propellers.
% BHPmax=max BHP per engine
% Typically BHPmax=300 hp
% Din in propeller diameter in inches
%      Typically for Nblades=2 for 2-bladed prop, Din=78 inches)
%                for Nblades=3 for 3-bladed prop, Din=76 inches)
% C in an input (typically C=.12)
% m is the slope of the prop[eller polar
% b is the intercept of the propeller polar.
% Typical numbers
%   if Nblades==2
%     m=1.34710448406077; b=-0.0424066524285002;  % 2 bladed propeller polar 2700 RPM
%     %m=1.39222; b=-.043472;  % 2 bladed propeller polar 2750RPM
%   elseif Nblades==3
%     m=1.76210205056895; b=-0.0884525793118286;  % 3 bladed propeller polar  2700RPM
%   else
%     m=1.34710448406077; b=-0.0424066524285002;  % default for other blade numbers
%   end
	
% 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*Temp);
sigma=rho/0.00237691267925741;
%C=.12; % C in an input
phi=(sigma-C)/(1-C);  % altitude drop-off factor
BHP=BHPmax*phi; % Compute hp from altitude drop-off factor

Dft=Din/12; 
Vtruefps=Vtrueknots*1.687808;  % Convert knots to ft/sec

n=RPM/60; 
Cpi=550*BHP./(rho.*n.^3*Dft^5); 
J=Vtruefps./(n*Dft); 

CtoverJ2=m*(Cpi./(J.*J))+b;    % Use propeller polar
Cti=CtoverJ2.*J.*J;
thrust=Cti.*rho.*n.^2*Dft^4;



