% sistema di condotti che collega tre serbatoi
% Determinare la portata di acqua Q in ciascun condotto
% si trascurano le perdite accidentali e si assume p1=p2=p3=0
% Q1+Q2+Q3=0
% hJ=zJ+pJ/(rho*g) HGL Hydraulic Grade Line o quota piezometrica del giunto
% z1=zJ+pJ/(rho*g)+delta_h1 teo bernoulli modificato tra 1 e J 
% z2=zJ+pJ/(rho*g)+delta_h2 teo bernoulli modificato tra 2 e J 
% z3=zJ+pJ/(rho*g)+delta_h3 teo bernoulli modificato tra 3 e J 
% delta_h1=f1 * L1/d1 * u1^2/(2*g) = z1-hJ
% delta_h2=f2 * L2/d2 * u2^2/(2*g) = z2-hJ
% delta_h3=f3 * L3/d3 * u3^2/(2*g) = z3-hJ
clear all;
close all;
clc;

g=9.81; % m/s^2 accelerazione di gravità

rho = 1000; % kg/m^3 densità acqua
nu = 1.02e-6; % m^2/s viscosità cinematica acqua

% pipe 1
L1=100; % m lunghezza
d1=8e-2; % m diametro
e1=0.24e-3; % m rugosità
ed1=e1/d1; % rugosità relativa
z1=20; % m quota del serbatoio

% pipe 2
L2=150; % m lunghezza
d2=6e-2; % m diametro
e2=0.12e-3; % m rugosità
ed2=e2/d2; % rugosità relativa
z2=100; % m quota del serbatoio

% pipe 3
L3=80; % m lunghezza
d3=4e-2; % m diametro
e3=0.20e-3; % m rugosità
ed3=e3/d3; % rugosità relativa
z3=40; % m quota del serbatoio

eps_Q=1e-2;

% quota piezometrica del giunto hJ=zJ+pJ/(rho*g) 
hJ=0:10:z2;
for i=1:length(hJ)
    
    % calcolo Q1
    delta_h1=z1-hJ(i);
    Q1 = calcolo_Q(ed1,L1,d1,delta_h1,nu,g,eps_Q);
    
    % calcolo Q2
    delta_h2=z2-hJ(i);
    Q2 = calcolo_Q(ed2,L2,d2,delta_h2,nu,g,eps_Q);
    
    % calcolo Q3
    delta_h3=z3-hJ(i);
    Q3 = calcolo_Q(ed3,L3,d3,delta_h3,nu,g,eps_Q);
    
    Q(i)=Q1+Q2+Q3; % deve essere zero
    
    % Q<0 hJ troppo alto
    % Q>0 hJ troppo basso
    
end

figure;
hold on;
plot(hJ,Q);
grid on;

% metodo di bisezione per determinare lo zero (Q=0)

% inizializzazione
hJa=20;
hJb=40;

Q1 = calcolo_Q(ed1,L1,d1,z1-hJa,nu,g,eps_Q);
Q2 = calcolo_Q(ed2,L2,d2,z2-hJa,nu,g,eps_Q);
Q3 = calcolo_Q(ed3,L3,d3,z3-hJa,nu,g,eps_Q);
Qa=Q1+Q2+Q3; % deve essere zero

Q1 = calcolo_Q(ed1,L1,d1,z1-hJb,nu,g,eps_Q);
Q2 = calcolo_Q(ed2,L2,d2,z2-hJb,nu,g,eps_Q);
Q3 = calcolo_Q(ed3,L3,d3,z3-hJb,nu,g,eps_Q);
Qb=Q1+Q2+Q3; % deve essere zero

hJc=(hJa+hJb)/2;
Q1 = calcolo_Q(ed1,L1,d1,z1-hJc,nu,g,eps_Q);
Q2 = calcolo_Q(ed2,L2,d2,z2-hJc,nu,g,eps_Q);
Q3 = calcolo_Q(ed3,L3,d3,z3-hJc,nu,g,eps_Q);
Qc=Q1+Q2+Q3; % deve essere zero

s=sprintf('hJ=%0.6f m   Q1=%0.6f m^3/s   Q2=%0.6f m^3/s   Q3=%0.6f m^3/s   Q=%0.6f m^3/s',hJc,Q1,Q2,Q3,Qc);
disp(s);

% ciclo while
eps=1e-4;
while abs(Qc)>eps
    if Qc>0
        hJa=hJc; Qa=Qc;
    else
        hJb=hJc; Qb=Qc;
    end
    hJc=(hJa+hJb)/2;
    Q1 = calcolo_Q(ed1,L1,d1,z1-hJc,nu,g,eps_Q);
    Q2 = calcolo_Q(ed2,L2,d2,z2-hJc,nu,g,eps_Q);
    Q3 = calcolo_Q(ed3,L3,d3,z3-hJc,nu,g,eps_Q);
    Qc=Q1+Q2+Q3; % deve essere zero
    s=sprintf('hJ=%0.6f m   Q1=%0.6f m^3/s   Q2=%0.6f m^3/s   Q3=%0.6f m^3/s   Q=%0.6f m^3/s',hJc,Q1,Q2,Q3,Qc);
    disp(s);
end

% il flusso va verso il serbatoio 1 bilanciato dal flusso dal serbatoio 2 e
% dal serbatoio 3