echo off
format compact
clear
disp(' ')
disp('Start new run')
disp('Propeller Analysis using Goldsteins Classical Vortex Theory')
disp('  This code works for two bladed propellers only.')
disp(' ')

% Input constants
Din=12;      % propeller diameter (inches)
Pin=12;      % pitch of the propeller (inches)
RPM=7720;    % RPM of propeller
Vmph=135*60/88 ;
%Vmph=0	% airspeed in mph
rho=.002378; % air density (slug/ft**2) (sea level)
%rho=.002065; % air density (slug/ft**2) at 4750 feet

% Pin gives geometric angle to the flat part of the 
%   rear of the propeller
aoldeg=-6;   % angle of zero lift of the propeller (degrees)
             %  measured from mean chord line (typically negative)
beta0deg=.5; % angle from flat part of the prop to mean chord line
a=2*pi;      % lift curve slope of propeller
Cd0=.00655;	 % 2-d minimum drag coefficient
k=.01;       % Cd = Cd0+k*Cl*Cl
B=2 ;        % number of blades (2 for standard type propeller)


% input nondimensional properties at each radial location 
% cR=c/R, x=r/R
x=[.3,.35,.4,.45,.5,.55,.6,.65,.7,.75,.8,.85,.9,.95,1.];
cR=.09*ones(size(x));
% END OF INPUTS

% Display inputs section
disp(['Diameter Din= ',num2str(Din),' inches'])
disp(['Pitch Pin= ',num2str(Pin),' inches'])
disp(['Angle of zero lift aoldeg= ',num2str(aoldeg),' degrees'])
disp(['Angle Flat side to mean chord beta0deg= ',num2str(beta0deg),' degrees'])
disp(['RPM= ',num2str(RPM),' rpm'])
disp(['2-d Lift curve slope a= ',num2str(a),' per rad'])
disp(['2-d Min drag coef Cd0= ',num2str(Cd0)])
disp(['k from 2-d drag polar= ',num2str(k)])
disp(['Number of blades B= ',num2str(B)])
disp(['air density rho= ',num2str(rho),' slug/ft^3'])
disp(['Airspeed Vmph= ',num2str(Vmph),' mph'])



% derived constants
V=Vmph*88/60; % airspeed in ft/sec
D=Din/12;     % Diameter in feet
R=D/2;        % Radius in feet
n=RPM/60;     % propeller frequency (rev/sec) or (hz)
omega=2*pi*n; % frequency of revolution of the propeller (rad/sec)
lamda=V/(omega*R);
r2d=180/pi;
Vt=omega*R;   % tip velocity (ft/sec)
J=V/(n*D);    % advance ratio

% Output scalar constants
disp(['Airspeed V= ',num2str(V),' ft/sec'])
disp(['Propeller Diameter D= ',num2str(D),' feet'])
disp(['Propeller Radius R= ',num2str(R),' feet'])
disp(['propeller RPS n= ',num2str(n),' hertz'])
disp(['omega= ',num2str(omega),' rad/sec'])
disp(['lamda= ',num2str(lamda)])
disp(['r2d= ',num2str(r2d),' deg/rad'])
disp(['Tip speed Vt= ',num2str(Vt),' ft/sec'])
disp(['Advance Ratio J= ',num2str(J)])
disp(' ')

%derived section constants
c=R*cR;    % chord in feet
cin=c*12;  % chord in inches
%disp(' ')
%disp('x in r/R and is nondimensional, cR=c/R, cin in the chord in inches')
%echo on
%disp([x',      cR',       cin'])
%echo off
%disp(' ')

beta1=atan(((Pin/Din)/pi)./x);
beta=beta1+(beta0deg-aoldeg)/r2d;
sigma=B*c/(pi*R);
r=x*R;
Vr=Vt*sqrt(x.*x+lamda*lamda);
phi=atan(lamda./x);
WtVt=.02*ones(size(c));  %initial guess
nr=length(c);
nr1=nr-1;
aiold=zeros(size(c));
for ii=1:40
	WaVt(1:nr1)=.5*(-lamda+sqrt(lamda*lamda+4*WtVt(1:nr1).*(x(1:nr1)-WtVt(1:nr1))));
	ai(1:nr1)=atan(WtVt(1:nr1)./WaVt(1:nr1))-phi(1:nr1);
	%ai(1:nr1)=atan((V+WaVt(1:nr1)*Vt)./(omega*r(1:nr1)-WtVt(1:nr1)*Vt))-phi(1:nr1);
	e=sum(abs(ai(1:nr1)-aiold(1:nr1)));
	iter=['Loop index= ',num2str(ii),'  error= ',num2str(e)];
	%disp(iter)
	if e<.0001 ; break; end
	aiold(1:nr1)=ai(1:nr1);
	Cl(1:nr1)=a*(beta(1:nr1)-ai(1:nr1)-phi(1:nr1));
	VeVt(1:nr1)=sqrt((lamda+WaVt(1:nr1)).^2+(x(1:nr1)-WtVt(1:nr1)).^2);
	gamma(1:nr1)=.5*c(1:nr1).*Cl(1:nr1).*VeVt(1:nr1)*Vt;
	sinphialp(1:nr1)=sin(phi(1:nr1)+ai(1:nr1));
	kappa(1:nr1)=kappa2(x(1:nr1),sinphialp(1:nr1));
	WtVt(1:nr1)=B*gamma(1:nr1)./(4*pi*Vt*r(1:nr1).*kappa(1:nr1));
end
Cl(nr)=0;
ai(nr)=beta(nr)-phi(nr);
VrVt=sqrt(lamda*lamda+1);
WaVt(nr)=VrVt*sin(ai(nr))*cos(ai(nr)+phi(nr));
WtVt(nr)=VrVt*sin(ai(nr))*sin(ai(nr)+phi(nr));
VeVt(nr)=sqrt((lamda+WaVt(nr))^2+(x(nr)-WtVt(nr))^2);
kappa(nr)=0;

Cd=Cd0+k*Cl.*Cl;
ZT=(pi/8)*(J*J+pi*pi*(x.*x)).*sigma;
ZP=pi*ZT.*x;
dCTdx=ZT.*(Cl.*cos(phi+ai)-Cd.*sin(phi+ai));
dCPdx=ZP.*(Cl.*sin(phi+ai)+Cd.*cos(phi+ai));

% Overall propeller performance
CT=trapi(dCTdx,x);
CP=trapi(dCPdx,x);
eta=CT*J/CP;
T=CT*rho*n^2*D^4;
P=CP*rho*n^3*D^5;
HP=P/550;
Pwatt=1.356*P;
torque=P/omega;
Clmax=max(Cl);
Toz=T*16;
PestWatts=1.31*D^4*(Pin/12)*(RPM/1000)^3;
% The above approximate formula works for 
% Top Flite, Zinger and Master Airscrews reasonably well.
% For Rev Up props subract .5 in from the pitch.
% For APC props use constant 1.11 instead of 1.31.
% For thin carbon fiber folding props use 1.18 instead of 1.31.
% Ref: Electric Motor Handbook, by Robert J. Boucher,
%    AstroFlight, Inc.

%echo on
%dat=[ x',     beta',     phi',     kappa',     ai',     WtVt']
%disp(' ')
%dat2=[x',     WaVt',     VeVt',    Cl',     dCTdx',    dCPdx']
%echo off

disp(' ')
disp(['Advance Ratio J= ',num2str(J)])
disp(['Thrust Coefficient CT= ',num2str(CT)])
disp(['Power Coefficient CP= ',num2str(CP)])
disp(['Propeller efficiency eta= ',num2str(eta)])
disp(['Speed V= ',num2str(V),' ft/sec'])
disp(['RPM= ',num2str(RPM),' rpm'])
disp(['Thrust T= ',num2str(T),' pounds'])
disp(['Thrust Toz= ',num2str(Toz),' ounces'])
disp(['Power used P= ',num2str(P),' ft*lbf/sec'])
disp(['Horsepower used HP= ',num2str(HP),' HP'])
disp(['Power used Pwatt= ',num2str(P),' watts'])
disp(['Torque used Q= ',num2str(torque),' ft*lbf'])
disp(['Torque used Q= ',num2str(torque*192),' in-oz'])
disp(['Clmax= ',num2str(Clmax)])
disp(' ')
disp(['Estimated power used, PestWatts= ',num2str(PestWatts),' watts, Ref: Boucher'])
disp(' ')

subplot(211)
plot(x,dCTdx)
z=axis;
axis([0,1,0,z(4)])
xlabel('nondimensional radial location')
ylabel('dCTdx')
subplot(212)
plot(x,dCPdx)
z=axis;
axis([0,1,0,z(4)])
xlabel('nondimensional radial location')
ylabel('dCPdx')

disp(' For model aircraft propellers this code underestimates') 
disp(' the power required. The underestimation is worse for RPM>10,000')
disp(' where it may underestimate by a factor of .5')
disp(' For RPM <10,000 the factor is about .75')
