
% Sean Costello et Sandra Brella, 2014
%--------------------------------------------------------------------------
% SD demo 4: Un exemple de simulation du reservoir de l'exercice 1, serie 4
% QUESTION: Simuler la vidange du reservoir quand on coupe l'alimentation
%--------------------------------------------------------------------------

% Nettoyage de l'espace de travail
clear all; close all; clc;

%--------------------------------------------------------------------------
% Definition des parametres du reservoir 
%--------------------------------------------------------------------------
S = 2;  % surface de Base du reservoir (m^2)

%--------------------------------------------------------------------------
% Definition des conditions initiales 
%--------------------------------------------------------------------------
h0=2.25; % hauteur initiale du liquide [m]

%--------------------------------------------------------------------------
% Debit entrant dans le reservoir 
%--------------------------------------------------------------------------
Qin=0;   %[m^3/s]

%--------------------------------------------------------------------------
% L'equation dynamique doit etre une fonction sous la forme dx/dt = f(t,x)
% Etats: x = h, Entree: u = Qin
%--------------------------------------------------------------------------
mydyn = @(t,h) (Qin-0.01*sqrt(h))/S; % h est l'etat du systeme

%--------------------------------------------------------------------------
% Integration Numerique de l'equation dynamique
%--------------------------------------------------------------------------
odeoptions = odeset('RelTol',1e-6,'AbsTol',1e-6,'Stats','on'); 
tspan = [0,600]; %horizon temporel de la simulation (s)
[tsim, hsim] = ode45(mydyn,tspan,h0,odeoptions); 

%--------------------------------------------------------------------------
% Affichage des resultats
%--------------------------------------------------------------------------
figure('Name','Reponse du Systeme');
% Calcul du temps pour vider le reservoir de moitie
thdemi=2*S/0.01*(h0)^(1/2)*(1-2^(-1/2));

plot(tsim,real(hsim),'b',[thdemi thdemi], [0 h0], 'r', [0 thdemi],[h0 h0], 'r',...
    [0 thdemi], [h0/2 h0/2], 'r')
xlabel('Temps t [s]', 'fontsize', 12);
ylabel('Hauteur h [m]', 'fontsize', 12);



