function Res = CSTR(Var,CSTRNo,Si,So,V,Mode)
% Reformer reactor model
% 02.05.02 (TL)

% Components
% Hydrogen   : H2 
% Gass       : Gs
% Paraffins  : Pf
% Naphthenes : Nt
% Aromatics  : Am

GlobalVars;

I  = ones(1,NC);             % Ones vector
II = speye(NC);              % Identity matrix

% Constants
R  = 8.314;                  % Ideal gas constant                      [kJ/(kmol*K)]
% rc = 0.8;                  % Volume fraction of gas in catalyst bed  [m3/m3]
dc = 4564;                   % Catalyst density                        [kg/m3]
% Vg = V*rc;                 % Gas vlume                               [m3]
mc = V;                      % Mass of catalyst in reactor             [kg]
kf = 0.016;                   % Outflow konstant                        [m^3/bar]
% Cc = 1e3;                  % Catalyst heat capacity                  [MJ/(kg*K)]
% kc = 1e-2;                 % Catalyst heat transfeer constant        [MJ/(K)]

% Stoichimetric matrix of reactions  (Smith)
N=zeros(4,NC);
%         H2 Gs Pf Nt Am
N(1,1:5)=[ 3  0  0 -1  1];   % N <=> A + 3H2     Conversion of naphthenes to aromatics
N(2,1:5)=[-1  0  1 -1  0];   % N + H2 <=> P      Conversion of naphthenes to paraffins
N(3,1:5)=[-2  2  0 -1  0];   % N + 2H2 <=> 2G    Hydrocracking of naphthenes
N(4,1:5)=[-1  2 -1  0  0];   % P + H2 <=> 2G     Hydrocracking of paraffins

% Heat of reaction from Bommannan
Hr=[ 2.2783 -0.44185 -0.51860 -0.56508].*1E5;           % [kJ/kmol]
%Hr=[ 1.7500 -0.44185 -0.51860 -0.56508].*1E5;           % [kJ/kmol]
Hr=-Hr;

Neq=NC+3;                                   % Number of equations
% S=speye(Neq);                               % Equation scaling matrix
% S(NC+2,NC+2)=0.00001;                       % EB scaling
% S(NC+3,NC+3)=0.1;                           % Pressure

if CSTRNo < 11
  S=sparse(diag([1.5481e+000 1.2399e+000 1.2377e+000 1.0272e+000 1.2300e+000 1.0000e+000 7.2553e-006 1.0000e-001])); % For Reactor #1
elseif (CSTRNo > 10) & (CSTRNo < 21)
  S=sparse(diag([1.4371e+000 1.0639e+000 1.0608e+000 8.4631e-001 1.0537e+000 1.0000e+000 5.5152e-006 1.0000e-001])); % For Reactor #2
elseif (CSTRNo > 20) & (CSTRNo < 31)
  S=sparse(diag([1.3713e+000 9.4437e-001 9.3497e-001 5.9944e+000 9.2738e-001 1.0000e+000 9.3445e-006 1.0000e-001])); % For Reactor #3
else
  S=sparse(diag([1.3426e+000 8.7871e-001 8.2977e-001 3.0896e+000 8.4437e-001 1.0000e+000 9.1220e-006 1.0000e-001])); % For Reactor #4
end

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

Ix=IxStr(Si);
Ixyi = Ix(1:NC);    yi = Var(Ixyi);         % Inlet stream komposition                            []
IxFi = Ix(NC+1);    Fi = Var(IxFi);         % Flow into this reactor element                      [kmol/s]
IxTi = Ix(NC+2);    Ti = Var(IxTi);         % Inlet stream temperature                            [K]
IxPi = Ix(NC+3);    Pi = Var(IxPi);         % Inlet stream pressure                               [bar a]

Ix=IxStr(So);
Ixyo = Ix(1:NC);   yo = Var(Ixyo);          % Outlet stream composition
IxFo = Ix(NC+1);   Fo = Var(IxFo);          % Outlet stream molar flow
IxTo = Ix(NC+2);   To = Var(IxTo);          % Outlet stream temperature
IxPo = Ix(NC+3);   Po = Var(IxPo);          % Outlet stream pressure  

Ix=IxCSTR(CSTRNo);
IxCa = Ix(1);      Ca = Var(IxCa);          % Catalyst activity parameter

if strcmp(Mode,'Mass')
    % Compute or define the mass matrix for this model
    M=spalloc(Neq,NVar,NC+2);
    Res=S*M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables reciding in this module

    Var(Ixyo) = [0.78 0.01 .07 .07 .07]';       % Outlet stream composition
    Var(IxFo) = 1.4;                            % Outlet stream molar flow
    Var(IxTo) = 650;                            % Outlet stream temperature
    Var(IxPo) =  10;                             % Outlet stream pressure  

    Var(IxCa) = 1;                              % Catalyst activity parameter

    Res    = Var;
    
elseif strcmp(Mode,'Jacobian')  % Compute Jacobian matrix
  % Reaction rate and it's derivatives
    J=spalloc(Neq,NVar,89);
    
    rr=Rr(yo,To,Po);
    drr_dy=dRr_dy(yo,To,Po);
    drr_dT=dRr_dT(yo,To,Po);
    drr_dP=dRr_dP(yo,To,Po);

    hi =h(Ti,1);                                     % [kJ/kmol]
    ho =h(To,1);                                     % [kJ/kmol]
    dhi_dT=dh_dT(Ti,1);                              % [kJ/kmol]
    dho_dT=dh_dT(To,1);                              % [kJ/kmol]

    % Component balance
    J(1:NC,Ixyi)= Fi*II;             J(1:NC,Ixyo)=-Fo*II+Ca*mc*N'*drr_dy; 
    J(1:NC,IxFi)= yi;                J(1:NC,IxFo)=-yo;
    J(1:NC,IxTo)= Ca*mc*N'*drr_dT;   J(1:NC,IxPo)= Ca*mc*N'*drr_dP;
    J(1:NC,IxCa)= mc*N'*rr;
    % Sum of components equal to one
    J(NC+1,Ixyo)= I;
    % Energy balance
    J(NC+2,IxFi)= yi'*hi;             J(NC+2,IxFo) =-yo'*ho;  
    J(NC+2,Ixyi)= Fi* hi';            J(NC+2,Ixyo) =-Fo*ho'+Ca*mc*Hr*drr_dy;  
    J(NC+2,IxTi)= Fi* yi'*dhi_dT;     J(NC+2,IxTo) =-Fo*yo'*dho_dT+Ca*mc*Hr*drr_dT;
    J(NC+2,IxPo)= Ca*mc*Hr*drr_dP;    J(NC+2,IxCa) = mc*Hr*rr;
    % Inlet pressure    
    J(NC+3,IxPi)=-1;                         J(NC+3,IxPo)=+1-2*(kf*R/100*Fo*To)^2/Po^3;
    J(NC+3,IxTo)= 2*(kf*R/100*Fo/Po)^2*To;   J(NC+3,IxFo)=   2*(kf*R/100*To/Po)^2*Fo;
    
    Res=S*J;

else
    % Compute dxdt
    rr =Rr(yo,To,Po);                               % [kmol/(s*kg)] Reaction rate
    hi=h(Ti,1);                                     % [kJ/kmol]
    ho=h(To,1);                                     % [kJ/kmol]
    
    r=[
      Fi*yi-Fo*yo+Ca*mc*N'*rr;                      % [kmol/s] (n)   Component balance
      I*yo-1                                        % [molfrac](Fo)  Sum of molfractions
      Fi*yi'*hi-Fo*yo'*ho+Ca*mc*Hr*rr;              % [kJ/s]   (To)  Energy balance
      Po-Pi+(kf*R/100*Fo*To/Po)^2                   % [bar]    (Po)  Pressure (0.1 bar pressure drop)
      ];
    Res=S*r;
end
  