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

% Neural network model of the separatio\n section
% Written by Julian Straus 20150805

% This function defines a simple isenthalpic, isobaric split.

% The state variables are defined as
% (flow, pressure, temperature, mole fractions)
%       x(1:k1):        Connection to the reaction section
%       x(k1+1:k2):     Liquid outlet
%       x(k2+1:k3):     Heat Exchanger inlet
%       x(k3+1:k4):     Gas outlet
%       x(k4+1:k5):     Heat Exchanger outlet

% 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

%       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;
k5 = 5*itg.k;

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

QLout = x(k1+1);
pLout = x(k1+2);
TLout = x(k1+3);
xLout = x(k1+4:k2);

QHin = x(k2+1);
pHin = x(k2+2);
THin = x(k2+3);
xHin = x(k2+4:k3);

QGout = x(k3+1);
pGout = x(k3+2);
TGout = x(k3+3);
xGout = x(k3+4:k4);

QHout = x(k4+1);
pHout = x(k4+2);
THout = x(k4+3);
xHout = x(k4+4:k5);

ratio = x(k5+1);
TRef  = x(k5+2:k5+4);

Qtrans = x(k5+5:k5+7);

% Definition of the reduced matrix size for the neural networks
if par.dep
    U = ones(15,1)*[Qin pin Tin xin(1)/xin(2) xin(3:5)' QHin pHin THin ratio TRef'];
else
    U = ones(15,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);

% Calculation of the latent variables
Cm{1} = U*par.xload{1};
Cm{2} = U*par.xload{2};
Cm{3} = U*par.xload{3};
Cm{4} = U*par.xload{4};
Cm{5} = U*par.xload{5};
Cm{6} = U*par.xload{6};
Cm{7} = U*par.xload{7};
Cm{8} = U*par.xload{8};
Cm{9} = U*par.xload{9};
Cm{10} = U*par.xload{10};
Cm{11} = U*par.xload{11};
Cm{12} = U*par.xload{12};
Cm{13} = U*par.xload{13};

% Calculation of the outlet values
p_G = par.net{1}(Cm{1}')*scl.x.std(1) + scl.x.mu(1);
T_G = par.net{2}(Cm{2}')*scl.x.std(2) + scl.x.mu(2);
p_L = par.net{3}(Cm{3}')*scl.x.std(3) + scl.x.mu(3);
T_L = par.net{4}(Cm{4}')*scl.x.std(4) + scl.x.mu(4);

T_H = par.net{5}(Cm{5}')*scl.x.std(5) + scl.x.mu(5);

alpha = [par.net{6}(Cm{6}')' par.net{7}(Cm{7}')' ...
         par.net{8}(Cm{8}')' par.net{9}(Cm{9}')' ...
         par.net{10}(Cm{10}')'].*(ones(15,1)*scl.x.std(6:10)) + ...
                                 (ones(15,1)*scl.x.mu(6:10));

Q_HEx8  = par.net{11}(Cm{11}')*scl.x.std(11) + scl.x.mu(11);
Q_HEx9  = par.net{12}(Cm{12}')*scl.x.std(12) + scl.x.mu(12);
Q_HEx10 = par.net{13}(Cm{13}')*scl.x.std(13) + scl.x.mu(13);

% Definition of the total mass balance
Ceq = zeros(27,1);
Ceq(1) = (Qin*(1-alpha(1,:))*xin) - QLout;
Ceq(9) = (Qin*alpha(1,:)*xin) - QGout;

% Definition of the pressure and temperature balance
if par.diff
    Ceq(2)  = pin - p_L(1) - pLout;
    Ceq(3)  = Tin - T_L(1) - TLout;
    Ceq(10) = pin - p_G(1) - pGout;
    Ceq(11) = Tin - T_G(1) - TGout;
    Ceq(18) = pHin - pHout;
    Ceq(19) = THin - T_H(1) - THout;
else
    Ceq(2)  = p_L(1) - pLout;
    Ceq(3)  = T_L(1) - TLout;
    Ceq(10) = p_G(1) - pGout;
    Ceq(11) = T_G(1) - TGout;
    Ceq(18) = pHin - pHout;
    Ceq(19) = T_H(1) - THout;
end

% Definition of the component mass balances
Ceq(4:8)   = Qin*xin.*(1-alpha(1,:)') - QLout*xLout;
Ceq(12:16) = Qin*xin.*alpha(1,:)' - QGout*xGout;

% Definition of the heat exchanger balances
Ceq(17)    = QHin - QHout;
Ceq(20:24) = QHin*xHin - QHout*xHout;

% Definition of the heat exchanger energy duties
Ceq(25) = Q_HEx8(1)  - Qtrans(1);
Ceq(26) = Q_HEx9(1)  - Qtrans(2);
Ceq(27) = Q_HEx10(1) - Qtrans(3);

% 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);
        On1j = ones(1,itg.j);
        On12 = ones(1,2);
        On14 = ones(1,4);
        On17 = ones(1,7);
        On117 = ones(1,17);

        % Definition of the Jacobian for the system
        dpLoutdx = (p_L(2:end)-p_L(1))./diag(dU)';
        dTLoutdx = (T_L(2:end)-T_L(1))./diag(dU)';
        dpGoutdx = (p_G(2:end)-p_G(1))./diag(dU)';
        dTGoutdx = (T_G(2:end)-T_G(1))./diag(dU)';
        dpHoutdx = 1;
%         dpHoutdx = (p_H(2:end)-p_H(1))./diag(dU)';
        dTHoutdx = (T_H(2:end)-T_H(1))./diag(dU)';
        
        dalphdx  = (alpha(2:end,:)-ones(14,1)*alpha(1,:))./(diag(dU)*On1j);
        
        dQ_HEx8dx  = (Q_HEx8(2:end)-Q_HEx8(1))./diag(dU)';
        dQ_HEx9dx  = (Q_HEx9(2:end)-Q_HEx9(1))./diag(dU)';
        dQ_HEx10dx = (Q_HEx10(2:end)-Q_HEx10(1))./diag(dU)';
        
        % Definition of the column values
        if par.dep
            dum1 = [1:k1 k2+1:k2+3 k5+1:k5+4];...         stream input and parameters
            
            % Rearrangement of the differentials
            dpLoutdx = [dpLoutdx(1:3) dpLoutdx(4)/xin(2) -dpLoutdx(4)*xin(1)/xin(2)^2 dpLoutdx(5:end)]; 
            dTLoutdx = [dTLoutdx(1:3) dTLoutdx(4)/xin(2) -dTLoutdx(4)*xin(1)/xin(2)^2 dTLoutdx(5:end)]; 
            dpGoutdx = [dpGoutdx(1:3) dpGoutdx(4)/xin(2) -dpGoutdx(4)*xin(1)/xin(2)^2 dpGoutdx(5:end)]; 
            dTGoutdx = [dTGoutdx(1:3) dTGoutdx(4)/xin(2) -dTGoutdx(4)*xin(1)/xin(2)^2 dTGoutdx(5:end)]; 
%             dpHoutdx = [dpHoutdx(1:3) dpHoutdx(4)/xin(2) -dpHoutdx(4)*xin(1)/xin(2)^2 dpHoutdx(5:end)]; 
            dTHoutdx = [dTHoutdx(1:3) dTHoutdx(4)/xin(2) -dTHoutdx(4)*xin(1)/xin(2)^2 dTHoutdx(5:end)]; 
            
            dalphdx  = [dalphdx(1:3,:); dalphdx(4,:)/xin(2); -dalphdx(4,:)*xin(1)/xin(2)^2; dalphdx(5:end,:)]'; 
            
            dQ_HEx8dx = [dQ_HEx8dx(1:3) dQ_HEx8dx(4)/xin(2) -dQ_HEx8dx(4)*xin(1)/xin(2)^2 dQ_HEx8dx(5:end)]; 
            dQ_HEx9dx = [dQ_HEx9dx(1:3) dQ_HEx9dx(4)/xin(2) -dQ_HEx9dx(4)*xin(1)/xin(2)^2 dQ_HEx9dx(5:end)]; 
            dQ_HEx10dx = [dQ_HEx10dx(1:3) dQ_HEx10dx(4)/xin(2) -dQ_HEx10dx(4)*xin(1)/xin(2)^2 dQ_HEx10dx(5:end)]; 
            str = 8;
            
        else
            dum1 = ones(3,1)*[1:2 4:k1 k2+1:k2+3 k5+1:k5+4];...     stream input and parameters
           	dalphdx = dalphdx';
            str = 7;
        end
        
        val = str+1;
        
        %% Definition of the column values
        C.Q_L = [dum1 k1+1];...                                 Ceq(1) 
        C.p_L = [dum1 k1+2];...                                 Ceq(2) 
        C.T_L = [dum1 k1+3];...                                 Ceq(3) 
        C.x_L = [Onj1*[dum1 k1+1] (k1+4:k2)'];...               Ceq(4:8)
        C.Q_G = [dum1 k3+1];...                                 Ceq(9) 
        C.p_G = [dum1 k3+2];...                                 Ceq(10) 
        C.T_G = [dum1 k3+3];...                                 Ceq(11) 
        C.x_G = [Onj1*[dum1 k3+1] (k3+4:k4)'];...               Ceq(12:16)
        C.Q_H = [k2+1 k4+1];...                                 Ceq(17)
        C.p_H = [k2+2 k4+2];...                                 Ceq(18) 
        C.T_H = [dum1 k4+3];...                                 Ceq(19) 
        C.x_H = [Onj1*[k2+1 k4+1] (k2+4:k3)' (k4+4:k5)'];...    Ceq(20:24)
        C.Q_8  = [dum1 k5+5];...                                Ceq(25)
        C.Q_9  = [dum1 k5+6];...                                Ceq(26)
        C.Q_10 = [dum1 k5+7];...                                Ceq(27)
        
        % Redefinition of the column values
        C.Q = [C.Q_L C.Q_G C.Q_H];
        C.p = [C.p_L C.p_G C.p_H];
        C.T = [C.T_L C.T_G C.T_H];
        C.x = [C.x_L(:)' C.x_G(:)' C.x_H(:)']; 
        C.Qh = [C.Q_8 C.Q_9 C.Q_10];

        %% Definition of the row values
        R.Q_L = ones(1,length(C.Q_L));...               Ceq(1)
        R.p_L = ones(1,length(C.p_L))*2;...             Ceq(2)
        R.T_L = ones(1,length(C.T_L))*3;...             Ceq(3)
        R.x_L = (4:8)'*On117;...                        Ceq(4:8)
        R.Q_G = ones(1,length(C.Q_G))*(k1+1);...        Ceq(9)
        R.p_G = ones(1,length(C.p_G))*(k1+2);...        Ceq(10)
        R.T_G = ones(1,length(C.T_G))*(k1+3);...        Ceq(11)
        R.x_G = (12:16)'*On117;...                      Ceq(12:16)
        R.Q_H = ones(1,length(C.Q_H))*(k2+1);...        Ceq(17)
        R.p_H = ones(1,length(C.p_H))*(k2+2);...        Ceq(18)
        R.T_H = ones(1,length(C.T_H))*(k2+3);...        Ceq(19)
        R.x_H = (20:24)'*On14;...                       Ceq(20:24)
        
        R.Q_8  = ones(1,length(C.Q_8))*(k3+1);...        Ceq(25)
        R.Q_9  = ones(1,length(C.Q_9))*(k3+2);...        Ceq(26)
        R.Q_10 = ones(1,length(C.Q_10))*(k3+3);...       Ceq(27)
        
        % Redefinition of the row values
        R.Q = [R.Q_L R.Q_G R.Q_H];
        R.p = [R.p_L R.p_G R.p_H];
        R.T = [R.T_L R.T_G R.T_H];
        R.x = [R.x_L(:)' R.x_G(:)' R.x_H(:)']; 
        R.Qh = [R.Q_8 R.Q_9 R.Q_10];
    
        %% Definition of the values of the derivatives
        
        % Definition of the values for the flow balances
        V.Q_L = [(1-alpha(1,:))*xin - Qin*xin'*dalphdx(:,1) ...
                 -Qin*xin'*dalphdx(:,2:3) ... 
                 Qin*(1-alpha(1,:)) - Qin*xin'*[zeros(5,8-str) dalphdx(:,4:str)] ...
                 -Qin*xin'*dalphdx(:,val:end) ...
                 -1];
        V.Q_G = [alpha(1,:)*xin + Qin*xin'*dalphdx(:,1) ...
                 Qin*xin'*dalphdx(:,2:3) ... 
                 Qin*alpha(1,:) + Qin*xin'*[zeros(5,8-str) dalphdx(:,4:str)] ...
                 Qin*xin'*dalphdx(:,val:end) ...
                 -1];
        V.Q_H = [1 -1];
        
        V.p_L = [dpLoutdx -1];
        V.p_G = [dpGoutdx -1];
        V.p_H = [dpHoutdx -1];
        
        V.T_L = [dTLoutdx -1];
        V.T_G = [dTGoutdx -1];
        V.T_H = [dTHoutdx -1];
        
        V.x_L = [xin.*(1-alpha(1,:)') - Qin*xin.*dalphdx(:,1) ...
                 -Qin*(xin*On12).*dalphdx(:,2:3) ...
                 -[zeros(5,8-str) dalphdx(:,4:str)]*Qin.*(xin*On1j) + ...
                 diag((1-alpha(1,:))*Qin) ...
                 -Qin*(xin*On17).*dalphdx(:,val:end) ...
                 -xLout -QLout*Onj1];
        V.x_G = [xin.*alpha(1,:)'+Qin*xin.*dalphdx(:,1) ...
                 Qin*(xin*On12).*dalphdx(:,2:3) ...
                 [zeros(5,8-str) dalphdx(:,4:str)]*Qin.*(xin*On1j) + ...
                 diag(alpha(1,:)*Qin) ...
                 Qin*(xin*On17).*dalphdx(:,val:end) ...
                 -xGout -QGout*Onj1];
        V.x_H = [xHin -xHout ...
                 QHin*Onj1 -QHout*Onj1];
             
        V.Q_8  = [dQ_HEx8dx -1];
        V.Q_9  = [dQ_HEx9dx -1];
        V.Q_10 = [dQ_HEx10dx -1];
        
        % Redefinition of the row values
        V.Q = [V.Q_L V.Q_G V.Q_H];
        V.p = [V.p_L V.p_G V.p_H];
        V.T = [V.T_L V.T_G V.T_H];
        V.x = [V.x_L(:)' V.x_G(:)' V.x_H(:)']; 
        V.Qh = [V.Q_8 V.Q_9 V.Q_10];
        
        if par.diff
            V.p([1:15 17:31]) = -V.p([1:15 17:31]);
            V.p([2 18]) = 1+V.p([2 18]);
            V.T([1:15 17:31 33:47]) = -V.T([1:15 17:31 33:47]);
            V.T([3 19 43]) = 1+V.T([3 19 43]);
        end
        
        %% Definition of the combined vectors
        C.tot = [C.Q C.p C.T C.x C.Qh];
        R.tot = [R.Q R.p R.T R.x R.Qh];
        V.tot = [V.Q V.p V.T V.x V.Qh]...
                 .*scl.Eq(R.tot)'.*scl.Var(C.tot);

        % Definition of the sparse Jacobian matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,k3+3,k5+7);

% dCeqdx  = sparse(24,44);
end