close all
clear all


global g h d l nu

d_meas = [.02 .04 .06 .08 .1];
Q_meas300 = [.3/40 0 1.3/20 0 1.9/10];
Q_meas200 = [.25/40 0 1.05/20 0 2.08/15];
Q_meas100 = [.285/60 0 .30/20 0 1.6/15];

%Q_meas300_45deg = [0 0 1.55/20 0 0];
%Q_meas200_45deg = [0 0 1.285/20 0 0];
%Q_meas100_45deg = [0 0 .9/20 0 0];

Q_meas400_radius = [0 0 1.08/10 0 0];
Q_meas300_radius = [0 0 0.875/10 1.7/10 0];
Q_meas200_radius = [0 0 0.7/10 1.375/10 1.725/10];
Q_meas100_radius = [0 0 .525/10 .98/10 1.6/10];


w_meas300_radius = [0 0 1500  ];


g = 9.81;  %m/s^2
%rho = 1680;  
rho = 1025;  %kg/m^3
rho_air = 1.25;  %kg/m^3
nu = 1.3E-6; %m^2/s
%nu = .5E-6;
nu_air = 1.4E-5;

h = 300/1.5; %40/1.5;  %m
l = 0.200*.0254;  %m
V_ideal = sqrt(2*g*h);


%Nozzle Behavior

diam_in = 0.001:.001:0.125;  %in
diam = diam_in*0.0254;  %m
for i = 1:length(diam)
d = diam(i);
Q_ideal(i) = V_ideal*pi*d^2/4;
V(i) = fzero(@mod_bernoulli,V_ideal);
Q(i) = V(i)*pi*d^2/4;
Re(i) = V(i)*d/nu;

eff(i) = (V(i)/V_ideal)^2;
end

BirdOrifices

paperwidth = 8.5;
paperheight = 11.0;
h_Orifices = figure('PaperPositionMode','manual','PaperSize',[paperwidth paperheight],'PaperPosition',[0 0 paperwidth paperheight],'Units','inches','Position',[2 -1 paperwidth paperheight],'Visible','on'); 

left = 1.0; bottom = 8.4; width = 6.5; height = 2.0;
h_VelVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 5.9; width = 6.5; height = 2.0;
h_EffVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 3.4; width = 6.5; height = 2.0;
h_QVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 0.9; width = 6.5; height = 2.0;
h_CdVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);


subplot(h_VelVsdiam);
plot(diam_in,V,d_in_Bird,V_Bird);
legend('Computed','From Bird Precision');
%xlabel('diam (in)');
ylabel('Vel (m/s)');

subplot(h_EffVsdiam);
plot(diam_in,eff,d_in_Bird,eff_Bird);
legend('Computed','From Bird Precision');
%xlabel('diam (in)');
ylabel('eff');

subplot(h_QVsdiam);
plot(diam_in,Q*1000,d_in_Bird,Q_Bird,d_meas,Q_meas300,'x',d_meas,Q_meas300_radius,'o')
legend('Computed','From Bird Precision');
%xlabel('diam (in)');
ylabel('Q (lps)');

A = pi*diam.^2/4;
subplot(h_CdVsdiam);
plot(diam_in,Q./(A*sqrt(2*g*h)),d_in_Bird,Cd_Bird)
legend('Computed','From Bird Precision');
xlabel('diam (in)');
ylabel('Cd');

print -depsc2 Nozzle.ps




%Pelton Wheel Behavior

d_pelton_in = 8.0;  %Pelton Wheel Diameter (in)
d_pelton = d_pelton_in*0.0254;  %Pelton Wheel Diameter (m)


paperwidth = 8.5;
paperheight = 11.0;
h_summary = figure('PaperPositionMode','manual','PaperSize',[paperwidth paperheight],'PaperPosition',[0 0 paperwidth paperheight],'Units','inches','Position',[2 -1 paperwidth paperheight],'Visible','on'); 

left = 1.0; bottom = 8.4; width = 6.5; height = 2.0;
h_wVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 5.9; width = 6.5; height = 2.0;
h_TVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 3.4; width = 6.5; height = 2.0;
h_PVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 0.9; width = 6.5; height = 2.0;
h_effVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);

for dp = d_pelton  
w = V/dp;  %Optimum Speed of runner
F = rho*Q.*V;
T = F.*dp/2;
eff_p = 0.8;    
subplot(h_wVsdiam);
plot(diam_in,w*60/(2*pi));
hold all;

subplot(h_TVsdiam);
plot(diam_in,eff_p*T);
hold all;

subplot(h_PVsdiam);
plot(diam_in,eff_p*T.*w);
hold all;

subplot(h_effVsdiam);
plot(diam_in,eff_p*T.*w./(rho*g*h.*Q));
hold all;

end

subplot(h_wVsdiam)
title('Torque and Speed at Optimum Speed relative to Jet - Perfectly Efficient Pelton Wheel');
%xlabel('diam (in)');
ylabel('w (rpm)');

subplot(h_TVsdiam)
%xlabel('diam (in)');
ylabel('T (N-m)');

subplot(h_PVsdiam)
%xlabel('diam (in)');
ylabel('P (W)');

subplot(h_effVsdiam)
xlabel('diam (in)');
ylabel('Eff');


print -depsc2 PerfectPelton.ps

%Pelton Wheel Losses
d_c_in = 1.0;  %inches
d_c = d_c_in*0.0254;  %meters
N = floor(pi*d_pelton/(2.5*d_c))

Cd_cup = 0.42;  %Assumes no interference.
Cd_cup = 0.8*Cd_cup;  %Assumes less than optimum interfence  (x/D = 1.5 or so).
Cd_cup = 1.1;  %From Measurements

figure('PaperPositionMode','manual','PaperSize',[paperwidth paperheight],'PaperPosition',[0 0 paperwidth paperheight],'Units','inches','Position',[2 -1 paperwidth paperheight],'Visible','on'); 
left = 1.0; bottom = 8.4; width = 6.5; height = 2.0;
h_wVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 5.9; width = 6.5; height = 2.0;
h_TVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 3.4; width = 6.5; height = 2.0;
h_PVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);
left = 1.0; bottom = 0.9; width = 6.5; height = 2.0;
h_effVsdiam  = subplot('Position',[left/paperwidth bottom/paperheight width/paperwidth height/paperheight]);

i = 0;
for dp = d_pelton 
i = i+1;
w = V/dp;
V_cup = V/2;  %Optimum efficiency speed.
Re_cup = V_cup*d_c/nu_air;
Drag_cup = Cd_cup*.5*rho_air*V_cup.^2*pi*d_c^2/4;  
T_cup = N(i)*Drag_cup*dp/2;

F = rho*Q.*V;
T = F.*dp/2;
eff_p = 0.80;    
figure(3)
subplot(h_wVsdiam);
plot(diam_in,w*60/(2*pi));
hold all;

subplot(h_TVsdiam);
plot(diam_in,(T-T_cup));
hold all;

%subplot(h_PVsdiam);
%plot(diam_in,100*(T_cup./T));
%hold all;

subplot(h_PVsdiam);
plot(diam_in,eff_p*(T-T_cup).*w);
hold all;


subplot(h_effVsdiam);
plot(diam_in,eff_p*(T-T_cup).*w./(rho*g*h.*Q));
hold all;

I = .0017;  %kg-m^2
figure(4)
plot(diam_in,2*I./(rho*Q*dp^2));
hold all;

figure(5)  %Specific Speed
plot(diam_in,w.*sqrt(Q)./(g*h)^.75,d_meas,w_meas300_radius.*sqrt(Q_meas300_radius)./(g*h)^.75);
hold all;

figure(6)  %Head Coefficient
plot(diam_in,g*h./((dp*w).^2));
hold all;

figure(7)  %Power Coefficient
plot(diam_in,Q./(w*dp^3));
hold all;

end
figure(3);
subplot(h_wVsdiam)
title('Torque and Speed at Optimum Speed relative to Jet - Including Windage Losses');
%xlabel('diam (in)');
ylabel('w (rpm)');

subplot(h_TVsdiam)
%xlabel('diam (in)');
ylabel('T (N-m)');

subplot(h_PVsdiam)
%foo = axis;
%foo(3) = 0;
%foo(4) = 100;
%axis(foo)
%xlabel('diam (in)');
%ylabel('Torque loss due to windage  (%)');
ylabel('P(W)');

subplot(h_effVsdiam)
foo = axis;
foo(3) = 0;
foo(4) = 1;
axis(foo)
xlabel('diam (in)');
ylabel('Eff');

print -depsc2 Pelton.ps


figure(4)
foo = axis;
foo(4) = 2;
axis(foo);
xlabel('diam(in)');
ylabel('Spin-up Time constant (s)');

figure(5)
xlabel('diam(in)');
ylabel('Specific Speed');

figure(6)
xlabel('diam(in)');
ylabel('Loading Coefficient');

figure(7)
xlabel('diam(in)');
ylabel('Flow Coefficient');