function xprime2 = colamod2(t,X2,U2)

%*******************************************************************
% DISTILLATION MODEL (LP COLUMN)

% Binary distillation model for LP column in multi-effect system
% H.K. Engelien, 14.03.2001
%*******************************************************************


%--------------------------------------------------------------------
%                           Assumptions
%--------------------------------------------------------------------

% Binary separation
% Linearised liquid flow dynamics
% Liquid activity coefficient from Wilson equation
% Ideal gas
% Vapour pressure calculated from Calusius-Clapeyron equation
% Variable relative volatility
% Vapour flows modelled using valve equation
% Hold-up of vapour is neglected

%--------------------------------------------------------------------
%                           Inputs
%--------------------------------------------------------------------

%   t -     time (min)
%   X -     State : the first NT states are compositions of light 
%           component with reboiler/btm stage as X(1) and condenser
%           stage as X(NT)
%           States (NT+1) to 2NT is holdup, state (NT+1) is the holdup 
%           in the reboiler/btm and state (2NT) is holdup in condenser.
%           States (2NT+1) to (3NT) are temperatures

%   U1(1) - reflux, LT (kmol/min) in HP column
%   U1(2) - boilup, VB (kmol/min) in HP column
%   U1(3) - top or distillate product flow, D (kmol/min) in HP column
%   U1(4) - btm product flow, B (kmol/min) in HP column
%   U1(5) - feedrate, F (kmol/min) to HP column
%   U1(6) - feed composition, zF to HP column
%   U1(7) - feed liq. fraction, qF to HP column
%   U1(8) - flow condensated in LP column (kmol/h)

%   U2(1) - reflux, LT (kmol/min) in LP column
%   U2(2) - boilup, VB (kmol/min) in LP column
%   U2(3) - top or distillate product flow, D (kmol/min) in LP column
%   U2(4) - btm product flow, B (kmol/min) in LP column
%   U2(5) - feedrate, F (kmol/min) to LP column
%   U2(6) - feed composition, zF, to LP column
%   U2(7) - feed liq. fraction, qF, to LP column
%   U2(8) - flow condensated in LP column (kmol/h)

% Outputs : xprime2 - vector with time derivatives of all the states

%---------------------------------------------------------------------

% NB : The following data should be altered for a new column:

%**********************************************************************
%                           LP Column (Column 2)
%**********************************************************************
global P_LP

% Components :      Methanol (L) / Water (H)

NT = 26;    % Number of stages (including reboiler and total condenser)
NF = 13;    % Location of feed stage (counted from btm)

% Nominal liquid holdups
M0(1) = 100;                % Nominal reboiler holdup (kmol)
i=2:NT-1;
M0(i) = 100*ones(1,NT-2);   % Nominal tray holdup (kmol)
M0(NT) = 100;               % Nominal condenser holdup (kmol)

% Data for linearized liquid flow dynamics (nominal rates from HYSYS):
tau = 0.063;                % Time constant for liquid dynamics (min)
F0  = 2598;                 % Nominal feedrate (kmol/min)
qF0 = 1;                    % Nominal liquid fraction in feed
L0  = 461;                  % Nominal reflux (from steady state data)
L0b = L0 + qF0*F0;          % Nominal liq. flow below feed (kmol/min)
lambda = 0;                 % Effect of vap.flow on liq.flow (K2 effect)
V0  = 1954;                 % Nominal vap. flow (only needed 
                            % if lambda is nonzero)
V0t = V0 + (1-qF0)*F0;      % Nominal vapour flow above feed (kmol/min)

% Valve equation constant
k = 0.1;                   % Value eq. constant

% Physical data (taken @ 1 bar, from HYSYS)
TBL = 337.63;               % Boiling point light component (K)
TBH = 373.15;               % Boiling point heavy component (K)
CpL = 120.2;                % Heat capacity light component (kJ/kmol*K)
CpH = 75.64;                % Heat capacity heavy component (kJ/kmol*K)
HvapL = 35290;              % Hvap for light component (kJ/kmol)
HvapH = 40650;              % Hvap for heavy component (kJ/kmol)
p0L = 1.01e5;               % Vap. pressure for pure light liq. comp.(Pa)
p0H = 1.01e5;               % Vap. pressure for pure heavy liq. comp.(Pa)
R   = 8.314;                % Univarsal gas constant (kJ/kmol*K)


% Wilson equation coefficients 
% (from Gmehling & Onken VLE Data Collection, p.47)

A11 = 1;
A12 = 82.9876;               % Reference pressure 1.013 bar
A21 = 520.6458;              % Reference pressure 1.013 bar    
A22 = 1;
V_water    = 18.07;          % Molar vol. of pure liq. component
V_methanol = 40.73;          % Molar vol. of pure liq. component 

%---------------------------------------------------------------------
% End of data that needs to be changed
% --------------------------------------------------------------------

% Splitting the states
x = X2(1:NT)';            % Liquid composition from btm to top
M = X2(NT+1:2*NT)';       % Liquid holdup from btm to top
T = X2(2*NT+1:3*NT)';     % Temperature

% Inputs and disturbances
LT = U2(1);               % Reflux
VB = U2(2);               % Boilup
D  =  U2(3);              % Distillate
B  =  U2(4);              % Bottoms

F  =  U2(5);              % Feedrate (= btm flow from HP col.)
zF = U2(6);               % Feed composition
qF = U2(7);               % Feed liquid fraction
VD = U2(8);               % Flow condensated

% Assume linearity between Hvap, Cp and composition (check this...)

Cp   = zF*CpL + (1-zF)*CpH;         % Heat capacity for feed
Hvap = zF*HvapL + (1-zF)*HvapH;     % Heat of vap. for feed

%**********************************************************************
%                           THE MODEL
%**********************************************************************

% Vapour pressure (Calusius-Calpeyron equation)
i=1:NT;
pL(i) = p0L*exp(-HvapL/R*(1./T(i) - 1/TBL));   % Vap. pres. light comp.
pH(i) = p0H*exp(-HvapH/R*(1./T(i) - 1/TBH));   % Vap. pres. heavy comp.
P(i)  = x(i).*pL(i) + (1-x(i)).*pH(i);       % Stage pres. (Raoult's law)
P_LP  = P;

R = 1.98721;                  % Universal gas constant (cal/kmol*K)

% Wilson equation (calculating liquid activity coefficients)

i=1:NT-1;
lambda_12(i) = V_water/V_methanol*exp(-(A12-A11)./(R.*T(i)));
lambda_21(i) = V_methanol/V_water*exp(-(A21-A22)./(R.*T(i)));
gammaL(i)    = exp(-log(x(i)+lambda_12(i).*(1-x(i)))...
               +(1-x(i)).*(lambda_12(i)./(x(i)+lambda_12(i).*(1-x(i)))...
               -lambda_21(i)./(lambda_21(i).*x(i)+(1-x(i)))));
gammaH(i)    = exp(-log((1-x(i))+lambda_21(i).*x(i))...
               -x(i).*(lambda_12(i)./(x(i)+lambda_12(i).*(1-x(i)))...
               -lambda_21(i)./(lambda_21(i).*x(i)+(1-x(i)))));


% Vapour liquid equilibria
i=1:NT-1;
y(i) = gammaL(i).*x(i).*pL(i)./P(i);        % Vap. fn of light component

% Vapour flows - modelled using valve equation
i=1:NT-1;
alpha = P(i+1)-P(i) < 1e-10;
V(i)  = (alpha).*k.*sqrt(P(i).^2 - P(i+1).^2) + (1-alpha).*0.001;


% Liquid flows - modelled using linearised tray hydraulics
i=2:NF;
L(i) = L0b + (M(i)-M0(i))./tau + lambda.*(V(i-1)-V0);
i=NF+1:NT-1;
L(i)  = L0 + (M(i)-M0(i))./tau + lambda.*(V(i-1)-V0t);
L(NT) = LT;

%--------------------------------------------------------------------
% Time derivatives for :
%   1) total holdup
%   2) component holdup
%   3) energy holdup (temperature)
%--------------------------------------------------------------------

% Column
i=2:NT-1;
dMdt(i)  = L(i+1) - L(i) + V(i-1) - V(i);                                 
dMxdt(i) = L(i+1).*x(i+1) - L(i).*x(i) + V(i-1).*y(i-1) - V(i).*y(i);        
MdTdt(i) = L(i+1).*(T(i+1)-T(i)) + V(i-1).*(T(i-1)-T(i)... 
           + Hvap/Cp)-V(i).*Hvap/Cp;                                                            

% Feed stage
dMdt(NF)  = dMdt(NF) + F;
dMxdt(NF) = dMxdt(NF) + F*zF;
MdTdt(NF) = MdTdt(NF) + F*(1-qF)*Hvap/Cp;


% Reboiler (assumed to be an equilibrium stage)
dMdt(1)  = L(2) - V(1) - B;
dMxdt(1) = L(2)*x(2) - V(1)*y(1) - B*x(1);
MdTdt(1) = L(2)*(T(2) - T(1)) + Hvap/Cp*(VB-V(1));

% Total condenser (no equilibrium stage)
dMdt(NT)  = V(NT-1) - LT - D;
dMxdt(NT) = V(NT-1)*y(NT-1) - LT*x(NT) - D*x(NT);
MdTdt(NT) = V(NT-1)*(T(NT-1) - T(NT)) + Hvap/Cp*(V(NT-1) - VD);

% Compute the derivatives for the model fractions from:
%       d(Mx) = x dM + Mdx
i=1:NT;
dxdt(i) = (dMxdt(i) - x(i).*dMdt(i))./M(i);

% Compute the derivatives for the temperature from :
%      dTdt = MdTdt/M
i=1:NT;
dTdt(i) = (MdTdt(i))./M(i);

% Output - matrix containing the time derivatives for x, M and T:

xprime2 = [dxdt';dMdt';dTdt'];     
