function [Ceq, dCeqdx] = ReacNN(x,par,scl,~,~,itg)

% Neural network model of the reaction section
% Written by Julian Straus 20150805

% This function defines the neural network model for the reaction section
% of the simple ammonia loop.

% The state variables are defined as
% (flow, pressure, temperature, mole fractions)
%       x(1:k1):        Dummy variables
%       x(k1+1:k2):     Connection to the separation section
%       x(k2+1:k3):     Connection to the Makeup section
%       x(k3+1:k4):     Dummy variables
%       x(k4+1:k4+3):   Parameters, Splitratios and outlet temperature of HEx5

% The necessary structures are defined as
%       par:        Parameters used in the calculation, must include:
%                   .xload      Cell structure containing multiplier, 
%                               if latent variables are used
%                   .net        Cell structure containing the neural
%                               network
%                   .dep        Logic value for dependency
%                   .diff       Logic value for differences
%                   .stoi       Stoichometry of the system

%       scl:        Scaling factors of the system, must include:
%                   .Var        Scaling factor for the variables
%                   .Eq         Scaling factor for the equations
%                   .x          Structure for the x values (mu and std)
%                   .U          Structure for the U values (mu and std)


%       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

% Definition of multiples of itg.k
k1 = itg.k;
k2 = 2*itg.k;
k3 = 3*itg.k;
k4 = 4*itg.k;

% Scaling and Reassignment of the variables
x = x.*scl.Var';
Qin  = x(k2+1);
pin  = x(k2+2);
Tin  = x(k2+3);
xin  = x(k2+4:k3);

Qout = x(k1+1);
pout = x(k1+2);
Tout = x(k1+3);
xout = x(k1+4:k2);

ratio = x(k4+1:k4+2);
TRef  = x(k4+3);

% Definition of the reduced matrix size for the neural networks
if par.dep
    U = ones(11,1)*[Qin pin Tin xin(1)/xin(2) xin(3:5)' ratio' TRef];
else
    U = ones(11,1)*[Qin pin Tin xin(2:5)' ratio' Tref];
end
dU = 1e-6*diag(diag(U(2:end,:)));
U(2:end,:) = U(2:end,:)+dU;
U = (U-ones(size(U,1),1)*scl.U.mu)./(ones(size(U,1),1)*scl.U.std);
Cm{1} = U*par.xload{1};
Cm{2} = U*par.xload{2};
Cm{3} = U*par.xload{3};

% Calculation of the outlet pressure, the outlet temperature, and the
% extent of reaction
p_NN   = par.net{1}(Cm{1}')*scl.x.std(1) + scl.x.mu(1);
T_NN   = par.net{2}(Cm{2}')*scl.x.std(2) + scl.x.mu(2);
EoR_NN = par.net{3}(Cm{3}')*scl.x.std(3) + scl.x.mu(3);

% Definition of the total mass balance
Ceq = zeros(8,1);
Ceq(1) = Qin - Qout + sum(par.stoi)*EoR_NN(1);

% Definition of the energy and temperature balance
if par.diff
    Ceq(2) = pin - p_NN(1) - pout;
    Ceq(3) = Tin - T_NN(1) - Tout;
else
    Ceq(2) = pout_NN(1) - pout;
    Ceq(3) = Tout_NN(1) - Tout;
end

% Definition of the component mass balances
Ceq(4:8) = Qin*xin - Qout*xout + par.stoi*EoR_NN(1);

% 

% Scaling of the equality constraints
Ceq = Ceq.*scl.Eq;

% Calculation of the gradient of the functions
switch nargout
    case 2
        % Definition of auxilliary values
        Onj1 = ones(itg.j,1);

        % Definition of the Jacobian for the system
        dpoutdx = (p_NN(2:end)-p_NN(1))./diag(dU)';
        dToutdx = (T_NN(2:end)-T_NN(1))./diag(dU)';
        dEoRdx  = (EoR_NN(2:end)-EoR_NN(1))./diag(dU)';

        % Definition of the column values
        if par.dep
            C.Q = [k1+1 k2+1:k3 k4+1:k4+3];...              Ceq(1) Qout, stream input, and parameters
            C.p = [k1+2 k2+1:k3 k4+1:k4+3];...              Ceq(2) pout, stream input,and parameters
            C.T = [k1+3 k2+1:k3 k4+1:k4+3];...              Ceq(3) Tout, stream input, Tout, and parameters
            dpoutdx = [dpoutdx(1:3) dpoutdx(4)/xin(2) -dpoutdx(4)*xin(1)/xin(2)^2 dpoutdx(5:end)]; 
            dToutdx = [dToutdx(1:3) dToutdx(4)/xin(2) -dToutdx(4)*xin(1)/xin(2)^2 dToutdx(5:end)]; 
            dEoRdx  = [dEoRdx(1:3)  dEoRdx(4)/xin(2)  -dEoRdx(4)*xin(1)/xin(2)^2  dEoRdx(5:end)];
            str = 8;
        else
            C.Q = [k1+1 k2+1:k2+3 k2+5:k3 k4+1:k4+3];... 	Ceq(1) Qout, stream input, and parameters
            C.p = [k1+2 k2+1:k2+3 k2+5:k3 k4+1:k4+3];...    Ceq(2) pout, stream input,and parameters
            C.T = [k1+3 k2+1:k2+3 k2+5:k3 k4+1:k4+3];...    Ceq(3) Tout, stream input, Tout, and parameters
            str = 7;
        end
        C.x = [Onj1*([k2+1:k2+3 k4+1:k4+3]) ...          Ceq(4:8) Qin, pin, Tin, and parameters
               Onj1*([k2+4:k3 k1+1]) ...             Ceq(4:8) xin and Qout
               (k1+4:k2)'];...                      Ceq(4:8) xout
        C.x = C.x(:)';
        val = str+1;

        % Definition of the row values
        R.Q = ones(1,length(C.Q));...     	Ceq(1)
        R.p = ones(1,length(C.p))*2;...     Ceq(2)
        R.T = ones(1,length(C.T))*3;...     Ceq(3)
        R.x = (4:8)'*ones(1,itg.k+5);...    Ceq(4:8)
        R.x = R.x(:)';

        % Definition of the values of the derivatives
        V.Q = [-1 eye(1,str)+sum(par.stoi(1:5))*dEoRdx(1:str) ...   Derivatives of Ceq(1) w.r.t to Qout and the stream input
                sum(par.stoi)*dEoRdx(val:end)];...                  Derivatives of Ceq(1) w.r.t to the parameters
        V.p = [-1 dpoutdx(1)   dpoutdx(2) dpoutdx(3:str) ...        Derivatives of Ceq(2) w.r.t to pout and the stream input 
               -dpoutdx(val:end)];...                               Derivatives of Ceq(2) w.r.t to the parameters
        V.T = [-1 dToutdx(1:2) dToutdx(3) dToutdx(4:str) ...        Derivatives of Ceq(3) w.r.t to Tout and the stream input
               dToutdx(val:end)];...                                Derivatives of Ceq(3) w.r.t to the parameters
        if par.diff
            V.p(2:end) = -V.p(2:end);
            V.p(3) = 1+V.p(3);
            V.T(2:end) = -V.T(2:end);
            V.T(4) = 1+V.T(4);
        end
        V.x = [xin+par.stoi*dEoRdx(1) ...                           Derivatives of Ceq(4:8) w.r.t. to Qin
               par.stoi*dEoRdx([2:3 val:end]) ...                   Derivatives of Ceq(4:8) w.r.t. to pin, Tin and parameters
               Qin*eye(itg.j)+ ...
               par.stoi*[zeros(8-str) dEoRdx(4:str)] ...            Derivatives of Ceq(4:8) w.r.t. to xin
              -xout ...                                             Derivatives of Ceq(4:8) w.r.t. to Qout
              -Onj1*Qout];...                                       Derivatives of Ceq(4:8) w.r.t. to xout
        V.x = V.x(:)';

        % Definition of the combined vectors
        R.tot = [R.Q R.p R.T R.x];
        C.tot = [C.Q C.p C.T C.x];
        V.tot = [V.Q V.p V.T V.x]...
                 .*scl.Eq(R.tot)'.*scl.Var(C.tot);

        % Definition of the sparse Jacobian matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k1,k4+3);
end