function r=Split(Var,Si,So1,So2,Mode)
% Stream split model
% Var   : Variables
% Si    : Inlet stream
% So1   : Outlet stream no. 1
% So2   : Outlet stream no. 2
% Mode  : What to return

GlobalVars;

% Constants
I  = ones(1,NC);

% Equation scaling matrix
Neq=2*NC+5;

S=sparse(diag([
  3.6010e+000
  2.2311e+000
  2.2311e+000
  2.2311e+000
  2.2311e+000
  1.0000e+000
  1.3694e-004
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  2.0000e-003
  1.2500e-001
  1.2500e-001
  ]));

if Unscaled
    S=speye(size(S));
end

% Variables
Ix=IxStr(Si);                 % Inlet stream no. 1
Ixxi = Ix(1:NC);  xi=Var(Ixxi);
IxFi = Ix(NC+1);  Fi=Var(IxFi);
IxTi = Ix(NC+2);  Ti=Var(IxTi);
IxPi = Ix(NC+3);  Pi=Var(IxPi);

Ix=IxStr(So1);                 % Inlet stream no. 2
Ixxo1 = Ix(1:NC);  xo1=Var(Ixxo1);
IxFo1 = Ix(NC+1);  Fo1=Var(IxFo1);
IxTo1 = Ix(NC+2);  To1=Var(IxTo1);
IxPo1 = Ix(NC+3);  Po1=Var(IxPo1);

Ix=IxStr(So2);                                  % Outlet stream
Ixxo2 = Ix(1:NC);  xo2=Var(Ixxo2);
IxFo2 = Ix(NC+1);  Fo2=Var(IxFo2);
IxTo2 = Ix(NC+2);  To2=Var(IxTo2);
IxPo2 = Ix(NC+3);  Po2=Var(IxPo2);


if strcmp(Mode,'Mass')
    % Compute or define the mass matrix for this model
    % Only algebraic equations in this model
    M=spalloc(Neq,NVar,0);
    r=S*M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables reciding in this module
    
    % Out Flow
    Var(Ixxo1)= [.99 .01 .0 .0 .0]';             % Composition
    Var(IxFo1)= 1;                               % Molar component flow      [kmol/s]
    Var(IxTo1)= 473;                             % Stream temperature        [C] 
    Var(IxPo1)= 10;                              % Outlet pressure           [bar]

    Var(Ixxo2)= [.99 .01 .0 .0 .0]';             % Composition
    Var(IxFo2)= 1;                               % Molar component flow      [kmol/s]
    Var(IxTo2)= 473;                             % Stream temperature        [C] 
    Var(IxPo2)= 10;                              % Outlet pressure           [bar]

    r=Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(Neq,NVar,29);
    hi =h(Ti ,1);                                 % [kJ/kmol]
    ho1=h(To1,1);                                 % [kJ/kmol]
    ho2=h(To2,1);                                 % [kJ/kmol]

    dhi_dT =dh_dT(Ti ,1);                         % [kJ/kmol]
    dho1_dT=dh_dT(To1,1);                         % [kJ/kmol]
    dho2_dT=dh_dT(To2,1);                         % [kJ/kmol]

    II=speye(NC);
    % Component balance
    J(1:NC,Ixxi)= Fi*II;             J(1:NC,Ixxo1)=-Fo1*II;               J(1:NC,Ixxo2)=-Fo2*II; 
    J(1:NC,IxFi)= xi;                J(1:NC,IxFo1)=-xo1;                  J(1:NC,IxFo2)=-xo2;  
    % Sum of components equal to one
    J(NC+1,Ixxo1)= I;
    % Energy balance
    J(NC+2,IxFi)= xi'*hi;            J(NC+2,IxFo1)=-xo1'*ho1;             J(NC+2,IxFo2) =-xo2'*ho2;  
    J(NC+2,Ixxi)= Fi*hi';            J(NC+2,Ixxo1)=-Fo1*ho1';             J(NC+2,Ixxo2) =-Fo2*ho2';  
    J(NC+2,IxTi)= Fi*xi'*dhi_dT;     J(NC+2,IxTo1)=-Fo1*xo1'*dho1_dT;     J(NC+2,IxTo2) =-Fo2*xo2'*dho2_dT;
    % Outlet composition
    J(NC+3:2*NC+2,Ixxo1)= II;        J(NC+3:2*NC+2,Ixxo2)=-II; 
    % Outlet temperature
    J(2*NC+3,IxTo1)= 1;              J(2*NC+3,IxTo2)=-1;
    % Inlet pressure
    J(2*NC+4,IxPi)=1;                J(2*NC+4,IxPo1)=-1;
    J(2*NC+5,IxPi)=1;                J(2*NC+5,IxPo2)=-1;

    r=S*J;

else
    r=spalloc(Neq,1,Neq);
    % Compute dxdt
    hi =h(Ti,1);                                 % [kJ/kmol]
    ho1=h(To1,1);                                 % [kJ/kmol]
    ho2=h(To2,1);                                 % [kJ/kmol]
    r =[
        Fi*xi-Fo1*xo1-Fo2*xo2                     % [kmol/s] (xo) Comonent balance
        I*xo1-1                                   % [molfrac](Fo) Sum of molfractions
        Fi*xi'*hi-Fo1*xo1'*ho1-Fo2*xo2'*ho2       % [kJ/s]   (U)  Energy balance
        xo1-xo2                                   % [molfrac]     Equal outlet composition
        To1-To2                                   % [C]           Equal outlet temperature
        Pi-Po1                                    % [bar]    (Pi1)Pressure inlet stream no. 1
        Pi-Po2                                    % [bar]    (Pi1)Pressure inlet stream no. 1
      ];               
    
    r=S*r;
end
return
