%% 
% codice cassa d'aria - come in slide 111 moto vario

n    = 5000;       % grandezza vettore tempo
dt   = 0.02;       % passo temporale (nota: n e dt vanno scelti adeguatamente in base a omega)
endTime = n*dt;    % intervallo di tempo n*dt osservato
t = dt:dt:endTime; % vettore tempo

% grandezze variabili 
z    = zeros(n,1); % altezza pelo libero nella cassa d'aria, sol numerica
z_an = zeros(n,1); % altezza pelo libero nella cassa d'aria, sol analitica 
v    = zeros(n,1); % velocitą del fluido nella condotta, sol numerica
p    = zeros(n,1); % pressione nella cassa d'aria, sol numerica
p_an = zeros(n,1); % pressione nella cassa d'aria, sol analitica
p_ta = zeros(n,1); % pressione nella cassa d'aria, sol approx Taylor

%***********************************************************
%********** DATI INPUT *************************************
%***********************************************************
lam  = 0.005;     % coeff lambda - formula Darcy-Weiesbac - provare es. 0.05 o 0.001

L    = 1000;     % lunghezza condotta
g    = 9.8;      % acc gravitą
rho  = 1000;     % densitą fluido 
patm = 101330;   % pressione atmosferica
kk   = 1.0;      % coeff trasf isoterma kk=1, adiabatica kk=1.4
Ac   = 0.5;      % area sezione cassa d'aria
zm   = 1.0;      % quota massima cassa d'aria 
h0   = 40.0;     % carico totale serbatoio, rispetto al sis. rif. posizionato in z(0)
Qp   = 0.05;     % portata indotta dalla pompa
D    = 0.2;      % diametro condotta

% cond iniziali per t=0
v(1) = 0.0;
z(1) = 0.0;
p(1) = 392400;   % p(1) = h0*rho*g in caso sia tutto in quiete

% **********************************************************


% Calcolo di alcune variabili utili

% area condotta
A     = pi*(D/2.0)^2;

% frequenza oscillazioni  
omega = sqrt( (9.8*A/(1000*Ac))*kk*(1 + (p(1)+patm)/(9.8*1000)) ) 

% periodo oscillazioni
T     = 2.0*pi/omega 

% parametri c1 e c2 
% [ calcolati tramite condizione iniziale al tempo t=0 
%   e l'eq. di continuitą Qp = A*V + Ac*z' ]
c1 = (Qp - A*v(1))/(Ac*omega); c2 = z(1);

% main loop : a partire dal passo k=2, calcolo J = lambda*|V|V / 2gD, poi v,z e p 
for k = 2 : 1 : n;
  % soluzione numerica
  J       = lam*v(k-1)*abs(v(k-1))/(2.0*g*D); 
  v(k)    = v(k-1) + dt*(g/L)*( z(k-1) + p(k-1)/(rho*g) - h0 - L*J);
  z(k)    = z(k-1) + dt*(Qp - A*v(k-1))/Ac;
  p(k)    = ((p(1) + patm)*(zm - z(1))^kk)/(zm - z(k))^kk - patm;
  % soluzione analitica
  z_an(k) =  c1*sin(omega*k*dt) + c2*cos(omega*k*dt);
  p_an(k) = ( (p(1) + patm)*(zm - z(1))^kk )/(zm - z_an(k))^kk - patm;  
  p_ta(k) =    p(1) + (kk/zm)*( p(1) + patm )*z_an(k) ;  
end
%
%
% plot 
figure(1) % figura: andamento di z nel tempo
    plot(t(1:n),z_an(1:n),"color","black") % sol analitica J=0
    hold on
    grid on
    plot(t(1:n),z(1:n),"color","red") % sol numerica, J diverso da 0
    legend('sol analitica','sol numerica')
    xlabel('t (s)')
    ylabel('z (m)')

figure(2) % figura: andamento di p nel tempo
    plot(t(1:n),p_an(1:n),"color","black") % sol analitica J=0
    hold on
    grid on
    plot(t(1:n),p(1:n),"color","red")  % sol numerica, J diverso da 0
    plot(t(1:n),p_ta(1:n),"color","blue")  % sol approx
    legend('sol analitica','sol numerica','sol approx')
    xlabel('t (s)')
    ylabel('p (Pa)')