function f = mod_bernoulli(Vj)
global g h d nu;

lj = 0.100*0.0254;
l1 = 0.7*0.0254; 
l2 = 4.7*0.0254;
l3 = 1.15*0.0254;

dj = d;
%d1 = .175*0.0254;
%d2 = 0.290*0.0254;
%d3 = .175*0.0254;
d1 = .290*0.0254;
d2 = 0.290*0.0254;
d3 = .290*0.0254;

Vj;
V1 = Vj*dj^2/d1^2;
V2 = Vj*dj^2/d2^2;
V3 = Vj*dj^2/d3^2;

Rej = Vj*d/nu;
Re1 = V1*d1/nu;
Re2 = V2*d2/nu;
Re3 = V3*d2/nu;


hl_fj = 0.5*(lj/dj)*frictionfactor(Rej).*Vj.^2;
hl_f1 = 0.5*(l1/d1)*frictionfactor(Re1).*V1.^2;
hl_f2 = 0.5*(l2/d2)*frictionfactor(Re2).*V2.^2;
hl_f3 = 0.5*(l3/d3)*frictionfactor(Re3).*V3.^2;

K_12 = (1-(d1/d2)^2)/((d1/d2)^4);
K_23 = 0.5*(1-(d3/d2)^2)/((d3/d2)^4);
%K_3j = 0.05;
K_3j = 0.5*(1-(dj/d3)^2)/((dj/d3)^4)*sqrt(sin(pi/4));
K_sharp = K_3j*V3^2/Vj^2;

hl_12 = (K_12*V2^2)/2;
hl_23 = (K_23*V2^2)/2;
%hl_3j = (K_3j*Vj.^2)/2;
hl_3j = (K_3j*V3.^2)/2;

hl_misc = 0%-155;

hl_total = hl_fj+hl_f1+hl_f2+hl_f3 + hl_12 + hl_23 + hl_3j + hl_misc;


f=.5*Vj.^2-g*h+hl_total;