function xprime=multicomAE_mat(t,X,U) 

% *************************************************************************************
% ***************** Program written by Stathis Skouras at March, 2001 *****************  
% *************************************************************************************

% multicomA- This is a nonlinear model of a continuous distillation column with
%            NT-1 theoretical stages including a reboiler (stage 1) plus a
%            total condenser ("stage" NT). 
%            Model assumptions: 
%            NC components (multicomponent mixture); Component NC is the heavy component
%            constant relative volatilities; 
%            no vapor holdup; 
%            one feed and two products;
%            constant molar flows (same vapor flow on all stages); 
%            Liquid flow dynamics modelled by Franci's Weir Formula.
%            total condenser.  
%
%            The model is based on column A in Skogestad and Postlethwaite
%            (1996). The model has NC*NT states.
%
% Inputs:    t    - time in [min].
%            X    - State, the first NT states are compositions of light
%                   component A, the next NT states are compositions of component B, etc,
%                   compositions of the heavy component are not reported directly  
%                   reboiler/bottom stage is stage (1) and condenser is stage (41). 
%                   The last NT states are liquid holdups in each stage. 
%            U(1) - reflux L,
%            U(2) - boilup V,
%            U(3) - top or distillate product flow D,
%            U(4) - bottom product flow B,
%            U(5) - feed rate F,
%            U(6) - feed liquid fraction, qF.
%             A   - feed compositions, zF.
%         U=[U A] - Disturbances vector
% Outputs:   xprime - vector with time derivative of all the states 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

% The following data need to be changed for a new column.
% These data are for "column A".

% Number of stages (including reboiler and total condenser: 
NT=32;
% Location of feed stage (stages are counted from the bottom):
NF=20;
% Number of components
% C02,C1,C2,C2,C3,iC4,nC4,iC5,nC5,C6,C7
% The total number of component was reduced to 5 by regrouping CO2,C1, C2
% in C2 and iC5,nC5,C6,C7 in C5
NC=5;
%% Relative volatilities: alpha=[alpha1 alpha2...alphaNC]
%%
%alpha=[2.98 2.177419355 1.9 1.223809524  1];%(C2/C3,C3/iC4,nC4/C5,iC4/nC4,C5+/C5+)
%alpha=[21.28 7.14 3.29 2.17  1];%(C2/iC5,C3/iC5,iC4/iC5,nC4/iC5,iC5/iC5)
alpha=[2.98 1 0.46 0.304 0.14]; % alpha=Ki/Kr where Kr is C3 and Ki any other component
% Splitting the states
x=X(1:(NC-1)*NT);                       % Liquid compositions from btm to top
M=X((NC-1)*NT+1:NC*NT);                 % Liquid hold up from btm to top
% Inputs and disturbances
% Inputs
LT = U(1);                             % Reflux
VB = U(2);                             % Boilup
D  = U(3);                             % Distillate
B  = U(4);                             % Bottoms
% Disturbances
F = U(5);                             % Feedrate
qF= U(6);                             % Feed liquid fraction
zF= U(7:size(U,2));                   % Feed composition
%zF= U(7:6+NC);                       % Feed composition

% Rearrange elements of composition vector (x) for later use
Iu=[1:NT]'*ones(1,NC-1)+NT*ones(NT,1)*[0:NC-2];    
x=(x(Iu))';                         


% THE MODEL

% Vapour-liquid equilibria (multicomponent ideal VLE, Stichlmair-Fair, 'Distillation', p. 36, 1998)
y=(alpha(1:NC-1)'*ones(1,NT).*x)./(ones(1,NC-1)'*(1+(alpha(1:NC-1)-1)*x));


% Vapor Flows assuming constant molar flows
i=1:NT-1;    V(i)=VB*ones(1,NT-1);
i=NF:NT-1;   V(i)=V(i) + (1-qF)*F;


% Liquid flows are given by Franci's Weir Formula L(i)=K*Mow(i)^1.5 
% Liquid flow L(i) dependent only on the holdup over the weir Mow(i)
% M(i)= Mow(i) + Muw(i) (Total holdup = holdup over weir + holdup below weir) 
                              
%% using Excel calculations "Elviras"
% Muw=pi*(d/2)^2*weir*rho/MW
% Kbf=L*(M(i)-Muw)^(2/3
%%
% Kuf=71.3734;                                    % Constant above feed min^-1
% Kbf=189.877;                                    % Constant below feed min^-1
% Muw=3.97;                                     % Liquid holdup under weir % (Kmol)
% Mow=3.18;
% i=2:NF;      L(i)= Kbf*(M(i)-Muw).^1.5;      % Liquid flows below feed (Kmol/min)    
% i=NF+1:NT-1; L(i)= Kuf*(M(i)-Muw).^1.5 ;      % Liquid flows above feed (Kmol/min) 
% i=NF+1:NT-1; L(i)= Kuf*(M(i)-Muw).^1.5 ;      % Liquid flows above feed (Kmol/min) 
% L(NT)=LT;                                       % Condenser's liquid flow (Kmol/min)  
%% using Excel calculations "Elviras" Methode 2
% Muw=(A_bubble+A_downcomerCenter)*weir*rho/MW
% M=(A_bubble+A_downcomerCenter)*level*rho/MW
% Kbf=L*(M(i)-Muw)^(2/3
%%
Kuf=272.9029964;                                    % Constant above feed min^-1
Kbf=87.0985927;                                    % Constant below feed min^-1
Muw=3.197907135;                                     % Liquid holdup under weir % (Kmol)
Mow=2.333537198;                                     % Liquid holdup under weir % (Kmol) 
i=2:NF;      L(i)= Kbf*(M(i)-Muw).^1.5;      % Liquid flows below feed (Kmol/min)    
i=NF+1:NT-1; L(i)= Kuf*(M(i)-Mow).^1.5 ;      % Liquid flows above feed (Kmol/min) 
L(NT)=LT;                                       % Condenser's liquid flow (Kmol/min) 
% Time derivatives from material balances for 
% 1) total holdup and 2) component holdup

% Column
j=2:NT-1;
dMdt(j) = L(j+1)         - L(j)       + V(j-1)         - V(j);

for i=1:NC-1;
for j=2:NT-1;
dMxdt(i,j) = L(j+1)*x(i,j+1)  - L(j).*x(i,j)   + V(j-1)*y(i,j-1)     - V(j)*y(i,j);
end 
end

% Correction for feed at the feed stage
% The feed is assumed to be mixed into the feed stage

dMdt(NF) = dMdt(NF)  + F;

dMxdt(:,NF)=dMxdt(:,NF)+F*zF(1:NC-1)';

% Reboiler (assumed to be an equilibrium stage)
dMdt(1) = L(2)          - V(1)        - B;

i=1:NC-1;
dMxdt(i,1)= L(2)*x(i,2) - V(1)*y(i,1) - B*x(i,1);

% Total condenser (no equilibrium stage)
dMdt(NT) = V(NT-1)           - LT         - D;

i=1:NC-1;
dMxdt(i,NT)= V(NT-1)*y(i,NT-1) - LT*x(i,NT) - D*x(i,NT);

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
dxdt=(dMxdt-x.*(ones(NC-1,1)*dMdt))./(ones(NC-1,1)*M');

% Rearrange elements of composition vector (dxdt) for later use
Ix=[1:(NC-1)*NT]';   
W=dxdt';
dxdt=W(Ix);                        

% Output
xprime=[dxdt;dMdt'];

