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

% Simple Split
% Written by Julian Straus 20150805

% This function defines a simple isenthalpic, isobaric split.

% The state variables are defined as (first inlet and then all outlets)
%       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:
%                   .ratio  Defines the ratios in the outlets, can be also
%                           empty
%                   .Num    Number of the Split

%       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


%       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(1:end)';

% Calculation of the number of outlets and defined splitratios
nratio = rem(numel(x),itg.k);
x_var = x(1:end-nratio);
if nratio > 0
    ratio = x(end-nratio+1:end);
else
    ratio = par.ratio;
end
itg.nout = length(x_var)/itg.k-1;

% Redefinition of the state variables
x_var = reshape(x_var,[itg.k,itg.nout+1]);

% Definition of Auxiliary Variables
RatioLength = length(ratio);
RatioNon0 = ratio(ratio>1);
RatioNon0Length = length(RatioNon0);
RatioSum  = sum(ratio);
Onj1 = ones(itg.j,1);
On1N = ones(1,itg.nout);

% Mass balance
if RatioNon0Length ~= itg.nout 
    Ceq = zeros(1+RatioNon0Length+(itg.k-1)*itg.nout,1);
    Ceq(1) = x_var(1,1) - sum(x_var(1,2:end));
end

% Definition of the inequality constraints imposed by the split ratio
if ~isempty(ratio)
    if RatioLength == itg.nout && RatioSum == 100
        Ceq = zeros(itg.k*itg.nout,1);
        Ceq(1:itg.nout) = x_var(1,2:end)' - x_var(1,1)*ratio/100;
        l = itg.nout+1;
    elseif RatioNon0Length == itg.nout && RatioSum ~= 100
        error(['The ratios defined in Split ',num2str(par.Num), ...
               ' is equal to the number of outlets, but do not sum up to 100.'])
                
    elseif RatioLength > itg.nout
        error(['The number of ratios defined in Split ',num2str(par.Num), ...
               ' is higher than the number of defined outlets.'])
           
%     elseif RatioLength < itg.nout
%         error(['The number of ratios defined in Split ',num2str(par.Num), ...
%                ' is smaller than the number of defined outlets.'])
           
    elseif numel(ratio) > numel(ratio(ratio>=0))
        error(['The ratios defined in Split ',num2str(par.Num), ...
               ' include negative values.'])
           
    else
        l = 2;
        for k = 1:RatioLength
            if ratio(k) > 0
                Ceq(l,1) =  x_var(1,k+1) - x_var(1,1)*ratio(k)/100;
                l = l+1;
            end
        end
    
    end
else 
    l = 2;
end
    
% Calculation of the pressure neutrality
R.p1 = l:itg.k-1:(itg.k-1)*itg.nout;
Ceq(R.p1) = (x_var(2,1)-x_var(2,2:end))';

% Calculation of the temperature neutrality
R.T1 = R.p1+1;
Ceq(R.T1) = (x_var(3,1)-x_var(3,2:end))';

% Calculation of the component mass balance
Xbal = (x_var(4:end,1)*On1N-x_var(4:end,2:end));
R.x1 = (l+2:l+itg.j+1)'*On1N+Onj1*(0:itg.k-1:(itg.nout-1)*(itg.k-1));
Ceq(R.x1) = Xbal;

% Scaling of the Equality Constraints
if isempty(ratio)
    scl.Eq = scl.Eq([1 itg.nout+1:end]);
else
    scl.Eq = scl.Eq([1:1+sum(ratio ~=0) itg.nout+1:end]);
end
Ceq = Ceq.*scl.Eq;

% Jacobian of the functions
switch nargout
    case 2
        if RatioNon0Length ~= itg.nout 
            % Definition of auxiliary variables
            aux = ones(1,size(x_var,2)-1);
            
            % Definition of the indices for the rows (R) and columns (C)
            R.Q = [1 aux];
            C.Q = 1:itg.k:numel(x_var);

            % Definition of the values at the positions
            V.Q = [1,...        Derivative of Ceq(1) w.r.t. Qin
                  -aux];      % Derivative of Ceq(1) w.r.t. Qouts
        end
         
        if ~isempty(ratio)
            % Definition of auxiliary variables
            NElCeqQR = length(Ceq)-(itg.k-1)*itg.nout;
            OnCeq = ones(1,NElCeqQR);
            
            if RatioLength == itg.nout && RatioSum == 100
                R.R = (1:NElCeqQR)'*ones(1,+(nratio>0));
                C.R = [OnCeq (itg.k+1:itg.k:numel(x_var)) numel(x_var)+1:numel(x_var)+nratio];
                V.R = [-ratio'/100 OnCeq -x_var(1,1)*OnCeq/100];
            else
                ZeCeq = zeros(1,NElCeqQR-1);
                n = 1;
                R.R = (2:NElCeqQR)'*ones(1,2+(nratio>0));
                C.R = [OnCeq(1:end-1) kron(ZeCeq,ones(1,1+(nratio>0)))];
                V.R = [ZeCeq  OnCeq(1:end-1)];
                
                for k = 1:RatioLength
                    if ratio(k) > 0
                        C.R(n+NElCeqQR-1)     =  k*itg.k+1;
                        V.R(n)                = -ratio(k)/100;
                        if nratio > 0
                            C.R(n+(NElCeqQR-1)*2) = numel(x_var)+k;
                            V.R(n+(NElCeqQR-1)*2) = -x_var(1,1)/100;
                        end
                        n = n+1;
                    end
                end
            end
            
            if RatioNon0Length ~= itg.nout 
                % Redefinition of the vectors
                R.QR = [R.Q R.R(:)'];
                C.QR = [C.Q C.R];
                V.QR = [V.Q V.R];
            else
                R.QR = R.R(:)';
                C.QR = C.R;
                V.QR = V.R;
            end
            
        else
            R.QR = R.Q;
            C.QR = C.Q;
            V.QR = V.Q;
        end
        
        % Definition of auxiliary values
        On1N = ones(1,itg.nout);
        On1Nj = ones(1,itg.nout*itg.j);
        
        % Definition of the indices for the rows (R) and columns (C) for
        % the pressure (.p), temperature (.T), and component mass (.x)
        % balances
        R.p = [R.p1 R.p1];
        R.T = R.p+1;
        R.x = [R.x1(:)' R.x1(:)'];
        
        C.p = [On1N*2 itg.k+2:itg.k:itg.k*(itg.nout+1)];
        C.T = C.p+1;
        R.x2 = (itg.k+4:2*itg.k)'*On1N+Onj1*(0:itg.k:(itg.nout)*(itg.k-1));
        C.x = [(On1N'*(4:itg.k))' R.x2];
        C.x = C.x(:)';
        
        % Definition of the values of the Jacobian  for the pressure (.p),
        % temperature (.T), and component mass (.x) balances
        V.p = [On1N -On1N]; 
        V.T = V.p;
        V.x = [On1Nj -On1Nj];
        
        % Combination of the vectors
        R.tot = [R.QR R.p R.T R.x];
        C.tot = [C.QR C.p C.T C.x];
        V.tot = [V.QR V.p V.T 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(x_var)+nratio);
end