function xprime=colmmod(t,X,flag,U,P) 
% function xprime=colmmod(t,X,flag,U,P) 
%
% colmmod - This is a nonlinear model of a 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 a simple linear relationship.
%           Model assumptions: Multi components (NC); constant
%           relative volatility; no vapor holdup; one feed and two products;
%           constant molar flows (same vapor flow on all stages); 
%           total condenser
%
%           PRIMARY PURPOSE: Dynamic simulation
%
% Change (NB): Boiler (M(1)) and condenser (M(NT)) holdup
%
%           The model is based on column A in Skogestad and Postlethwaite
%           (1996). The model has NT*NC states.
%
% Inputs:    t    - time in [min].
%            X    - State, the first NT states are holdup of light
%                   component 1 with reboiler/bottom stage as X(1) and 
%                   condenser as X(NT). State X(NT+1) is holdup in reboiler/
%                   bottom stage of component 2 and X(2*NT) is hold-up in 
%                   condenser of component 2.... and X(NT*(NC-1)+1:NT*NC)
%                   are holdup of heavy component NC in the column.
%            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 composition, zF.
%            U(7) - feed liquid fraction, qF.
%
% Outputs:   xprime - vector with time derivative of all the states 
%
% Modified by: msg20000106 Multicomponent column (NC)
%                          Extern parameters 
%              msg20000202 round to get NC, NT, NF as integer
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------------------------------------------
% Parameters
NC=round(P(1)); % Number of components
NT=round(P(2)); % Number of stages (including reboiler and total condenser: 
NF=round(P(3)); % Location of feed stage (stages are counted from the bottom):
alpha=P(4:NC+3);alpha=alpha(:); % Relative volatility
% Nominal liquid holdups
    M00=P(NC+4);                     % Nominal liquid holdups
    M0=zeros(NT,1);
    M0(1)=M00;                       % Nominal reboiler holdup (kmol)
    i=2:NT-1; M0(i)=M00*ones(NT-2,1);% Nominal stage (tray) holdups (kmol)
    M0(NT)=M00;                      % Nominal condenser holdup (kmol)
% Data for linearized liquid flow dynamics (does not apply to reboiler and condenser):
    taul=P(NC+5);     % time constant for liquid dynamics (min)
    F0=P(NC+6);       % Nominal feed rate (kmol/min) 
    qF0 =P(NC+7);     % Nominal fraction of liquid in feed 
    L0=P(NC+8);       % Nominal reflux flow (from steady-state data)
    L0b=L0 + qF0*F0;  % Nominal liquid flow below feed (kmol/min)
    lambda=P(NC+9);   % Effect of vapor flow on liquid flow ("K2-effect")
    V0=P(NC+10);
    V0t=V0+(1-qF0)*F0;% Nominal vapor flows - only needed if lambda is nonzero 
% End data which need to be changed
%------------------------------------------------------------



%Index
I0=ones(NT,1)*[1:NT:NC*NT]+[0:1:NT-1]'*ones(1,NC); %Different comp. in
                                                   %diff. rows 
Im=(1:NT)'*ones(1,NC); Im=Im(1:NC*NT)';            %Full length vector:
                                                   %NT -> NC*NT
Im2=ones(NT,1)*(1:NC); Im2=Im2(1:NC*NT)';          %Full length vector:
                                                   %NC -> NC*NT

% Splitting the states
Mi=X;                     % Component liquid hold up 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
zF = U(6);                            % Feed composition
qF = U(7);                            % Feed liquid fraction

% THE MODEL

%Total liquid holdup from btm to top
M=sum(Mi(I0),2);          % Total liquid hold up from btm to top
Mm=M(Im);

%Liquid fraction
x=Mi./Mm;

% Vapor-liquid equilibria
alpham=alpha(Im2); 
alphaxm=alpham.*x;
sumalphax=sum(alphaxm(I0),2);
sumalphaxm=sumalphax(Im);
y=alphaxm./sumalphaxm;

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

% Liquid flows assuming linearized tray hydraulics with time constant taul
% Also includes coefficient lambda for effect of vapor flow ("K2-effect").
L=zeros(NT,1);
i=2:NF;      L(i) = L0b + (M(i)-M0(i))/taul + lambda*(V(i-1)-V0);
i=NF+1:NT-1; L(i) = L0  + (M(i)-M0(i))/taul + lambda*(V(i-1)-V0t);
L(NT)=LT;
Lm=L(Im);

% Time derivatives from  material balances for component holdup 
dMxdt=zeros(NC*NT,1);

% Column
i=2:NT-1; Ic=I0(2:NT-1,:); Ic=Ic(1:NC*(NT-2))';Icp=I0(3:NT,:); ...
    Icp=Icp(1:NC*(NT-2))'; Icm=I0(1:NT-2,:); Icm=Icm(1:NC*(NT-2))';
dMxdt(Ic)=Lm(Icp).*x(Icp)-Lm(Ic).*x(Ic)+Vm(Icm).*y(Icm)-Vm(Ic).*y(Ic);


% Correction for feed at the feed stage
% The feed is assumed to be mixed into the feed stage
i=NF; Ic=I0(i,:); Ic=Ic(1:NC)';
zF=[zF;1-zF];
dMxdt(Ic)= dMxdt(Ic) + F*zF;

 
% Reboiler (assumed to be an equilibrium stage)
i=1; Ic=I0(i,:); Ic=Ic(1:NC)'; Ic2=I0(2,:); Ic2=Ic2(1:NC)';
dMxdt(Ic)= L(2)*x(Ic2) - V(1)*y(Ic) - B*x(Ic);


% Total condenser (no equilibrium stage)
i=NT; Ic=I0(i,:); Ic=Ic(1:NC)'; Ic2=I0(NT-1,:); Ic2=Ic2(1:NC)';
dMxdt(Ic)= V(NT-1)*y(Ic2) - LT*x(Ic) - D*x(Ic);


% Output
xprime=dMxdt(:);


