function residue=colamodSS(x,U) 
%
%           This is a nonlinear steady state model of a distillation column with
%           NT-1 theoretical stages including a reboiler (stage 1) plus a
%           total condenser ("stage" NT).
%
%           Model assumptions:
%
%           Two components (binary separation); 
%           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). The model has NT states.
%
% Inputs:    x    - states, the NT compositions of light component A with 
%                   reboiler/bottom stage as x(1) and condenser as x(NT). 
%            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.
%            U(8) - number of stages, NT.
%            U(9) - location of feed stage, NF.  
%            U(10)- relative volatility, alpha.
%
% Outputs:   [residue] = f(x), residue=0 if x is a solution of the
%            system of nonlinear equations

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

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

% 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

NT = U(8);                            % Number of stages (including reboiler and total condenser:
NF = U(9);                            % Location of feed stage (stages are counted from the bottom):
alpha = U(10);                        % Relative volatility


%Preallocation
y=ones(NT-1,1);
dMxdt = ones(NT,1);
dMdt = ones(2,1);

% THE MODEL

% Vapor-liquid equilibria
i=1:NT-1;    y(i)=alpha*x(i)./(1+(alpha-1)*x(i));
% --------------------Column--------------------------
%
% Component balances
% =============================================
% Reboiler (assumed to be an equilibrium stage)
dMxdt(1)= (LT+qF*F)*x(2) - VB*y(1) - B*x(1);
% Stripping section trays
i=2:NF-1; 
dMxdt(i)= (LT+qF*F)*x(i+1) - (LT+qF*F)*x(i) + VB*y(i-1) - VB*y(i);
% Feed tray
dMxdt(NF)= LT*x(NF+1) - (LT+qF*F)*x(NF) + VB*y(NF-1) - (VB+(1-qF)*F)*y(NF) + F*zF;
% Enrichment section trays
i=NF+1:NT-1; 
dMxdt(i)= (LT)*x(i+1) - LT*x(i) + (VB+(1-qF)*F)*y(i-1) - (VB+(1-qF)*F)*y(i);
% Total condenser (no equilibrium stage)
dMxdt(NT)= (VB+(1-qF)*F)*y(NT-1) - LT*x(NT) - D*x(NT);
%===============================================
%
% Mass balances
%===============================================
% Reboiler
dMdt(1) = LT+qF*F      - VB      - B;
% Condenser
dMdt(2) = VB+(1-qF)*F  - LT      - D;
% ==============================================
%------------------------------------------------------
%
% Output
residue=[dMxdt; dMdt];

