function xprime=colamod_nonideal(t,X,U) 

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

% colamod_nonideal - 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). The liquid flow dynamics are
%           modelled by Franci's Weir Formula.
%           Model assumptions: Three components (ternary mixture); constant
%           relative volatility; no vapor holdup; one feed and two products;
%           constant molar flows (same vapor flow on all stages); 
%           total condenser
%
%           The model is based on column A in Skogestad and Postlethwaite
%           (1996). 
%
% Inputs:    t    - time in [min].
%            X    - States: (NC-1)*NT mole fractions, NT liquid holdups, NT temperatures
%            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) - U(8) - feed composition, zFA, zFB, zFC.
%            U(9) - feed liquid fraction, qF.
%
% Outputs:   xprime - vector with time derivative of all the states 

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

%------------------------------------------------------------
% START: 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=41; 
% Location of feed stage (stages are counted from the bottom):
    NF=21;
% END: The following data need to be changed for a new column.
    
% Extract the states
xA=X(1:NT)';                           % Liquid composition of A from btm to top
xB=X(NT+1:2*NT)';                      % Liquid composition of B from btm to top
M=X(2*NT+1:3*NT)';                     % Liquid hold up from btm to top
T=X(3*NT+1:4*NT);                      % Temperature from btm to top 

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

F = U(5);                              % Feedrate
zFA =U(6);                             % Feed composition of component A
zFB =U(7);                             % Feed composition of component B
zFC =U(8);                             % Feed composition of component C
qF = U(9);                             % Feed liquid fraction


% THE MODEL

% Vapor-liquid equilibria by Wilson or UNIQUAC equation (your choice!)

Y=wilson(xA,xB,T,NT);                  % Calls 'Wilson' VLE model 
% Y=uniquac(xA,xB,T,NT);                 % Calls 'UNIQUAC' VLE model

% Splitting the states
yA=Y(1,1:NT);                          % Vapor composition for component A
yB=Y(2,1:NT);                          % Vapor composition for component B 
T=Y(3,1:NT);                           % Temperature estimation in each stage
DP=Y(4,1:NT);                          % Pressure deviation at each stage (used for the bubble point calculation)


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

Kuf=21.65032;                               % Constant above feed
Kbf=29.65032;                               % Constant below feed
Muw=0.25;                                   % 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)-Muw).^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
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);

% Correction for feed at the feed stage
% The feed is assumed to be mixed into the feed stage
dMdt(NF) = dMdt(NF)  + F;
dMxAdt(NF)= dMxAdt(NF) + F*zFA;
dMxBdt(NF)= dMxBdt(NF) + F*zFB;


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

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

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
i=1:NT;
dxAdt(i) = (dMxAdt(i) - xA(i).*dMdt(i) )./M(i);
dxBdt(i) = (dMxBdt(i) - xB(i).*dMdt(i) )./M(i);

% Output - Note the algebraic equations introduced by DP' (see mass matrix)
xprime=[dxAdt';dxBdt';dMdt';DP'];
