% *****           Example 3.6 - Tank Draining               *****
% *****                  from Marlin 1994                   *****
% COPYRIGHT McMASTER UNIVERSITY, 1994
% WRITTEN BY : FRANCOIS BOUDREAU
% FILENAME : EXAM3_6.M
% LAST UPDATED : NOVEMBER 1994
% VERSION : 2.0
% DESCRIPTION OF FILE : TANK DRAINING SIMULATION, EXAMPLE 3.6.
%    SIMULATES A NON-LINEAR & LINEARIZED.
% ***************************************************************
clear;                       % CLEARS WORKSPACE
clf;                         % CLEARS GRAPH
clc;                         % CLEARS SCREEN
axis('auto');

% ***************************************************************
% PROCESS DATA
% ***************************************************************
A = 7.0;                     % CROSS-SECTIONAL AREA (m)
Linit = 7.0;                 % INITIAL LEVEL (m)
Finit = 100.0;               % INITIAL IN/OUTPUT FLOWS (m^3/h)
kF1 = Finit/sqrt(Linit);     % CONST. OF RESTRICTION (m^(5/2)/h)
                             % FOR NON-LINEAR SYSTEM
tau = A/0.5/kF1*sqrt(Linit); % LINEARIZED PROCESS TIME CONST. (h)
Kp  =  (Linit^0.5)/(0.5*kF1);% LINEARIZED PROCESS GAIN

% ***************************************************************
% SIMULATION PARAMETERS
% ***************************************************************
tstart = 0.0;                % TIME TO BEGIN SIMULATION (min)
tstep = 0.5;                 % TIME TO INTRODUCE STEP (min)
step = -60.0;                % STEP IN INLET FLOW (m^3/h)
tend = 5.0;                  % TIME THE SIMULATION END (min)
delt = 0.005;                % EULER INTEGRATION STEP SIZE (min)
n = round((tend-tstart)/delt)+1;
                             % TOTAL NO. OF DATA POINTS
                             % (INCLUDING START AND END POINTS)
t = tstart:delt:tend;        % TIME VECTOR USED FOR PLOTTING

% ***************************************************************
% INITIAL CONDITIONS AND STORAGE VECTORS
% ***************************************************************
% INITIALIZATION OF STORAGE VECTORS:
Lnlin = zeros(1,n);          % LEVEL RESPONSE, N-L SYSTEM
Llind = zeros(1,n);          % LEVEL RESPONSE, LINEARIZED
F0 = zeros(1,n);             % INLET FLOW

Lnlin(1) = Linit;            % INITIALIZATION OF LEVEL
F0(1) = Finit;               % INITIALIZATION OF INLET FLOW

% ***************************************************************
% SIMULATION
% ***************************************************************
expat = exp(-delt/tau);      % CONSTANT USED IN ANALYTICAL SOL'N

for k = 1:n-1                % LOOP ON TIME INDEX
   F0(k) = Finit;            % SET INLET FLOW
   if t(k) >= tstep          % INTRODUCE STEP IF TIME>=TSTEP
      F0(k) = F0(k) + step;
   end

   % ANALYTICAL SOLUTION (LINEARIZED SYSTEM)
   Llind(k+1) = Linit;
   if t(k) >= tstep
       Llind(k+1) = Linit + step*Kp*(1-exp(-(t(k)-tstep)));
   end

   % EULER INTEGRATION (N-L SYSTEM)
   % CALCULATE DERIVATIVE:
   Lnlindot = -kF1*sqrt(Lnlin(k))/A + F0(k)/A;
   % EULER FORMULA:
   Lnlin(k+1) = Lnlin(k) + Lnlindot*delt;
end

% SET INPUT TO LAST SAMPLING TIME:
F0(n) = F0(n-1);

% ***************************************************************
% PLOTTING DATA
% ***************************************************************

subplot(211); plot(t,Llind,'-',t,Lnlin,'--');
axis([0 tend -2 8]);
xlabel('Time (h)');
ylabel('Level (m)');
title('Example 3.6,  solid: linearized;  dashed: non-linear');

subplot(212); plot(t,F0,'-');
axis([0 tend 0 110]);
xlabel('Time (h)');
ylabel('Inlet flow (m3/h)');

figure (1) ; pause ; close all