function xprime=rectifier(t,X,U) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%  Program written by Stathis Skouras, April 2002  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rectifier - This is a nonlinear model of a rectifier batch distillation column with a condenser,
%    Nr theoretical stages and a reboiler. Reboiler is considered as an additional stage. 
%    The number of real stages is NT=Nr+2. For practical reasons we consider the condenser (top vessel) 
%    as an additional stage. Stages numbered from the top to the bottom. 
%    First 'stage' is condenser, then Nr real stages in the rectifying section and finally the reboiler
%    Model assumptions: 
%    1. Staged distillation column sections
%    2. Perfect mixing and equilibrium in all stage and vessels
%    3. No vapor holdup on the stages
%    4. Constant liquid holdup on all stages (i.e., neglecting liquid dynamics)
%    5. Constant molar vapor flows in all stages (i.e., simplified energy balances) 
%    6. Constant pressure
%    7. Total condenser
%    8. Ideal behaviour in the vapor phase
%    9. VLE modelled by NRTL equation 
%      
%
% Inputs and controller parameters 
%            t    - time in [min].
%            X    - Contains initial states for all the stages.
%                   The first NT molar fractions are for light component A
%                   The next NT molar fractions are for intermediate component B
%                   NT molar liquid holdups in all stages
%                   NT temperatures in all stages
%            U(1) - Boilup, VB
%            U(2) - Nominal liquid flows in the column, L0
% Outputs:   xprime - vector with time derivative of all the states 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Number of stages in the rectifying section: 
Nr=51;
% Total number of 'stages'  
NT=Nr+2;
    
% Splitting the states
xA=X(1:NT)';                            % Liquid composition of A from top to btm
xB=X(NT+1:2*NT)';                       % Liquid composition of B from top to btm
M=X(2*NT+1:3*NT)';                      % Liquid hold up from top to btm
T=X(3*NT+1:4*NT);                       % Temperature from top to btm

% Inputs and disturbances
VB=U(1);                                % Boilup
L0=U(2);                                % Nominal liquid flows in the column (achieved only at steady state)


% THE MODEL

% Vapor-liquid equilibria by NRTL
Y=NRTL(xA,xB,T,NT);                     % Calls 'NRTL' VLE model and returns vapor compositions, temperature and stage pressure

% 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 difference vector (for BP calculations)

% Vapor flows assuming constant molar flows
i=1:NT;    V(i)=VB*ones(1,NT);

% Liquid flows in column sections
% Top section
i=1:NT;   
L(i) = L0;                              % Start-up period under total reflux  
       

% Material balances for the column 1. Overall material balances 2. Component material balances

% Total condenser (not considered as an equilibrium stage)
dMdt = V(2) - L(1);
dMxAdt= V(2).*yA(2) - L(1).*xA(1); 
dMxBdt= V(2).*yB(2) - L(1).*xB(1);

% Rectifying section (not the last stage)
i=2:NT-1;
dMdt(i) = L(i-1)           - L(i)        + V(i+1)          - V(i);                   %=0
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);

% Reboiler (assumed to be an equilibrium stage)
dMdt(NT) = L(NT-1)      - V(NT);
dMxAdt(NT)= L(NT-1)*xA(NT-1) - V(NT)*yA(NT);
dMxBdt(NT)= L(NT-1)*xB(NT-1) - V(NT)*yB(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
xprime=[dxAdt';dxBdt';dMdt';DP'];