function [sys,x0] = tank1(t,x,u,flag)
%
% Simulink interface to ph_model.m     
%
% Inputs:    t    - time [s]
%            X    - State, X(1) = c, X(2) = V
%
%            U(1) - reagent flow, qr [l/s]
%            U(2) - Excess H+ in reagent, cr [mol/l]
%            U(3) - Influent flow rate, qi [l/s]
%            U(4) - Excess H+ in influent, ci [mol/l]
%
% Outputs:   x0:    Column vector of initial conditions
%            sys:   The main vector of results requested by Simulink.
%                   Will hold different information, depending on the flag.
%                   
%                   When flag=0:
%                       x0 contains initial conditions
%                       sys = [a, b, c, d, e, f, g]
%                       a = # of continuous time states
%                       b = # of discrete time states
%                       c = # of outputs
%                       d = # of inputs
%                       e = required to be 0
%                       f = 0(no) or 1(yes) for direct algebraic feed
%                           through of input to output.
%                       g = # of sample times (for continuous process, we
%                           set this equal to 1)
%
%                   When flag = 1:
%                       sys = column vector of the derivatives
%
%                   When flag = 3:
%                       sys  = column vector of the output variables:
%                       y(1): effluent concentration c in terms of excess H+
%                       y(2): tank volume V
%
%                   Else: sys is a null vector


if abs(flag) == 1
  % Return state derivatives.
  sys = ph_model(t,x,u);
  
elseif abs(flag) == 3
  % Return system outputs.
  sys(1,1) = x(1);          % Outlet concentration, c [mol/l]
  sys(2,1) = x(2);          % Tank volume, V [l]

elseif flag == 0
  % Initialize the system
  x0 = [ 0.001 ; 1000*2*40.7/3];         % col.vector of init.values of c and V
  sys = [2, 0, 2, 5, 0, 0];
  
else
  sys = [];  

end
  




