close all
clear all
rho_w =1025;
rho = 1800;
nu = 1e-6;   %m2/s  Kinematic viscosity of water
sigma = 72e-3  %  N/m   Surface tension of water
g = 9.81;


Cv = .9;
Cc = .61;

H = 50:50:200;  % Meters
Do = [3 2 1]*.0254;      % Turbine Diam (m)
%Do = [5]*.0254;      % Turbine Diam (m)
d = [.05:.001:.2]*.0254;  % Nozzle Diam (m)


K2 = 0.242738;   %For round orifice
for D=Do  %Individual plot for each turbine diameter

for i = 1:length(H)
  for j = 1:length(d)
    Hprime = rho_w*H(i)/rho;
    Re(i,j) = sqrt(2*g*Hprime)*d(j)/nu;
    Cv(i,j) = 0.6; %sqrt(1-K2*(Cc^1.5)/sqrt(Re(i,j)));
    Pressure = rho_w*g*H(i);
    Vj(i,j) = Cv(i,j)*sqrt(2*Pressure/rho);  % Ideal jet velocity
    %Vj(i,j) = Cv(i,j)*sqrt(2*g*Hprime);  % Ideal jet velocity
    w(i,j) = (Vj(i,j)/2)/(D/2);
    Q(i,j) = Cc*Vj(i,j)*0.25*pi*d(j)^2;
    Ns(i,j) = SpecificSpeed(w(i,j),g,Q(i,j),Hprime);    % SpecificSpeed(w(i,j),Q(i,j),g,H(i));
    P(i,j) = Q(i,j)*Pressure;
    T(i,j) = P(i,j)/w(i,j);
  end
end


NumPlots = 5; k = 1;
figure('Position',[1000 228 815 1270]);

ah(k) = subplot(NumPlots,1,k); k = k+1;
hold all
for i = 1:length(H)
  plot(d/.0254,Ns(i,:));
  leg_str{i} = ['H = ' num2str(H(i)) ' m'];
end
xlabel('Jet Diameter (in)')
ylabel('Specific Speed');
legend(leg_str,'Location','NorthWest');
title(['Turbine Diam = ' num2str(D/.0254) ' in'])
grid on;

ah(k) = subplot(NumPlots,1,k); k = k+1;
hold all
for i = 1:length(H)
  plot(d/.0254,w(i,:)*60/(2*pi));
  leg_str{i} = ['H = ' num2str(H(i)) ' m'];
end
xlabel('Jet Diameter (in)')
ylabel('RPM');
legend(leg_str,'Location','NorthWest');
title(['Turbine Diam = ' num2str(D/.0254) ' in'])
grid on;

ah(k) = subplot(NumPlots,1,k); k = k+1;
hold all
for i = 1:length(H)
  plot(d/.0254,T(i,:));
  leg_str{i} = ['H = ' num2str(H(i)) ' m'];
end
xlabel('Jet Diameter (in)')
ylabel('Torque (N-m)');
legend(leg_str,'Location','NorthWest');
title(['Turbine Diam = ' num2str(D/.0254) ' in'])
grid on;


ah(k) = subplot(NumPlots,1,k); k = k+1;
hold all;
for i = 1:length(H)
  plot(d/.0254,P(i,:));
  leg_str{i} = ['H = ' num2str(H(i)) ' m'];
end
xlabel('Jet Diameter (in)')
ylabel('Power (W)');
legend(leg_str,'Location','NorthWest');
title(['Turbine Diam = ' num2str(D/.0254) ' in'])
grid on;


ah(k) = subplot(NumPlots,1,k); k = k+1;
hold all;
for i = 1:length(H)
  plot(d/.0254,(1/1000)./Q(i,:));
  leg_str{i} = ['H = ' num2str(H(i)) ' m'];
end
xlabel('Jet Diameter (in)')
ylabel('Seconds per Liter');
legend(leg_str,'Location','NorthWest');
title(['Turbine Diam = ' num2str(D/.0254) ' in'])
grid on;

end