clear all
close all

g = 9.81;
rho = 1025;  %kg/m^3
nu = 0.000001;  %kinematic viscosity

h = 200;  %m

L = .25;  %Nozzle length (cm)
L = L/100;  %m


d = 1:.25:5;  %mm
d = d/1000;  %m

D = 6*.0254;  %Pelton wheel diameter

f = 0.04;  %Pipe flow friction coefficient  %Need to make a function of reynolds number, makes Vi computation implicit.
K = 0.04;  %0.4;  r/d of entrance > 0.15;

Pi = 0;  
Po = rho*g*h;

dP = Po-Pi;


Vi = sqrt(2*dP./(rho*(1+f*L./d+K)));

Q = pi*d.^2.*Vi/4;
Power = rho*Q.*Vi.^2/2;

Rey = Vi.*d/nu;

eta = Power./(rho*g*h*Q);  % !!!!!!!!!! NEED TO CHECK !!!!!!!!!!!

Vavg = 4*Q./(pi*d.^2);

w = Vavg/(2*pi*D);


figure('PaperType','usletter','PaperPositionMode','manual','PaperPosition',[.75,.75,8.5-.75,11-.75])

h1 = subplot('Position',[.1 .54 .8 .4]);
h2 = subplot('Position',[.1 .04 .8 .4]);

subplot(h1)
plot(d*100,100*100*100*Q/1000)
xlabel('Jet Diam (cm)');
ylabel('Flow Rate  (L/s)');

subplot(h2)
plot(d*100,Power)
xlabel('Jet Diam (cm)');
ylabel('Power  (W)');


figure
plot(d*100,1000./(100*100*100*Q))
xlabel('Jet Diam (cm)');
ylabel('Flow Rate  (seconds per liter)');



figure
plot(d*100,w*60/(2*pi));
xlabel('Jet Diam (cm)');
ylabel('Wheel Angular Velocity (rpm)');




figure
plot(d*100,eta);
xlabel('Jet Diam (cm)');
ylabel('Efficiency');
