function r=Compressor(Var,Cmp,Si,So,Mode)
% Compressor model 
% Var   : Variables
% Cmpr  : Compressor 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)]

% Equation scaling matrix
Neq=NC+4;
S=sparse(diag([
  2.1606e+000
  2.0533e+000
  2.0533e+000
  2.0533e+000
  2.0533e+000
  1.0000e+000
  5.2896e-005
  1.0000e-004
  5.8826e-002
  ]));

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=IxCmp(Cmp);  % Compressor
IxW = Ix(1);      W = Var(IxW);               % Compressor work
IxnC= Ix(2);      nC= Var(IxnC);              % Compressor efficiency
IxTs= Ix(3);      Ts= Var(IxTs);              % Compressor isentropic outlet temperature

if strcmp(Mode,'Mass')
    %Compute or define the mass matrix for this model
    % Steady state model, all elements eq. zero
    M=spalloc(Neq,NVar,0);
    r=S*M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables reciding in this module
    
    Var(Ixyo)= [.97 .03 .0 .0 .0]';              % Composition
    Var(IxFo)= 0.8;                              % Molar component flow      [kmol/s]
    Var(IxTo)= 333;                              % Stream temperature        [C] 
    Var(IxPo)= 10;                               % Outlet pressure           [bar]

    Var(IxW) = 600 ;                             % Compressor work           [kW]
    Var(IxnC)= 0.75;                             % Compressor efficiency
    Var(IxTs)= 320;                              % Compressor isentropic outlet temperature [C]

    r=Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(Neq,NVar,59);
    
    II=speye(NC);
    hi=h(Ti,1);
    ho=h(To,1);
    hs=h(Ts,1);
    dhi_dT=dh_dT(Ti,1);
    dho_dT=dh_dT(To,1);
    dhs_dT=dh_dT(Ts,1);
    
    si=s(Ti,Pi,1);
    so=s(Ts,Po,1);
    dsi_dT=ds_dT(Ti,Pi,1);
    dso_dT=ds_dT(Ts,Po,1);
    dsi_dP=ds_dP(Ti,Pi,1);
    dso_dP=ds_dP(Ts,Po,1);

    
    % 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 (Reversible)
    J(NC+2,IxFi) = yi'*hi;            J(NC+2,IxFo) =-yo'*hs;  
    J(NC+2,Ixyi) = Fi*hi';            J(NC+2,Ixyo) =-Fo*hs';  
    J(NC+2,IxTi) = Fi*yi'*dhi_dT;     J(NC+2,IxTs) =-Fo*yo'*dhs_dT;
    J(NC+2,IxW)  = nC;                J(NC+2,IxnC) = W;
    % Energy balance (Irreversible)
    J(NC+3,IxFi) = yi'*hi;            J(NC+3,IxFo) =-yo'*ho;  
    J(NC+3,Ixyi) = Fi*hi';            J(NC+3,Ixyo) =-Fo*ho';  
    J(NC+3,IxTi) = Fi*yi'*dhi_dT;     J(NC+3,IxTo) =-Fo*yo'*dho_dT;
    J(NC+3,IxW)  = 1;
    % Entropy balance
    J(NC+4,IxFi) = yi'*si;            J(NC+4,IxFo) =-yo'*so;  
    J(NC+4,Ixyi) = Fi*si';            J(NC+4,Ixyo) =-Fo*so';  
    J(NC+4,IxTi) = Fi*yi'*dsi_dT;     J(NC+4,IxTs) =-Fo*yo'*dso_dT;
    J(NC+4,IxPi) = Fi*yi'*dsi_dP;     J(NC+4,IxPo) =-Fo*yo'*dso_dP;
    
    r=S*J;

else
    r=spalloc(Neq,1,Neq);
    hi=h(Ti,1);                                   % [kJ/kmol]
    ho=h(To,1);                                   % [kJ/kmol]
    hs=h(Ts,1);                                   % [kJ/kmol]
    si=s(Ti,Pi,1);                                % [kJ/kmol/K]
    so=s(Ts,Po,1);                                % [kJ/kmol/K]
    
    % 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'*hs+W*nC                  % [kJ]     (Ts) Energy balance  (Reversible)
        Fi*yi'*hi-Fo*yo'*ho+W                     % [kJ]     (To) Energy balance  (Irreversible)
        Fi*yi'*si-Fo*yo'*so                       % [kJ/K]   (Po) Entropy balance (Isentropic process)
       ];               
    
    r=S*r;
end
return
