function [BHP,thrust,Mtip]=Gthrust(Vtrueknots,Hpress,RPM,Temp,Nblades)
% function [BHP,thrust,Mtip]=Gthrust(Vtrueknots,Hpress,RPM,Temp,Nblades)
% Function to compute PER ENGINE brake horsepower and thrust.
%
% Vtrueknots, Hpress, RPM, Temp, BHP, thrust and Mtip are scalars.
% Vtrueknots, knots 
% Hpress, ft 
% 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 a generic six cylinder 300HP aircraft engine
% and either 2-bladed or 3-bladed generic propellers.
% The engine is assumed to be operating at maximum manifold pressure (MaxMP=1).
% (BHP=Gengine(Trpm,Tmp,Talt,Temp,MaxMP,ploton).
% RPM is engine revolutions per minute 
% Mtip is the mach number of the blade tip.
% Nblades=2 for 2-bladed prop, (Din=78 inches)
% Nblades=3 for 3-bladed prop, (Din=76 inches)

if length(Vtrueknots)>1
	disp('ERROR in Gthrust. Vtrueknots must be a scalar.')
	BHP=NaN;
	thrust=NaN;
	Mtip=NaN;
else
	
	BHP=Gengine(RPM,0,Hpress,Temp,1,0); % Use maximum manifold pressure and do not plot.
	
	% 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);
	
	n=RPM/60; 
	if Nblades==2
		Din=78; % Inches in diameter of 2-bladed propeller
	else
		Din=76; % Inches in diameter of 3-bladed propeller
	end
	
	Dft=Din/12; 
	Vtruefps=Vtrueknots*1.687808;  % Convert knots to ft/sec
	
	Cpi=550*BHP./(rho.*n.^3*Dft^5);
	Ji=Vtruefps./(n*Dft); 
	
	Vtip1=n*Dft*pi;
	Vtip=sqrt(Vtip1.*Vtip1+Vtruefps.*Vtruefps);
	aknots=sound(temp);
	afps=aknots*1.687808;
	Mtip=Vtip./afps;
	
	Cti=zeros(size(RPM));
	CtiHigh=Cti;
	CtiLow=Cti;
	
	
	if Nblades==2                       % For 2-bladed prop
		i7=find(Mtip<=.7);
		if length(i7)>0
			Cti(i7)=prop2bladeM7(Ji(i7),Cpi(i7)); %   M=.7, 2-bladed prop
		end
	
		i8=find(Mtip<=.8 & Mtip>.7);
		if length(i8)>0
			CtiHigh(i8)=prop2bladeM8(Ji(i8),Cpi(i8)); %   M=.8, 2-bladed prop
			CtiLow(i8) =prop2bladeM7(Ji(i8),Cpi(i8)); %   M=.7, 2-bladed prop
			Cti(i8)=CtiLow(i8)+(Mtip(i8)-.7).*(CtiHigh(i8)-CtiLow(i8))/.1;
		end
	
		i9=find(Mtip<=.9 & Mtip>.8);
		if length(i9)>0
			CtiHigh(i9)=prop2bladeM9(Ji(i9),Cpi(i9)); %   M=.9, 2-bladed prop
			CtiLow(i9) =prop2bladeM8(Ji(i9),Cpi(i9)); %   M=.8, 2-bladed prop
			Cti(i9)=CtiLow(i9)+(Mtip(i9)-.8).*(CtiHigh(i9)-CtiLow(i9))/.1;
		end
	
		i95=find(Mtip<=.95 & Mtip>.9);
		if length(i95)>0
			CtiHigh(i95)=prop2bladeM95(Ji(i95),Cpi(i95)); %   M=.95, 2-bladed prop
			CtiLow(i95) =prop2bladeM9(Ji(i95),Cpi(i95)); %   M=.8, 2-bladed prop
			Cti(i95)=CtiLow(i95)+(Mtip(i95)-.9).*(CtiHigh(i95)-CtiLow(i95))/.05;
		end
	
		i98=find(Mtip<=.98 & Mtip>.95);
		if length(i98)>0
			CtiHigh(i98)=prop2bladeM98(Ji(i98),Cpi(i98)); %   M=.98, 2-bladed prop
			CtiLow(i98) =prop2bladeM95(Ji(i98),Cpi(i98)); %   M=.95, 2-bladed prop
			Cti(i98)=CtiLow(i98)+(Mtip(i98)-.95).*(CtiHigh(i98)-CtiLow(i98))/.03;
		end
	
		i1=find(Mtip>.98);
		if length(i1)>0
			Cti(i1)=prop2bladeM98(Ji(i1),Cpi(i1)); %   M=.98, 2-bladed prop
			disp('Tip mach number is >.98')
		end
	
		
	elseif Nblades==3                  % For 3-bladed  prop
		i7=find(Mtip<=.7);
		if length(i7)>0
			Cti(i7)=prop3bladeM7(Ji(i7),Cpi(i7)); %    M=.7, 3-bladed prop
		end
	
		i8=find(Mtip<=.8 & Mtip>.7);
		if length(i8)>0
			CtiHigh(i8)=prop3bladeM8(Ji(i8),Cpi(i8)); %    M=.8, 3-bladed prop
			CtiLow(i8) =prop3bladeM7(Ji(i8),Cpi(i8)); %    M=.7, 3-bladed prop
			Cti(i8)=CtiLow(i8)+(Mtip(i8)-.7).*(CtiHigh(i8)-CtiLow(i8))/.1;
		end
	
		i9=find(Mtip<=.9 & Mtip>.8);
		if length(i9)>0
			CtiHigh(i9)=prop3bladeM9(Ji(i9),Cpi(i9)); %    M=.9, 3-bladed prop
			CtiLow(i9) =prop3bladeM8(Ji(i9),Cpi(i9)); %    M=.8, 3-bladed prop
			Cti(i9)=CtiLow(i9)+(Mtip(i9)-.8).*(CtiHigh(i9)-CtiLow(i9))/.1;
		end
	
		i95=find(Mtip<=.95 & Mtip>.9);
		if length(i95)>0
			CtiHigh(i95)=prop3bladeM95(Ji(i95),Cpi(i95)); %    M=.95, 3-bladed prop
			CtiLow(i95) =prop3bladeM9(Ji(i95),Cpi(i95));  %    M=.8,  3-bladed prop
			Cti(i95)=CtiLow(i95)+(Mtip(i95)-.9).*(CtiHigh(i95)-CtiLow(i95))/.05;
		end
	
		i98=find(Mtip<=.98 & Mtip>.95);
		if length(i98)>0
			CtiHigh(i98)=prop3bladeM98(Ji(i98),Cpi(i98)); %    M=.98, 3-bladed prop
			CtiLow(i98) =prop3bladeM95(Ji(i98),Cpi(i98)); %    M=.95, 3-bladed prop
			Cti(i98)=CtiLow(i98)+(Mtip(i98)-.95).*(CtiHigh(i98)-CtiLow(i98))/.03;
		end
	
		i1=find(Mtip>.98);
		if length(i1)>0
			Cti(i1)=prop3bladeM98(Ji(i1),Cpi(i1)); %    M=.98, 3-bladed prop
			disp('Tip mach number is >.98')
		end
	else
		disp('ERROR Number of blades is not 2 or 3')
	end
	
	thrust=Cti.*rho.*n.^2*Dft^4;

end



