function [xprime B D LT VB xB xD e] = multicomA(t,varargin) %#ok<STOUT>

% ***************** Program written by Stathis Skouras in March, 2001 *****************
% ***************** This version: March 21, 2008 by Antonio Araujo ********************
% 
% 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 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

parameters
global config number
%Now application of this file is insensitive to whether we choose to use
%fsolve or ode15s to solve it
switch nargin
    case 3 % Steady state
        X  = t ;
        U  = varargin{1} ;
        Dd = varargin{2} ;
    case 4 % Dynamic
        X  = varargin{1};
        U  = varargin{2};
        Dd = varargin{3};
end


% 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

% Rearrange elements of composition vector (x) for later use
Iu = reshape(1:(NC-1)*NT,NT,NC-1);
x  = x(Iu);                          

xB = [x(1,:) 1-sum(x(1,:))] ;
xD = [x(end,:) 1-sum(x(end,:))] ;

% Inputs
switch config
    case('LV')
        LT = U(1);                      % Reflux
        VB = U(2);                      % Boilup
    case('BV')
        B  = U(1);                      % Bottoms flowrate
        VB = U(2);                      % Boilup
    case('XdXb')
        xDs = U(1) ;                    % Set point for top composition
        xBs = U(2) ;                    % Set point for bottom composition
        eint = X(end-1:end)    ;
        if number == 1 % First column
            e(1) = xD(2) - xDs ; %Control B at top of column 1
            e(2) = xB(1) - xBs ; %Control A at bottom of column 1
        elseif number == 2 
            e(1) = xD(3) - xDs ; %Control C at top of column 2
            e(2) = xB(2) - xBs ; %Control B at bottom of column 2  
        end
end

% Disturbances
F  = Dd(1);                             % Feedrate
qF = Dd(2);                             % Feed liquid fraction
zF = Dd(3:end);                         % Feed compositions

% The distillation model
flag = 0; % This flag signals the inclusion of equations for liquid levels. If flag~=0, Simulink is being used and no liquid level is calculated in the code.
model


% Output
switch config
    case('LV')
        e = [];
        xprime = [dxdt; dMdt];
    case('BV') 
        e = [] ;
        xprime = [dxdt; dMdt];
    case('XdXb')
        e = e(:) ;
        xprime = [dxdt; dMdt ; e] ; % Include error 
end


