function xprime1 = colamod1(t,X1,U1)

%*******************************************************************
% DISTILLATION MODEL (HP COLUMN)

% Binary distillation model for HP 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 : xprime1 - vector with time derivatives of all the states

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

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

%**********************************************************************
%                           HP Column (Column 1)
%**********************************************************************
global P_HP

% Components :      Methanol/Water

NT = 16;    % Number of stages (including reboiler and total condenser)
NF = 4;     % 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:
tau = 0.063;               % Time constant for liquid dynamics (min)
F0  = 4320;                % Nominal feedrate (kmol/min)
qF0 = 1;                   % Nominal liquid fraction in feed
L0  = 778;                 % 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  = 2513;                % 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 @ 5.2 bar, from HYSYS)
TBL = 385.85;              % Boiling point light component (K)
TBH = 426.45;              % Boiling point heavy component (K)
CpL = 130.2;               % Heat capacity light component (kJ/kmol*K)
CpH = 77.91;               % Heat capacity heavy component (kJ/kmol*K)
HvapL = 31630;             % Hvap for light component (kJ/kmol)
HvapH = 37870;             % Hvap for heavy component (kJ/kmol)
p0L = 5.2e5;               % Vap. pressure for pure light liq. comp.(Pa)
p0H = 5.2e5;               % Vap. pressure for pure heavy liq. comp.(Pa)
R   = 8.314;               % Universal gas constant (kJ/kmol*K)

% Wilson equation coefficients :
% (from Gmehling & Onken VLE Data Collection, p.51)
    
A11 = 1;
A12 = -102.7135;           % Reference pressure : 5.066 bar
A21 = 695.4902;            % Reference pressure : 5.066 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 = X1(1:NT)';              % Liquid composition from btm to top
M = X1(NT+1:2*NT)';         % Liquid holdup from btm to top
T = X1(2*NT+1:3*NT)';       % Temperature

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

F =  U1(5);                 % Feedrate
zF = U1(6);                 % Feed composition
qF = U1(7);                 % Feed liquid fraction
VD = U1(8);                 % Flow condensated

% Assume linearity between Hvap, Cp and composition

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_HP  = P;
P_HPtop = P_HP(NT);

% Activity coefficients from Wilson equation
R = 1.98721;                  % Universal gas constant (cal/kmol*K)

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. fracn. 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

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