clear all
close all

% Definizione delle costanti del modello
A=1;                        % Sezione del serbatoio (m^2)
a=.25;                      % Sezione uscita (m^2)
g=9.8;                      % Accelerazione di gravita' (m/s^2)

% Definizione dei parametri temporali
Tend=30;                    % Tempo di simulazione (s) 
Ts=0.01;                    % Step di simulazione (s)
N=floor(Tend/Ts);           % Numero di step
T=linspace(0,N*Ts,N+1);     % Vettore tempo (s)

% Definizione dei parametri dell'ingresso
bias=2;                     % Bias (m^3/s)
amp=1;                      % Ampiezza (m^3/s)
omega=0.33;                 % Frequenza (rad/s)
U=bias+amp*sin(omega*T');   % Vettore flusso di ingresso (m^3/s)

% Simulazione del sistema
h0=1;                       % Altezza iniziale (m)
H=h0;                       % Inizializzazione del vettore altezza (m) 
for k=1:N
    H=[ H ; H(k)+Ts/A*(-a*sqrt(2*g*H(k))+U(k)) ];
end
Q=a*sqrt(2*g*H);            % Vettore flusso di uscita (m^3/s)

% Visualizzazione
subplot(121)
plot(T,H)
grid
xlabel('Tempo (s)')
ylabel('Altezza h(t) (m)')

subplot(122)
plot(T,[ Q U ])
grid
xlabel('Tempo (s)')
ylabel('Flussi di ingresso e uscita (m^3/s)')