function [Ceq, dCeqdx] = Sep(xs,~,scl,~,ther,itg)

% Simple Separator Unit
% Written by Julian Straus 20160317

% This function is a simple separator, in which the inlet stream is
% separated into a gas and liquid stream. It operates at exactly the same
% pressure and temnperature as the inlet stream, no heat transfer is performed.
%
% The equilibrium constant is calculated using Raoult's law for ammonia and
% Henry's law for the other ocmponents

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

%       ther:       Thermodynamical values to calculate henry coefficients of
%                   the supercritical components, must include:
%                   .Hi          Henry coefficient, i = 1-3

%       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

% Defintion of multiples of itg.k
k1 = itg.k;
k2 = 2*itg.k;

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

% Redefinition of the state variables
xs = reshape(xs,[itg.k,3]);

F = xs(1,1); D = xs(1,2); B = xs(1,3); 
p = xs(2,:);
T = xs(3,:) + 273.15;
z = xs(4:end,1); y = xs(4:end,2); x = xs(4:end,3);

% Calulation of the equilibrium constant, Henry's law for supercritical
% components, Raoult's law for subcritical ammonia
H = 101325*exp(ther.H1 + ther.H2/T(1) + ther.H3/T(1)^2);
ps = 101325*(-0.114397*10^3 + 1.24673*T(1) - 0.353366*10^-2*T(1)^2 - ...
               0.304684*10^-5*T(1)^3 + 0.186446*10^-7*T(1)^4);
Keq = [H(1:2);ps;H(3:4)];

% Predefinition of the equality vector
Ceq = zeros(2*itg.j+2,1);

% Calculation of the overall mass balance
Ceq(1) = (F - D - B);

% Calculation of the pressure neutrality
Ceq(2) = (p(1)-p(2));
Ceq(4) = (p(1)-p(3));

% Calculation of the temperature neutrality
Ceq(3) = (T(1)-T(2));
Ceq(5) = (T(1)-T(3));

% Calculation of the component mass balance
Ceq(6:itg.j+4) = (F*z(1:end-1) - D*y(1:end-1) - B*x(1:end-1));

% Calculation of the liquid-vapour equilibrium
Ceq(itg.j+5:2*itg.j+4) = (y*p(1) - Keq.*x);

% Calculation of the sum of mole fractions
Ceq(2*itg.j+5) = 1-sum(y);   	% Gas phase
Ceq(2*itg.j+6) = 1-sum(x);      % Liquid phase

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

% Calculation of the Jacobians of the functions
switch nargout
    case 2
        % Definition of auxiliary values        
        dHdT = (-ther.H2/T(1)^2 - 2*ther.H3/T(1)^3).*H;
        dpsdT = 101325*(1.24673 - 2*0.353366*10^-2*T(1) - ...
                3*0.304684*10^-5*T(1)^2 + 4*0.186446*10^-7*T(1)^3);
        dKeqdT = [dHdT(1:2)' dpsdT dHdT(3:4)'];
        M1 = 1:itg.k:3*itg.k;
        On14 = ones(1,itg.j-1);
        On1j = ones(1,itg.j);

        % Definition of the indices for the rows (R) and columns (C)
        R.Q = ones(1,3);                        % Ceq(1)
        R.p = ones(2,1)*(2:2:4); R.p = R.p(:)'; % Ceq(2:3)
        R.T = R.p+1;                            % Ceq(4:5)
        R.x = (6:itg.j+4)'*ones(1,3);           % Ceq(5:itg.j+4)
        R.eq = (itg.j+5:2*itg.j+4);             % Ceq(itg.j+5:2*itg.j+4)
        R.Sy = On1j*(2*itg.j+5);                % Ceq(2*itg.j+5)
        R.Sx = On1j*(2*itg.j+6);                % Ceq(2*itg.j+6)
        R.x = R.x(:)';

        C.Q = M1;                                       % Ceq(1) w.r.t. the flowrates
        C.p = [0 k1 0 k2]+2;                            % Ceq(2:3) w.r.t. the pressures
        C.T = C.p+1;                                    % Ceq(4:5) w.r.t. the temperatures
        C.xQ = On14'*M1;                             	% Ceq(5:itg.j+4) w.r.t. the flowrates
        C.xx = C.xQ +(3:itg.k-2)'*ones(1,3);            % Ceq(5:itg.j+4) w.r.t. the mole fractions
        C.eqp = On1j*2;                                 % Ceq(itg.j+5:2*itg.j+4) w.r.t. the pressure
        C.eqT = On1j*3;                                 % Ceq(itg.j+5:2*itg.j+4) w.r.t. the temperature
        C.eqx = On1j'*M1(2:end)+(3:itg.k-1)'*ones(1,2);	% Ceq(itg.j+5:2*itg.j+4) w.r.t. the mole fractions
        C.Sy = (itg.k+4:2*itg.k);                       % Ceq(2*itg.j+5) w.r.t. the mole fractions
        C.Sx = (2*itg.k+4:3*itg.k);                     % Ceq(2*itg.j+6) w.r.t. the mole fractions
        C.xQ = C.xQ(:)'; C.xx = C.xx(:)'; C.eqx = C.eqx(:)';

        % Definition of the values of the Jacobian at the positions
        V.Q = [1 -ones(1,2)];                           % Derivatives of Ceq(1) w.r.t. the flowrates
        V.p = [1 -1 1 -1];                              % Ceq(2:3) w.r.t. the pressures
        V.T = V.p;                                      % Ceq(4:5) w.r.t. the temperatures
        V.xQ = [z(1:end-1)' -y(1:end-1)' -x(1:end-1)']; % Derivatives of Ceq(5:itg.j+4) w.r.t. the flowrates
        V.xx = [F*On14 -D*On14 -B*On14];                % Derivatives of Ceq(5:itg.j+4) w.r.t. the mole fractions
        V.eqp = y';                                     % Derivatives of Ceq(itg.j+5:2*itg.j+4) w.r.t. the pressure
        V.eqT = -dKeqdT.*x';                            % Derivatives of Ceq(itg.j+5:2*itg.j+4) w.r.t. the temperature
        V.eqx = [On1j*p(1) -Keq'];                      % Derivatives of Ceq(itg.j+5:2*itg.j+4) w.r.t. the mole fractions
        V.Sy = -On1j;                                   % Ceq(2*itg.j+5) w.r.t. the mole fractions
        V.Sx = -On1j;                                   % Ceq(2*itg.j+6) w.r.t. the mole fractions
       
        % Definition of the combined vectors
        R.tot = [R.Q R.p R.T R.x  R.x  R.eq  R.eq  R.eq R.eq R.Sy R.Sx];
        C.tot = [C.Q C.p C.T C.xQ C.xx C.eqp C.eqT C.eqx     C.Sy C.Sx];
        V.tot = [V.Q V.p V.T V.xQ V.xx V.eqp V.eqT V.eqx     V.Sy V.Sx]...
                .*scl.Eq(R.tot)'.*scl.Var(C.tot);

        % Definition of the sparse Jacobian matrix
        dCeqdx = sparse(R.tot,C.tot,V.tot,2*itg.j+6,itg.k*3);
end