function r=ConstraintCalc(Var,Mode)
% Calculation of constraint values

GlobalVars;

Neq=7;                             % Number of equations
I=ones(1,NC);

% Scaling
% S=speye(Neq); 
% S(1,1)         =0.01;              % RON scaling
% S(2,2)         = 1;                % H2/HC scaling
% S(3,3)         =0.001;             % Feed massflow
% S(4,4)         =0.001;             % Recyclegas massflow
% S(5,5)         =0.001;             % Gas product massflow
% S(6,6)         =0.001;             % Reformate product massflow
% S(7,7)         =0.001;             % H2 product massflow

S=sparse(diag([
  1.0000e-002
  1.2500e-001
  2.5000e-003
  5.0000e-003
  1.0000e-002
  1.2500e-002
  3.3333e-001
  ]));

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

% Variables

% Feed stream
Ix=IxStr(1);
Ixxi = Ix(1:NC);   xi=Var(Ixxi);
IxFi = Ix(NC+1);   Fi=Var(IxFi);

% G+H2 product stream
Ix=IxStr(53);
IxxG = Ix(1:NC);   xG=Var(IxxG);
IxFG = Ix(NC+1);   FG=Var(IxFG);

% reformate product stream
Ix=IxStr(52);
IxxR = Ix(1:NC);   xR=Var(IxxR);
IxFR = Ix(NC+1);   FR=Var(IxFR);

% recycle gas stream
Ix=IxStr(55);
IxxC = Ix(1:NC);   xC=Var(IxxC);
IxFC = Ix(NC+1);   FC=Var(IxFC);

% Constraint values
Ix=IxConstr(1);
IxRON  = Ix(1);      RON  = Var(IxRON);    % Reformate RON
IxH2HC = Ix(2);      H2HC = Var(IxH2HC);   % H2HC ratio in reactor feed
IxmFi  = Ix(3);      mFi  = Var(IxmFi);    % Feed mass flow
IxmFC  = Ix(4);      mFC  = Var(IxmFC);    % Compressor mass flow
IxmFG  = Ix(5);      mFG  = Var(IxmFG);    % Gas product mass flow
IxmFR  = Ix(6);      mFR  = Var(IxmFR);    % Reformate mass flow
IxmFH  = Ix(7);      mFH  = Var(IxmFH);    % Hydrogen mass flow

if strcmp(Mode,'Mass')
    %Compute or define the mass matrix for this model
    M=spalloc(Neq,Nx,0);
    r = M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables receding in this module
    Var(IxRON)  = 100;
    Var(IxH2HC) = 3.3;
    Var(IxmFi)  = 85;
    Var(IxmFC)  =  8;
    Var(IxmFG)  =  4;
    Var(IxmFR)  = 80;
    Var(IxmFH)  =  3;
    
    r = Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(Neq,NVar,1+NC);
    II=speye(NC);
    H2=[1 0 0 0 0]';
    G =[0 1 0 0 0]';
    J(1,IxRON) =1; J(1,IxxR(5))=-51.3848041;
    J(2,IxH2HC)=1; J(2,IxFC)=-xC(1)/Fi;   J(2,IxxC(1))=-FC/Fi; J(2,IxFi)= FC*xC(1)/Fi^2;
    J(3,IxmFi) =1; J(3,IxFi)=-3.6*MW(xi); J(3,Ixxi) =-3.6*dMW(xi)*Fi;
    J(4,IxmFC) =1; J(4,IxFC)=-3.6*MW(xC); J(4,IxxC) =-3.6*dMW(xC)*FC;
    J(5,IxmFG) =1; J(5,IxFG)=-3.6*MW(xG); J(5,IxxG) =-3.6*dMW(G) *FG;
    J(6,IxmFR) =1; J(6,IxFR)=-3.6*MW(xR); J(6,IxxR) =-3.6*dMW(xR)*FR;
    J(7,IxmFH) =1; J(7,IxFG)=-3.6*MW(H2)*xG(1); J(7,IxxG) =-3.6*dMW(H2)*FG.*H2';
    
    r = S*J;

else
    % Compute dxdt
     H2=[1 0 0 0 0]';
     G =[0 1 0 0 0]';

    r= [ ...
       RON-63.64667849-51.3848041*xR(5) % RON of reformat product
       H2HC-FC*xC(1)/Fi           % H2/HC ratio
       mFi-3.6*MW(xi)*Fi          % Feed massflow
       mFC-3.6*MW(xC)*FC          % Recyclegas massflow
       mFG-3.6*MW(xG)*FG          % D-1405 vapor product massflow
       mFR-3.6*MW(xR)*FR          % Reformate product massflow
       mFH-3.6*MW(H2)*xG(1)*FG    % Hydrogen product massflow
       ];
    r = full(S*r);
    
end

return
