function [c,ceq,DC,DCeq]=nonlcon_colm(ux,p)
%
% nonlcon_colm - 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: Optimization %Dynamic simulation
%
%           NB: Manually change necessary when change no. components
%
%           The model is based on column A in Skogestad and Postlethwaite
%           (1996). The model has NT*NC states.
%
% Inputs:    U=UX(1:NC+6); X=UX(NC+7:NC*NT+NC+6)
%            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.
%            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.
%
%            P    - parameter vector - explained below 
%
% Outputs:   c    - inequality constraints 
%            ceq  - equality constraints
%            dc   - gradient to inequality constraints
%            dceq - gradient to equality constraints
%
% Modified by: 
%         msg20000106 Multicomponent column (NC)
%                          Extern parameters 
%         msg20010109 dydx=dydx*dxdx
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------------------------------------------
% Parameters
NC=p(1); % Number of components
NT=p(2); % Number of stages (including reboiler and total condenser: 
NF=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
%------------------------------------------------------------

if 1==0
M0=[1.75157700000000
   0.24999996319485
   0.24999999592758
   0.25000000488171
   0.25000000180141
   0.25000000459198
   0.24999998034159
   0.24999994036494
   0.25000005371721
   0.25000007294201
   0.24999973230857
   0.25000015004953
   0.25000065349722
   0.24999951624868
   0.24999984101420
   0.24999999004301
   0.25000009529941
   0.24999987648075
   0.25000006610026
   0.25000003254892
   0.24999999988138
   0.25000001836307
   0.24999999750362
   0.25000002523682
   0.25000003638026
   0.25000013717866
   0.25000240487654
   0.24999756568588
   0.25000314494121
   0.25000314494121
   0.64999990000000]*1e2;
M0=M0(:);
end
%M0(1)=275;
%M0(NT)=185;

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

U=ux(1:NC+6);
U=U(:);
X=ux(NC+7:NC*NT+NC+6);
X=X(:);

%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:NC+5); zF=zF(:);             % Feed composition
qF = U(NC+6);                         % 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);

%NB! No. of eye = no. of comp.
dMdx=[zeros(NT,NC+6) eye(NT) eye(NT)];        % NC=2
%dMdx=[zeros(NT,NC+6) eye(NT) eye(NT) eye(NT)]; % NC=3
dMmdx=dMdx(Im,:);

%Liquid fraction
x=Mi./Mm;
dxdx=[zeros(NC*NT,NC+6) diag(1./Mm)]-diag(Mi./(Mm.*Mm))*dMmdx;

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

%NB! No. of eye = no. of components
dSdx=[alpha(1)*eye(NT) alpha(2)*eye(NT)];   %NC=2 
%dSdx=[alpha(1)*eye(NT) alpha(2)*eye(NT) alpha(3)*eye(NT)]; %NC=3

dSmdx=dSdx(Im,:);
%dSmdx(:,9:20),size(dSmdx(:,9:20)),pause

dydx=[diag(alpham./sumalphaxm)];
%dydx(:,9:20),size(dydx(:,9:20)),pause

dydSm=-diag(alphaxm./(sumalphaxm.*sumalphaxm));
%dydSm,size(dydSm),pause

dydx=dydSm*dSmdx+dydx;
dydx=dydx*dxdx; 
%dydx(:,9:20),size(dydx(:,9:20)),pause


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

dVdx=zeros(NT,NT*NC+NC+6);
dVdx(:,2)=1;
dVdx(NF:NT-1,5)=1-qF;
dVdx(NF:NT-1,NC+6)=-F;
dVmdx=dVdx(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);

dLdx=zeros(NT,NT*NC+NC+6);
dLdx(2:NT-1,:)=1/taul*dMdx(2:NT-1,:)+lambda*dVdx(1:NT-2,:);
dLdx(NT,1)=1;
dLmdx=dLdx(Im,:);

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


% 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);
%dMxdt(Ic)=-Vm(Ic).*y(Ic);

IC=Ic*ones(1,NT*NC+NC+6);ICP=Icp*ones(1,NT*NC+NC+6);ICM=Icm*ones(1,NT*NC+NC+6);
dfdx(Ic,:)=x(ICP).*dLmdx(Icp,:)+Lm(ICP).*dxdx(Icp,:) ...
          -x(IC) .*dLmdx(Ic,:) -Lm(IC) .*dxdx(Ic,:) ...
	  +y(ICM).*dVmdx(Icm,:)+Vm(ICM).*dydx(Icm,:) ...
	  -y(IC) .*dVmdx(Ic,:) -Vm(IC) .*dydx(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)';
dMxdt(Ic)= dMxdt(Ic) + F*zF;

dfdx(Ic,5)     =dfdx(Ic,5)     +zF;
dfdx(Ic,6:NC+5)=dfdx(Ic,6:NC+5)+F*eye(NC);


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

IC=Ic*ones(1,NT*NC+NC+6);IC2=Ic2*ones(1,NT*NC+NC+6);
dfdx(Ic,:)=dfdx(Ic,:)+x(IC2).*dLmdx(Ic2,:)+Lm(IC2).*dxdx(Ic2,:) ...
                     -y(IC) .*dVmdx(Ic,:) -Vm(IC) .*dydx(Ic,:)  ...
		     -B      *dxdx(Ic,:);
dfdx(Ic,4)=dfdx(Ic,4)-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);
%dMxdt(Ic),pause

IC=Ic*ones(1,NT*NC+NC+6);IC2=Ic2*ones(1,NT*NC+NC+6);
dfdx(Ic,:)=dfdx(Ic,:)+y(IC2).*dVmdx(Ic2,:)+Vm(IC2).*dydx(Ic2,:) ...
                     -x(IC) .*dLmdx(Ic,:) -Lm(IC) .*dxdx(Ic,:)  ...
		     -D*dxdx(Ic,:);
dfdx(Ic,3)=dfdx(Ic,3)-x(Ic);


% Output
c=[];
ceq=dMxdt;
DC=[];
DCeq=dfdx';






