clear all;
close all;

rho = 1025;  %kg/m^3
g = 9.81;
PaPerPsi = 6894.75729; 
MPerIn = .0254;

Po = 300;  %psi
Pi = 25;  %psi
Po = Po*PaPerPsi;
Pi = Pi*PaPerPsi;

f = 0.035;

lj = .25; %inches
dj = .1:.05:.5;  %inches
lj = lj*MPerIn;
dj = dj*MPerIn;

Aj = pi*dj.^2/4;

do = 1.5;  %inches
lc = 1.0;  %inches
do = do*MPerIn;
lc = lc*MPerIn;


th = atan2(do-dj/2,lc);

k = 0.5*(1-(dj/do).^2).*(sin(th).^.5)./(dj/do).^4;  %Applies for th > 45deg

Vj = ((2*(Po-Pi)/rho)./(1+f*lj./dj+k)).^.5;


Q = Aj.*Vj;

hl = (Vj.^2/2).*(f*lj./dj+k);



[AX,H1,H2] = plotyy(dj/MPerIn,Vj,dj/MPerIn,Q*1000)
xlabel('jet diameter (in)');
set(get(AX(1),'Ylabel'),'String','Jet Velocity (m/s)') 
set(get(AX(2),'Ylabel'),'String','Flow Rage (L/s)') 




figure
plot(dj/MPerIn,hl/g)







dp = 2:1:6;  %Diameter of pelton wheel (in).
dp = dp*MPerIn;

Vp = Vj/2;  %

hold on;
for i = 1:length(dp);
wp{i} = 2*Vp/dp(i);
end
figure
plot(dj/MPerIn,wp{1}*60/(2*pi),dj/MPerIn,wp{2}*60/(2*pi),dj/MPerIn,wp{3}*60/(2*pi),dj/MPerIn,wp{4}*60/(2*pi),dj/MPerIn,wp{5}*60/(2*pi));  
xlabel('jet diameter (in)');
ylabel('rpm');
title('Wheel speed versus jet diamter for a range of wheel diameters (2in->6in)');
legend('Wheel Diam = 2in','Wheel Diam = 3in','Wheel Diam = 4in','Wheel Diam = 5in','Wheel Diam = 6in','Location','NorthWest')