function xprime1 = colamod1(t,X1,U1)

%*******************************************************************
%                       DISTILLATION MODEL 
%
% Binary distillation model for column with pressure variations
%
% H.K. Engelien, 23.10.2001
%*******************************************************************

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

% Binary separation
% LLiquid flows - Francis weir formula (or linearised liq. flow dynamics)
% Liquid activity coefficient from Wilson equation
% Ideal gas
% Vapour pressure calculated from Antoine equation
% Vapour flows modelled using valve equation
% Hold-up of vapour is neglected
% Heat of vapourisation from equation of Haggenmacher
% Heat capacity - linear relationship with temp. for pure components 
%
%--------------------------------------------------------------------
%                           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)
%   U1(9) - feed temperature (K)
%
% Outputs : xprime1 - vector with time derivatives of all the states

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

% NB : The following data should be altered for a new column:
% (+ relationship for heat capacity, Cp)
%**********************************************************************
%                           HP Column 
%**********************************************************************
global P_HP Hvap

% 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/60;           % Time constant for liquid dynamics (h)
%F0  = 4320;               % Nominal feedrate (kmol/h)
%qF0 = 1;                  % Nominal liquid fraction in feed
%L0  = 5747;               % Nominal reflux (from steady state data)
%L0b = L0 + qF0*F0;        % Nominal liq. flow below feed (kmol/h)
%lambda = 0;               % Effect of vap.flow on liq.flow (K2 effect)
%V0  = 7522;               % Nominal vap. flow (only needed 
                           % if lambda is nonzero)
%V0t = V0 + (1-qF0)*F0;    % Nominal vapour flow above feed (kmol/h)

% Data for Francis Weir formula (from steady state)
F0  = 4320;                % Nominal feedrate (kmol/h)
L0  = 5747;                % Nominal reflux (kmol/h)
Muw = 50;                  % Liq. holdup under weir (kmol/h)
M0 = 100;                  % Nominal liquid holdup

% Valve equation constant
k = 0.121;                 % Value eq. constant (calculated based on HYSYS data)

% Physical data (taken @ 5.79 bar, from HYSYS)
PcL = 81e5;                % Critical pressure of light component (Pa)
PcH = 220.5e5;             % Critical pressure of heavy component (Pa)
TcL = 512.6;               % Critical temperature of light component (K)
TcH = 647.3;               % Critical temperature of heavy component (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 

% Antoine Coefficients
A1 = 18.5875; B1 = 3626.55; C1 = -34.29;
A2 = 18.3036; B2 = 3816.44; C2 = -46.13;


%---------------------------------------------------------------------
% 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
TF = U1(9);                 % Feed temperature (K)

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

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

i=1:NT;
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 pressure fra Antoine equation (Pa)
i=1:NT;
pL(i) = exp(A1 - B1./(T(i)+C1))*133.32;           
pH(i) = exp(A2 - B2./(T(i)+C2))*133.32;           
P(i)  =  x(i).*pL(i).*gammaL(i) + (1-x(i)).*pH(i).*gammaH(i);
P_HP  = P;                      
P_HPtop = P_HP(NT);

% Heat of Vapourisation (from Haggenmacher)
i=1:NT;
PrL(i) = P(i)./PcL; PrH = P(i)./PcH;
TrL(i) = T(i)./TcL; TrH = T(i)./TcH;
Delta_zL(i) = (1-PrL./(TrL.^3)).^0.5; Delta_zH = (1-PrH./(TrH.^3)).^0.5;
HvapL(i) = 8.32*B1.*T(i).^2.*Delta_zL(i)./((T(i)+C1).^2);
HvapH(i) = 8.32*B2.*T(i).^2.*Delta_zH(i)./((T(i)+C2).^2);

Hvap(i) = x(i).*HvapL(i) + (1-x(i)).*HvapH(i);   % Heat of vapourisation (kJ/kmol)

% Heat capacity (Cp) - take pure Cp as a linear fn. of temperature
% - the linear relationship is found using HYSYS at different temperatures (pure component)
i=1:NT;
CpL(i) = 0.3555.*T(i) + 91.583;                  % Heat capacity light component
CpH(i) = 0.0394.*T(i) + 71.835;                  % Heat capacity light component
Cp(i)  = x(i).*CpL(i) + (1-x(i)).*CpH(i);        % Heat capacity (kJ/kmol K)

% Vapour liquid equilibria
i=1:NT-1;
y(i) = gammaL(i).*x(i).*pL(i)./P(i);             % Vap. fraction 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;

% Liquid flows - modelled Franci's Weir Formula: L(i) = K*Mow(i)^1.5
% Liquid flow L(i) is dependent only on the holdup over the weir Mow(i)
% M(i) = Mow(i) + Muw(i)  (Total holdup = holdup over weir + holdup below weir)

Kaf = L0/((M0-Muw)^1.5);            % Constant K above feed
Kbf = (L0+F)/((M0-Muw)^1.5);        % Constant k below feed

i=2:NF;
L(i)  = Kbf*(M(i)-Muw).^1.5;
i=NF+1:NT-1;
L(i)  = Kaf*(M(i)-Muw).^1.5;
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(i)./Cp(i))-V(i).*Hvap(i)./Cp(i);                                                            

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


% 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(1)/Cp(1)*(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(NT)/Cp(NT)*(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'];
