
function xprime=ph_model(t,X,U) 
%
% ph_model - This is a model of a single pH neutralization tank, where an
%            acid (influent) is mixed with a base (reagent) to yield a
%            product with pH = 7. Concentrations are in terms of excess H+,
%            e.g. c = c_H+ - c_OH-
%
% Inputs:    t      - Time [s].
%            X      - States, 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:   xprime - vector with time derivative of the states
%                     (states: X(1) = c, X(2) = V)


 

% Splitting the states
c=X(1);                              % Excess H+ in effluent
V=X(2);                              % Volume of tank

% One manipulated input:
qr = U(1);                            % Reagent flow rate
q = U(2);

% Three disturbances:
cr = U(3);                            % Excess H+ in reagent, cr [mol/l]
qi = U(4);                            % Influent flow rate
ci = U(5);                            % Excess H+ in influent, ci [mol/l]


% Total mass balance:
dVdt = qi + qr - q;


% Component mass balance (in terms of excess H+):
%dcVdt = ci*qi + cr*qr - c*q;

dcdt = (1/V)*((ci-c)*qi + (cr-c)*qr);


% Compute the derivative for the concentration from d(cV) = c dV + V dc                        
%dcdt = (dcVdt - c*dVdt)/V;                                             

% Output
xprime = [dcdt; dVdt];      