function xprime1 = colamod1(t,X1,U1,start)

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

% Ternary distillation model for HP column in multi-effect system
% H.K. Engelien, 14.03.2001, edited by Magnus H. Nilsen 18/01/2002
% edited H. K. Engelien 06/02/2002
%*******************************************************************

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

% Ternary separation
% liquid flows calculated from Francis Weir formula
% Liquid activity coefficient from Wilson equation
% Ideal gas
% Vapour pressure calculated from Antoine equation
% Variable relative volatility
% Vapour flows modelled using valve equation
% Hold-up of vapour is neglected

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

%   t -     Time (h)
%   X -     State : the first NT states are compositions of  
%           component A with reboiler/btm stage as X(1) and condenser
%           stage as X(NT), States (NT+1) to 2NT are compositions of component B
%           with reboiler/btm stage as X(NT+1) and condenser stage as X(2NT).
%           States (2NT+1) to 3NT is holdup, state (2NT+1) is the holdup 
%           in the reboiler/btm and state (3NT) is holdup in condenser.
%           States (3NT+1) to (4NT) are temperatures.
%
%   U1(1)  - reflux, LT (kmol/h) in HP column
%   U1(2)  - boilup, VB (kmol/h) in HP column
%   U1(3)  - top or distillate product flow, D (kmol/h) in HP column
%   U1(4)  - btm product flow, B (kmol/h) in HP column
%   U1(5)  - feedrate, F (kmol/h) to HP column
%   U1(6)  - feed composition, zAF to HP column
%   U1(7)  - feed composition, zBF to HP column
%   U1(8)  - feed composition, zCF to HP column
%   U1(9)  - feed liq. fraction, qF to HP column
%   U1(10) - flow condensated in HP column (kmol/h)
%   U1(11) - Temperature in feed to HP column

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

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

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

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

global P_HP

% Components :      Methanol/Water/Ethanol

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

% Nominal liquid holdups
%M0(1) = 100;               % Nominal reboiler holdup (kmol)  % These nominal values are only for use with 
%i=2:NT-1;                                                    % linearized liquid dynamics, and must be changed
%M0(i)=100*ones(1,NT-2);    % Nominal tray holdup (kmol)      % for the process
%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; %4320                 % Nominal feedrate (kmol/h)
qF0 = 1;                    % Nominal liquid fraction in feed
L0 = start(4);                 % 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 = start(10);                % Nominal vap. flow (only needed % if lambda is nonzero)
V0t = V0 + (1-qF0)*F0;      % Nominal vapour flow above feed (kmol/h)
TF_HP = 406;                % Feed temperature(K)

% Valve equation constant
k = 0.1;%0.125;                  % Value calculated from HYSYS data

% Physical data (taken @ 8 bar, from HYSYS)
TBA = 401.25;               % Boiling point, comp. A (K)
TBB = 443.57;               % Boiling point, comp. B (K)
TBC = 414.92;               % Boiling point, comp. C (K)				 
CpA = 134.9;                % Heat capacity, comp. A( kJ/kmol*K)    Liquid, at boiling point
CpB = 79.056;               % Heat capacity, comp. B (kJ/kmol*K)    Liquid, at boiling point
CpC = 188.28;               % Heat capacity, comp. C (kJ/kmol*K)    Liquid, at boiling point
HvapA = 30384;              % Hvap for component A (kJ/kmol)
HvapB = 36876;              % Hvap for component B (kJ/kmol)
HvapC = 32396;			    % Hvap for component C (kJ/kmol)
p0A = 8e5;                  % Vap. pressure for pure liq. comp. A (Pa)
p0B = 8e5;                  % Vap. pressure for pure liq. comp. B( Pa)
p0C = 8e5;                  % Vap. pressure for pure liq. comp. C (Pa)
R = 8.314;                  % Universal gas constant (kJ/kmol*K)


% Splitting the states
xA = X1(1:NT)';  			% Liquid composition of A from btm to top
xB = X1(NT+1:2*NT)';		% Liquid composition of B from btm to top
M  = X1(2*NT+1:3*NT)';      % Liquid holdup from btm to top
T  = X1(3*NT+1:4*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 
zAF = U1(6);                % Feed composition, A
zBF = U1(7);  			    % Feed composition, B
zCF = U1(8);  			    % Feed composition, C
qF  = U1(9);                % Feed liquid fraction
VD  = U1(10);               % Flow condensated
%TF_HP = U1(11);	        % Temperature in feed

i = 1:NT;
xC(i) = 1 - xA(i) - xB(i);

% Heat capacity (kJ/kmol K)
CpAA = 21.152;CpAB=70.924e-3;CpAC=25.870e-6;CpAD=-2.852e-8;   % Ref. C&R, Vol. 6
CpBA = 32.243;CpBB=19.238e-4;CpBC=10.555e-6;CpBD=-3.596e-9;   % Ref. C&R, Vol. 6
CpCA = 9.014;CpCB=21.407e-2;CpCC=-8.390e-5;CpCD=13.733e-10;   % Ref. C&R, Vol. 6

i=1:NT;
CpA = CpAA+CpAB.*T(i)+CpAC.*T(i).^2+CpAD.*T(i).^3;
CpB = CpBA+CpBB.*T(i)+CpBC.*T(i).^2+CpBD.*T(i).^3;
CpC = CpCA+CpCB.*T(i)+CpCC.*T(i).^2+CpCD.*T(i).^3;
Cp = xA(i).*CpA + xB(i).*CpB + (1-xA(i)-xB(i)).*CpC;  % Mixture heat capacity (kJ/kmol K)

% Heat of vapourisation (kJ/kmol)
i =1:NT;
Hvap(i) = xA(i)*HvapA + xB(i)*HvapB + xC(i)*HvapC;

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

% ******************************************************************************
% ************ Program written by Stathis Skouras, September 2001 **************
% ****************************************************************************** 

% Program for ternary VLE by using Wilson equation
%*******************************************************************************
% Ternary system
%                (1): Methanol
%                (2): Water
%                (3): Ethanol
%*******************************************************************************
% Pressure is 1 atm = 1.013 bar = 760 mmHg
% Composition vectors xA, xB, are line vectors of NT columns
%*******************************************************************************

%*******************************************************************************
% Data that have to be changed for different mixtures

% Pure component constants for Antoine equation (Gmehling et al, 1977-1990)
A1=8.08097; B1=1582.271; C1=239.726;  % (T-range: 15C - 84C, Component 1)
A2=8.07131; B2=1730.630; C2=233.426;  % (T-range: 1C - 100C, Component 2)
A3=8.11220; B3=1592.864; C3=226.184;  % (T-range: 20C - 93C, Component 3)

% Partial molar volumes of pure components (Gmehling et al, 1977-1990)
V1=40.73;     % Component 1 (cm^3/mol)
V2=18.07;     % Component 2 (cm^3/mol)
V3=58.68;     % Component 3 (cm^3/mol)

% Binary interaction parameters Aij for Wilson equation (Gmehling et al. 1977-1990)
A11=0;                       % Pure component 1 (cal/mol)
A22=0;                       % Pure component 2 (cal/mol)
A33=0;                       % Pure component 3 (cal/mol)
A12=-102.7135;%19.2547;                 % Binary mixture 1-2 (cal/mol) 
A21=695.4909;%554.0494;                % Binary mixture 1-2 (cal/mol) 
A13=-240.152;               % Binary mixture 1-3 (cal/mol) 
A31=637.320;               % Binary mixture 1-3 (cal/mol) 
A23=478.0720;%904.717;                 % Binary mixture 2-3 (cal/mol) 
A32=875.1566;%447.171;               % Binary mixture 2-3 (cal/mol) 

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

% Gas constant
R=1.98721;     % cal/(mol*K)
i=1:NT;

% Binary interaction parameters for Wilson equation 
L11=V1/V1*exp(-A11./(R*T));            % Pure component 1 
L22=V2/V2*exp(-A22./(R*T));            % Pure component 2 
L33=V3/V3*exp(-A33./(R*T));            % Pure component 3 
L12=V2/V1*exp(-A12./(R*T));            % Binary mixture 1-2 
L21=V1/V2*exp(-A21./(R*T));            % Binary mixture 1-2 
L13=V3/V1*exp(-A13./(R*T));            % Binary mixture 1-3  
L31=V1/V3*exp(-A31./(R*T));            % Binary mixture 1-3  
L23=V3/V2*exp(-A23./(R*T));            % Binary mixture 2-3  
L32=V2/V3*exp(-A32./(R*T));            % Binary mixture 2-3 

%Vapor pressures by Antoine equation
P1s=10.^(A1-B1./(T-273.15+C1))*(1e5/760);       % Component 1 (Pa)
P2s=10.^(A2-B2./(T-273.15+C2))*(1e5/760);       % Component 2 (Pa)
P3s=10.^(A3-B3./(T-273.15+C3))*(1e5/760);       % Component 3 (Pa)

% Activity coefficients by Wilson equation
G1(i) = exp(1-log(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xA(i).*L11(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L21(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L31(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 1

G2(i) = exp(1-log(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...
          - xA(i).*L12(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L22(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L32(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 2

G3(i) = exp(1-log(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i))...
          - xA(i).*L13(i)./(xA(i).*L11(i)+xB(i).*L12(i)+(1-xA(i)-xB(i)).*L13(i))...
          - xB(i).*L23(i)./(xA(i).*L21(i)+xB(i).*L22(i)+(1-xA(i)-xB(i)).*L23(i))...    
          - (1-xA(i)-xB(i)).*L33(i)./(xA(i).*L31(i)+xB(i).*L32(i)+(1-xA(i)-xB(i)).*L33(i)));    % Component 3


% Total pressure in each stage (bar)
P(i) = (xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i)+(1-xA(i)-xB(i)).*G3(i).*P3s(i));
P_HP = P;

% Vapor Liquid Equilibria (VLE)
yA(i)=xA(i).*G1(i).*P1s(i)./(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));            % Component 1
yB(i)=xB(i).*G2(i).*P2s(i)./(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));            % Component 2
yC(i)=(1-xA(i)-xB(i)).*G3(i).*P3s(i)./(xA(i).*G1(i).*P1s(i) + xB(i).*G2(i).*P2s(i) + (1-xA(i)-xB(i)).*G3(i).*P3s(i));  % Component 3


% 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)

Muw_af = 14.4;                  % Liq. holdup under weir (kmol/h) above feed
M0_af = 28.8;                  % Nominal liquid holdup above feed
Muw_bf = 17.75;                  % Liq. holdup under weir (kmol/h) below feed
M0_bf = 35.5;                  % Nominal liquid holdup below feed
Kaf = L0/((M0_af-Muw_af)^1.5);            % Constant K above feed
Kbf = (L0+F)/((M0_bf-Muw_bf)^1.5);        % Constant k below feed

i=2:NF;      L(i)= Kbf*(M(i)-Muw_bf).^1.5;     % Liquid flows below feed (Kmol/h)    
i=NF+1:NT-1; L(i)= Kaf*(M(i)-Muw_af).^1.5;     % Liquid flows above feed (Kmol/h) 
L(NT)=LT;                                   % Condenser's liquid flow (Kmol/h)  

% 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;

%--------------------------------------------------------------------
% 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);                                 
dMxAdt(i) = L(i+1).*xA(i+1) - L(i).*xA(i) + V(i-1).*yA(i-1) - V(i).*yA(i);
dMxBdt(i) = L(i+1).*xB(i+1) - L(i).*xB(i) + V(i-1).*yB(i-1) - V(i).*yB(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;
dMxAdt(NF) = dMxAdt(NF) + F*zAF;
dMxBdt(NF) = dMxBdt(NF) + F*zBF;
MdTdt(NF)  = MdTdt(NF)  + F*((1-qF)*Hvap(NF)/Cp(NF));

% note, the feed energy balance assumes temperature of feed is the same as 
% the temperature of feed stage, otherwise should include + (TF_HP - T(NF))
% in the last term: F*((1-qF)*Hvap(NF)/Cp+ (TF_HP - T(NF)))

% Reboiler (assumed to be an equilibrium stage)
dMdt(1)   = L(2) - V(1) - B;
dMxAdt(1) = L(2)*xA(2) - V(1)*yA(1) - B*xA(1);
dMxBdt(1) = L(2)*xB(2) - V(1)*yB(1) - B*xB(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;
dMxAdt(NT) = V(NT-1)*yA(NT-1) - LT*xA(NT) - D*xA(NT);
dMxBdt(NT) = V(NT-1)*yB(NT-1) - LT*xB(NT) - D*xB(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;
dxAdt(i) = (dMxAdt(i) - xA(i).*dMdt(i))./M(i);
dxBdt(i) = (dMxBdt(i) - xB(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 = [dxAdt';dxBdt';dMdt';dTdt'];     





   
