function r=Heater(Var,Htr,Si,So,Mode)
% Heater model
% Var   : Variables
% Htr   : Heater number
% Si    : Inlet stream number
% So    : Outlet stream number
% Mode  : What to return

GlobalVars;

% Constants
I  = ones(1,NC);
R  = 8.314;                                   % Ideal gas constant        [kJ/(kmol*K)]
kf = 63;                                      % Outflow konstant          [m^3/bar]
V  = 5;                                       % Volume of heater pipe     [m^3]

% Equation scaling matrix
Neq=NC+3;
S=sparse(diag([1.4796e+000 6.2661e+001 1.3925e+001 8.3547e+000 1.2532e+001 1.4286e+000 1.3751e-005 1.0000e-001]));
if Htr < 4
S=sparse(diag([1.6330e+000 1.3583e+000 1.3583e+000 1.3583e+000 1.3583e+000 1.0000e+000 1.3098e-005 1.0000e-001]));
else
S=sparse(diag([1.3473e+000 8.9697e-001 8.9697e-001 8.9697e+000 8.9697e-001 1.0000e+000 1.2672e-005 1.0000e-001]));
end

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

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

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

Ix=IxHtr(Htr);                                 % Heater
IxQ = Ix(1);      Q = Var(IxQ);

if strcmp(Mode,'Mass')
    %Compute or define the mass matrix for this model
    M=spalloc(Neq,NVar,1);
    M(NC+2,IxU)=1;
    r=S*M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables reciding in this module
    
    % Out Flow
    Var(Ixyo)= [0.78 0.01 .07 .07 .07]';        % Composition
    Var(IxFo)= 1.4;                             % Molar component flow      [kmol/s]
    Var(IxTo)= 650;                             % Stream temperature        [C] 
    Var(IxPo)=  10;                             % Outlet pressure           [bar]
    % Heater
    Var(IxQ) =  8000;                           % Heater duty               [kW]

    r=Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(1,NVar,NC+4);
    
    hi=h(Ti,1);                                 % [kJ/kmol]  , Vapor
    ho=h(To,1);                                 % [kJ/kmol]  , Vapor

    dhi_dT=dh_dT(Ti,1);                         % [kJ/kmol/K], Vapor
    dho_dT=dh_dT(To,1);                         % [kJ/kmol/K], Vapor
    
    II=speye(NC);
    % Component balance
    J(1:NC,Ixyi)= Fi*II;             J(1:NC,Ixyo)=-Fo*II; 
    J(1:NC,IxFi)= yi;                J(1:NC,IxFo)=-yo;  
    % 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';  
    J(NC+2,IxTi)= Fi* yi'*dhi_dT;     J(NC+2,IxTo) =-Fo*yo'*dho_dT;
    J(NC+2,IxQ )= 1;
    % Inlet pressure    
    J(NC+3,IxPi)=1;                   J(NC+3,IxPo)=-1;
    
    r=S*J;

else
    r=spalloc(Neq,1,Neq);
    hi=h(Ti,1);                                   % [MJ/kmol]
    ho=h(To,1);                                   % [MJ/kmol]
    % Compute dxdt
    r =[ ...
        Fi*yi-Fo*yo                               % [kmol/s] (yo) Comonent balance
        I*yo-1                                    % [molfrac](Fo) Sum of molfractions
        Fi*yi'*hi-Fo*yo'*ho+Q                     % [kJ/s]   (U)  Energy balance
        Pi-Po-0.3                                 % [bar]    (Pi) Pressure 
      ];               
    
    r=S*r;
end
