function [Ceq, dCeqdx] = CSTR(xs,par,scl,con,~,itg)

% Simple CSTR Cascade
% Written by Julian Straus 20150730

% This function defines a simple CSTR cascade to replicate a PFR in
% equation orentated mode.

% The state variables are defined as (first inlet and outlet)
%       xs(1):      Flow [mol/s]
%       xs(2):      Pressure [bar]
%       xs(3):      Temperature [C]
%       xs(4:j+3):  Molar fractions of the j components, e(0,1)    

% The necessary structures are defined as

%       par:        Parameters used in the calculation, must include:
%                   .n      Number of sequential CSTR reactors
%                   .Mr     Total Holdup of the reactor sequence in [m^3]
%                   .stoi   Stoichometry of the components

%                   .for.A  Prexponential factor of the forward reaction
%                   .for.Ea Activation energy of the forward reaction
%                   .for.b  Factor for the temperature dependency

%                   .bac.A  Prexponential factor of the backward reaction
%                   .bac.Ea Activation energy of the backward reaction
%                   .bac.b  Factor for the temperature dependency

%       scl:        Scaling factors of the system, must include:
%                   .Q      Scaling factor for the flow
%                   .p      Scaling factor for the pressure
%                   .T      Scaling factor for the temperature

%       con:        Natural constants, must include:
%                   .R      Ideal gas constant in [J/k/mol]

%       itg:        Intiger values, which should include
%                   .k    	Total number of variables per stream
%                   .j     	Total number of components in each stream

% The output is defined
%   	Ceq:        Equality Constraint for the optimizer/solver
%       dCeqdx:     Jacobian of the equality constraints

% Predefinition of the derivatives
dNdt = sparse(itg.j,par.n);  % Derivative of the mole fraction

% Rescaling of the variables
xs = xs.*scl.Var';

% Redefinition of the state variables
xs = reshape(xs,[itg.k,par.n+1]);

Q = xs(1,:);
p = xs(2,:);
T = xs(3,:)+273.15;
x = xs(4:end,:);

%% Calculation of the reaction terms
% Calculation of the reaction rate coefficient
k = par.for.A*exp(-par.for.Ea./(con.R*T(2:end))).*T(2:end).^par.for.b;
kp = par.bac.A*exp(-par.bac.Ea./(con.R*T(2:end))).*T(2:end).^par.bac.b;

% Calculation of the reaction rates and introduction into the differential
% equations
Mn = par.Mr/par.n;
r = (k.*x(2,2:end).*p(2:end) - kp.*x(3,2:end).^2.*x(1,2:end).^-3./p(2:end))*Mn;
dNdt = dNdt + par.stoi*r;

%% Calculation of mass, energy and pressure balances

Qx = ones(itg.j,1)*Q.*x(:,1:end);

% Calculation of input flows of each reactor
dNdt = dNdt + Qx(:,1:end-1);

% Calculation of output flows of each reactor
dNdt = dNdt - Qx(:,2:end);

% Calculation of the flowrate at the end of each reactor
Qvar = Q(1:par.n) - Q(2:par.n+1) + sum(par.stoi)*r;

% Calculation of the pressure at the end of each reactor 
pvar = par.k^2*par.n*(p(1:par.n).^2-p(2:par.n+1).^2) - Q(2:par.n+1).^2; % p(2:par.n+1)-p(1:par.n); %

% Calculation of the energy balance
dEdt = - (T(2:par.n+1).*Q(2:par.n+1)-T(1:par.n).*Q(1:par.n))*par.cp...
       - par.dH*r - par.cool/par.n;

% Combination of all equality equations
Ceq = [pvar;...
       Qvar;...
       dEdt;...
       dNdt];

Ceq = reshape(Ceq,[],1);

% Scaling of the Equality Constraints
Ceq = Ceq.*scl.Eq;

% Calculation of the gradient of the functions
switch nargout
case 2
% Calculation of the derivatives of the reaction rate
RT2 = con.R*T(2:end).^2;
Tm1 = 1./T(2:end);
drdp = (k.*x(2,2:end) + kp.*x(3,2:end).^2.*x(1,2:end).^-3./p(2:end).^2)*Mn;
drdT = ((par.for.Ea./RT2-Tm1).*k.*x(2,2:end).*p(2:end) -...
       (par.bac.Ea./RT2-Tm1).*kp.*x(3,2:end).^2.*x(1,2:end).^-3./p(2:end)) *...
        Mn;
drdH2 = 3* kp.*x(3,2:end).^2.*x(1,2:end).^-4./p(2:end)*Mn;
drdN2 = k.*p(2:end)*Mn;
drdNH3 = - 2*kp.*x(3,2:end).*x(1,2:end).^-3./p(2:end)*Mn;

% Definition of matrixes necesary for the introduction of the
% derivatives of the reaction rate
Stoi.drdp = par.stoi*drdp;
Stoi.drdT = par.stoi*drdT;
Stoi.drdH2 = par.stoi*drdH2;
Stoi.drdN2 = par.stoi*drdN2;
Stoi.drdNH3 = par.stoi*drdNH3;

%% Jacobian of the mass balances
% Definition of the rows
nk = itg.k*par.n;
nj = itg.j*par.n;
Rx.Nin = (4:itg.k)'*ones(1,par.n)+ones(itg.j,1)*(0:itg.k:nk-itg.k);
Rx.N = Rx.Nin(:)*ones(1,4);

% Definition of the Jacobian of the Mass Balances with respect to the
% flow rate
Cx.Qin = ones(itg.j,1)*(1:itg.k:par.n*itg.k); Cx.Qin = Cx.Qin(:)';
Cx.Qout = Cx.Qin+itg.k;
Vx.Qin = x(1:nj);
Vx.Qout = -x(6:(par.n+1)*itg.j);

% Definition of the Jacobian of the Mass Balances with respect to the
% pressure and the temperature
Cx.p = Cx.Qout + 1;
Cx.T = Cx.Qout + 2;
Vx.p = Stoi.drdp(:)';
Vx.T = Stoi.drdT(:)';

% Definition of the Jacobian of the Mass Balances with respect to the
% inlet concentration
Rx.Rout = Rx.N(:,1:3);
Rx.Rin = [4:itg.k:nk 5:itg.k:nk 6:itg.k:nk];
Rx.I = [7:itg.k:nk 8:itg.k:nk];

Cx.H2in = 4:itg.k:par.n*itg.k;
Cx.N2in = Cx.H2in + 1;
Cx.NH3in = Cx.H2in + 2;
Cx.HI1in = Cx.H2in + 3;
Cx.HI2in = Cx.H2in + 4;

Cx.H2out = Cx.Qout + 3;
Cx.N2out = Cx.Qout + 4;
Cx.NH3out = Cx.Qout + 5;
Cx.HI1out = (itg.k+7:itg.k:(par.n+1)*itg.k);
Cx.HI2out = (itg.k+8:itg.k:(par.n+1)*itg.k);

Vx.H2out = Stoi.drdH2(:)';
Vx.H2out(1:itg.j:end) = Vx.H2out(1:itg.j:end) - Q(2:end);
Vx.N2out = Stoi.drdN2(:)';
Vx.N2out(2:itg.j:end) = Vx.N2out(2:itg.j:end) - Q(2:end);
Vx.NH3out = Stoi.drdNH3(:)';
Vx.NH3out(3:itg.j:end) = Vx.NH3out(3:itg.j:end) - Q(2:end);
Vx.xQin = Q(1:end-1)'*ones(1,5);

% Definition of the matrices as column vectors
Rx.N = Rx.N(:)';
Rx.Rout = Rx.Rout(:)';
Vx.xQin = Vx.xQin(:)';

% Definition of the input column indices and the values
R.x = [Rx.N Rx.Rin Rx.I Rx.Rout Rx.I];
C.x = [Cx.Qin Cx.Qout ...                                   Column of the flow
       Cx.p Cx.T ...                                        Column of the pressure and temperature
       Cx.H2in Cx.N2in Cx.NH3in Cx.HI1in Cx.HI2in ...       Column of the inlet mole fractions  
       Cx.H2out Cx.N2out Cx.NH3out Cx.HI1out Cx.HI2out];... Column of the outlet mole fractions

V.x = [Vx.Qin Vx.Qout ...                                   Derivative w.r.t. the flow
       Vx.p Vx.T ...                                        Derivative w.r.t. the pressure and temperature
       Vx.xQin ...                                          Derivative w.r.t. inlet mole fractions
       Vx.H2out Vx.N2out Vx.NH3out -Q(2:end) -Q(2:end)];... Derivative w.r.t. outlet mole fractions

%% Jacobian of the flow balances
% Definition of the rows
R.n = 1:itg.k:nk;
R.tot = (ones(8,1)*R.n)';
R.tot = R.tot(:)';
R.Q = R.tot(1:7*par.n)+1;

aux1 = 1:itg.j:par.n*itg.j;

% Calculation and scaling of the Jacobian parts
VQ.Q = ones(1,par.n);
VQ.p = sum(Stoi.drdp);
VQ.T = sum(Stoi.drdT);
VQ.H2 = sum(Stoi.drdH2);
VQ.N2 = sum(Stoi.drdN2);
VQ.NH3 = sum(Stoi.drdNH3);

% Stacking of the vectors
R.Q = R.Q;
C.Q = [Cx.Qin(aux1) Cx.Qout(aux1) ...                       Column of the flow
       Cx.p(aux1) Cx.T(aux1) ...                          	Column of the pressure and temperature
       Cx.H2out(aux1) Cx.N2out(aux1) Cx.NH3out(aux1)];...	Column of the mole fractions

V.Q = [VQ.Q -VQ.Q ...                                       Derivative w.r.t. the flow
       VQ.p VQ.T ...                                        Derivative w.r.t. the pressure and temperature
       VQ.H2 VQ.N2 VQ.NH3];...                              Derivative w.r.t. the mole fractions
  
%% Gradient of the pressure balances
% Definition of the rows
R.p = R.tot(1:3*par.n);

% Calculation and scaling of the Jacobian parts
Vp.Q = -2*Q(2:end);
Vp.pin = 2*par.k^2*par.n*p(1:end-1);
Vp.pout = -2*par.k^2*par.n*p(2:end);

% Stacking of the vectors
C.p = [Cx.Qout(aux1) ...                                    Column of the flow
       Cx.p(aux1)-itg.k Cx.p(aux1)];...                     Column of the pressure

V.p = [Vp.Q ...                                             Derivative w.r.t. the flow
       Vp.pin Vp.pout];...                                  Derivative w.r.t. the pressure
  
%% Gradient of the temperature balances
% Definition of the rows
R.E = R.tot+2;

% Calculation and scaling of the Jacobian parts
VE.Qin = T(1:end-1)*par.cp;
VE.Qout = -T(2:end)*par.cp;
VE.p = -drdp*par.dH;
VE.Tin = Q(1:end-1)*par.cp;
VE.Tout = (-Q(2:end)*par.cp-par.dH*drdT);
VE.H2 = -par.dH*drdH2;
VE.N2 = -par.dH*drdN2;
VE.NH3 = -par.dH*drdNH3;

% Stacking of the vectors
C.E = [Cx.Qin(aux1) Cx.Qout(aux1) ...                       Column of the flow
       Cx.p(aux1) ...                                       Column of the pressure
       Cx.T(aux1)-itg.k Cx.T(aux1) ...                      Column of the temperature
       Cx.H2out(aux1) Cx.N2out(aux1) Cx.NH3out(aux1)];...	Column of the mole fractions

V.E = [VE.Qin VE.Qout ...                                   Derivative w.r.t. the flow
       VE.p ...                                             Derivative w.r.t. the pressure
       VE.Tin VE.Tout ...                                   Derivative w.r.t. the temperature
       VE.H2 VE.N2 VE.NH3];...                              Derivative w.r.t. the mole fractions

   
% Definition of the combined vectors
R.tot = [R.Q R.p R.E R.x];
C.tot = [C.Q C.p C.E C.x];
V.tot = [V.Q V.p V.E V.x]...
        .*scl.Eq(R.tot)'.*scl.Var(C.tot);
   
% Definition of the sparse matrix
dCeqdx = sparse(R.tot,C.tot,V.tot,length(Ceq),numel(xs));

end