%close all
%clear all

h = 200;
rho = 1025;
g = 9.81;
SG = rho/1000;

V_ideal = sqrt(2*g*h);

dP = rho*g*h;  %N/m^2
dP_psi = dP/6895;  %psi

d_in  = [0.0008 0.0012 0.0016 0.0020 0.0024 0.0028 0.0031 0.0035 0.0040 0.0050 0.0060 0.0070 0.0080 0.0090 0.0100 0.0110 0.0120 0.0130 0.0140 0.0150 0.0160 0.0170 0.0180 0.0190 0.0200 0.0210 0.0220 0.0230 0.0240 0.0250 0.030 0.035 0.040 0.043 0.047 0.052 0.060 0.063 0.067 0.070 0.073 0.079 0.081 0.086 0.089 0.094 0.096 0.100 0.104 0.109 0.120 0.125];
Cv = [0.00001 0.00003 0.000053 0.000090 0.00012 0.00017 0.00022 0.00028 0.00035 0.00061 0.00086 0.0012 0.0015 0.0019 0.0025 0.0028 0.0034 0.0038 0.0043 0.0050 0.0055 0.0067 0.0073 0.0080 0.0088 0.0096 0.011 0.012 0.013 0.014 0.020 0.028 0.036 0.041 0.048 0.059 0.081 0.088 0.100 0.110 0.120 0.140 0.150 0.170 0.180 0.200 0.210 0.23 0.25 0.27 0.34 0.37];

d = d_in*0.0254;

Q_gpm = Cv*sqrt(dP_psi/SG);
Q = Q_gpm*3.785/60;  %L/s

V = (Q/1000)./(0.25*pi*d.^2);  %m/s

eff = 0.5*V.^2/(g*h);

Cd = (Q/1000)./(0.25*pi*d.^2*sqrt(2*dP/rho));


figure
plot(d_in,V,d_in,V_ideal*ones(1,length(d_in)));
title('From Bird Precision Cv Values');

figure
plot(d_in,Q_gpm)
title('From Bird Precision Cv Values');
xlabel('diam (in)');
ylabel('Q (gpm)');


figure
plot(d_in,eff^2,d_in,Cd);
title('From Bird Precision Cv Values');
legend('eff','C_discharge');
xlabel('diam (in)');