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

% Simple Mixer
% Written by Julian Straus 20160317

% This function defines a simple isenthalpic, isobaric mixer. It assume
% constant heat capcity of all streams at the inlet. If necessary, it can
% also be modified to include a real enthalpy balance with varying cp/cv value.

% The state variables are defined as (first all inlets and then the outlet)
%       x(1):       Flow of the system [mol/s]
%       x(2):       Pressure of the system [bar]
%       x(3):       Temperature of the system [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 inlets

%       scl:        Scaling factors of the system

%       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

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

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

% Switch between concentration calculations or not
Ceq = zeros(itg.j+par.n+2,1);        
        
% Mass balances
Ceq(1) = sum(x(1,1:end-1)) - x(1,end);              % Total flow

% Pressure Neutrality
for i = 1:par.n
    Ceq(1+i) = x(2,i) - x(2,i+1);
end

% Energy balance
Ceq(par.n+2) = x(1,1:par.n)*x(3,1:par.n)' - x(1,end)*x(3,end);

% Mass balances
Ceq(par.n+3:par.n+itg.j+1) = x(4:end-1,1:par.n)*x(1,1:par.n)' - ... % Inlets components
                             x(4:end-1,end)*x(1,end);               % Outlet components

% Stoichometric balance for the outlet stream
Ceq(end) = 1 - sum(x(4:end,end));
Ceq = Ceq.*scl.Eq;

%% Jacobian of the functions
switch nargout
    case 2
        % Definition of auxiliary values
        M1 = linspace(0,par.n,par.n+1)*itg.k+1;                
        On1N = ones(1,par.n+1);
        aux = itg.k*par.n;

        % Definition of the indices for the rows (R) and columns (C)
        R.Q = On1N;                                     % Ceq(1)
        R.p = [1 1]'*(2:par.n+1);                       % Ceq(2:par.n+1)
        R.E = [On1N On1N]*(par.n+2);                    % Ceq(par.n+2)
        R.x = (par.n+3:par.n+itg.j+1)'*ones(1,par.n+1); % Ceq(par.n+3:par.n+itg.j+1)
        R.Sx = ones(5,1)*(par.n+itg.j+2);               % Ceq(end = par.n+itg.j+2)
        R.p = R.p(:)'; R.x = R.x(:)'; R.Sx = R.Sx(:)'; 

        C.Q = M1;                                       % Ceq(1)
        C.p = ones(2,1)*(2:itg.k:(par.n+1)*itg.k);      % Ceq(2:par.n+1)
        C.p = C.p(2:end-1); 
        C.E = [M1 M1+2];                                % Ceq(par.n+2)
        C.xQ = ones(itg.j-1,1)*M1;                      % Ceq(par.n+3:par.n+itg.j+1) w.r.t. the flowrates
        C.xx = C.xQ +(3:itg.k-2)'*ones(1,par.n+1);      % Ceq(par.n+3:par.n+itg.j+1) w.r.t. the mole fractions
        C.Sx = aux+4:aux+8;                             % Ceq(end = par.n+itg.j+2)
        C.xQ = C.xQ(:)'; C.xx = C.xx(:)';

        % Definition of the values of the Jacobian at the positions
        V.Q = [On1N(1:end-1) -1];                       % Derivatives of Ceq(1) w.r.t. the flowrates
        V.p = [1; -1]*ones(1,par.n);                    % Derivatives of Ceq(2:par.n+1) w.r.t. the pressure
        V.E = [x(3,1:par.n) -x(3,end)...                  Derivatives of Ceq(par.n+2) w.r.t. the flowrates
               x(1,1:par.n) -x(1,end)];                 % Derivatives of Ceq(par.n+2) w.r.t. the temperatures
        V.xQ = [x(4:end-1,1:par.n) -x(4:end-1,end)];    % Derivative of Ceq(par.n+3:par.n+itg.j+1) w.r.t. the flowrates
        V.xx = ones(itg.j-1,1)*[x(1,1:par.n) -x(1,end)];% Derivative of Ceq(par.n+3:par.n+itg.j+1) w.r.t. the mole fractions
        V.Sx = -ones(1,itg.j);                          % Derivative of Ceq(par.n+itg.j+2) w.r.t. mole fractions
        V.p = V.p(:)'; V.xQ = V.xQ(:)'; V.xx = V.xx(:)';
        
        % Definition of the combined vectors
        R.tot = [R.Q R.p R.E R.x  R.x  R.Sx];
        C.tot = [C.Q C.p C.E C.xQ C.xx C.Sx];
        V.tot = [V.Q V.p V.E V.xQ V.xx V.Sx]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);
        
        % Definition of the sparse Jacobian matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,par.n+itg.j+2,itg.k*(par.n+1));
end