function r=SSSpecs(Var,Mode)
% Used for testing of unit models
%

GlobalVars;
Tk=273.15;

% Equation scaling matrix
Neq = NC+54;  % All eq. NC+54
% S   = speye(Neq);
%     i=37;
%     S(i,i)   = 0.1;            i=i+1;    % Plo
%                                i=i+NC;   % xi
%     S(i,i)   = 0.1;            i=i+1;    % Fcw
%     S(i,i)   = 0.01;           i=i+1;    % Tcwi
%                                i=i+4;    % Ca19,29,39,49
%     S(i,i)   = 0.01;           i=i+1;    % U2
%     S(i,i)   = 0.1;            i=i+1;    % U3
%     S(i,i)   = 0.01;           i=i+1;    % mFi
%     S(i,i)   = 0.001;          i=i+1;    % Ti
%     S(i,i)   = 0.01;           i=i+1;    % nC
%     S(i,i)   = 0.1;            i=i+1;    % Pvo
%     S(i,i)   = 0.1;            i=i+1;    % mFC
%     S(i,i)   = 0.001;          i=i+1;    % Th1
%     S(i,i)   = 0.001;          i=i+1;    % Th2     
%     S(i,i)   = 0.001;          i=i+1;    % Th3
%     S(i,i)   = 0.001;          i=i+1;    % Th4

S=sparse(diag([
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.2500e-001
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  5.0000e-001
  3.4722e-003
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  4.1667e-003
  1.0000e-002
  2.5000e-003
  2.5445e-003
  2.0000e-002
  1.2500e-001
  5.0000e-003
  1.2937e-003
  1.2937e-003
  1.2937e-003
  1.2937e-003
  ]));

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


% Inlet stream no. 1
Ix=IxStr(1);
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);

% Outlet vapor stream 
Ix=IxStr(53);
Ixxvo = Ix(1:NC);   xvo=Var(Ixxvo);
IxFvo = Ix(NC+1);   Fvo=Var(IxFvo);
IxTvo = Ix(NC+2);   Tvo=Var(IxTvo);
IxPvo = Ix(NC+3);   Pvo=Var(IxPvo);

% Heater 1 outlet stream 
Ix=IxStr(5);
Ixxh1 = Ix(1:NC);   xh1=Var(Ixxh1);
IxFh1 = Ix(NC+1);   Fh1=Var(IxFh1);
IxTh1 = Ix(NC+2);   Th1=Var(IxTh1);
IxPh1 = Ix(NC+3);   Ph1=Var(IxPh1);

% Heater 2 outlet stream 
Ix=IxStr(16);
Ixxh2 = Ix(1:NC);   xh2=Var(Ixxh2);
IxFh2 = Ix(NC+1);   Fh2=Var(IxFh2);
IxTh2 = Ix(NC+2);   Th2=Var(IxTh2);
IxPh2 = Ix(NC+3);   Ph2=Var(IxPh2);
% 
% % Heater 3 outlet stream 
Ix=IxStr(27);
Ixxh3 = Ix(1:NC);   xh3=Var(Ixxh3);
IxFh3 = Ix(NC+1);   Fh3=Var(IxFh3);
IxTh3 = Ix(NC+2);   Th3=Var(IxTh3);
IxPh3 = Ix(NC+3);   Ph3=Var(IxPh3);
% 
% % Heater 4 outlet stream 
Ix=IxStr(38);
Ixxh4 = Ix(1:NC);   xh4=Var(Ixxh4);
IxFh4 = Ix(NC+1);   Fh4=Var(IxFh4);
IxTh4 = Ix(NC+2);   Th4=Var(IxTh4);
IxPh4 = Ix(NC+3);   Ph4=Var(IxPh4);

% Compressor inlet stream
% Ix=IxStr(52);
% Ixyr = Ix(1:NC);   yr=Var(Ixyr);
% IxFr = Ix(NC+1);   Fr=Var(IxFr);
% IxTr = Ix(NC+2);   Tr=Var(IxTr);
% IxPr = Ix(NC+3);   Pr=Var(IxPr);

% Compressor outlet stream 
Ix=IxStr(55);
IxxCmp = Ix(1:NC);   xCmp=Var(IxxCmp);
IxFCmp = Ix(NC+1);   FCmp=Var(IxFCmp);
IxTCmp = Ix(NC+2);   TCmp=Var(IxTCmp);
IxPCmp = Ix(NC+3);   PCmp=Var(IxPCmp);

% Condenser 1 liquid outlet stream 
Ix=IxStr(52);
Ixxlo = Ix(1:NC);   xlo=Var(Ixxlo);
IxFlo = Ix(NC+1);   Flo=Var(IxFlo);
IxTlo = Ix(NC+2);   Tlo=Var(IxTlo);
IxPlo = Ix(NC+3);   Plo=Var(IxPlo);


% HX Feed-Reactor effluent
Ix=IxHex2(1);         
IxQ2  = Ix(1);       Q2  = Var(IxQ2);
IxU2  = Ix(2);       U2  = Var(IxU2);

% HX Coolingwater-Reactor effluent
Ix=IxHex3(1);         
IxQ3   = Ix(1);       Q3   = Var(IxQ3);
IxFcw  = Ix(2);       Fcw  = Var(IxFcw);
IxTcwi = Ix(3);       Tcwi = Var(IxTcwi);
IxTcwo = Ix(4);       Tcwo = Var(IxTcwo);
IxU3   = Ix(5);       U3   = Var(IxU3);

% Reactor 1 catalyst efficiency
Ix=IxCSTR(1);
IxCa10 = Ix(1);      Ca10=Var(IxCa10);
IxCa11 = Ix(1)+1;    Ca11=Var(IxCa11);
IxCa12 = Ix(1)+2;    Ca12=Var(IxCa12);
IxCa13 = Ix(1)+3;    Ca13=Var(IxCa13);
IxCa14 = Ix(1)+4;    Ca14=Var(IxCa14);
IxCa15 = Ix(1)+5;    Ca15=Var(IxCa15);
IxCa16 = Ix(1)+6;    Ca16=Var(IxCa16);
IxCa17 = Ix(1)+7;    Ca17=Var(IxCa17);
IxCa18 = Ix(1)+8;    Ca18=Var(IxCa18);
IxCa19 = Ix(1)+9;    Ca19=Var(IxCa19);

% Reactor 2 catalyst efficiency
Ix=IxCSTR(11);
IxCa20 = Ix(1);      Ca20=Var(IxCa20);
IxCa21 = Ix(1)+1;    Ca21=Var(IxCa21);
IxCa22 = Ix(1)+2;    Ca22=Var(IxCa22);
IxCa23 = Ix(1)+3;    Ca23=Var(IxCa23);
IxCa24 = Ix(1)+4;    Ca24=Var(IxCa24);
IxCa25 = Ix(1)+5;    Ca25=Var(IxCa25);
IxCa26 = Ix(1)+6;    Ca26=Var(IxCa26);
IxCa27 = Ix(1)+7;    Ca27=Var(IxCa27);
IxCa28 = Ix(1)+8;    Ca28=Var(IxCa28);
IxCa29 = Ix(1)+9;    Ca29=Var(IxCa29);

% Reactor 3 catalyst efficiency
Ix=IxCSTR(21);
IxCa30 = Ix(1);      Ca30=Var(IxCa30);
IxCa31 = Ix(1)+1;    Ca31=Var(IxCa31);
IxCa32 = Ix(1)+2;    Ca32=Var(IxCa32);
IxCa33 = Ix(1)+3;    Ca33=Var(IxCa33);
IxCa34 = Ix(1)+4;    Ca34=Var(IxCa34);
IxCa35 = Ix(1)+5;    Ca35=Var(IxCa35);
IxCa36 = Ix(1)+6;    Ca36=Var(IxCa36);
IxCa37 = Ix(1)+7;    Ca37=Var(IxCa37);
IxCa38 = Ix(1)+8;    Ca38=Var(IxCa38);
IxCa39 = Ix(1)+9;    Ca39=Var(IxCa39);

% Reactor 4 catalyst efficiency
Ix=IxCSTR(31);
IxCa40 = Ix(1);      Ca40=Var(IxCa40);
IxCa41 = Ix(1)+1;    Ca41=Var(IxCa41);
IxCa42 = Ix(1)+2;    Ca42=Var(IxCa42);
IxCa43 = Ix(1)+3;    Ca43=Var(IxCa43);
IxCa44 = Ix(1)+4;    Ca44=Var(IxCa44);
IxCa45 = Ix(1)+5;    Ca45=Var(IxCa45);
IxCa46 = Ix(1)+6;    Ca46=Var(IxCa46);
IxCa47 = Ix(1)+7;    Ca47=Var(IxCa47);
IxCa48 = Ix(1)+8;    Ca48=Var(IxCa48);
IxCa49 = Ix(1)+9;    Ca49=Var(IxCa49);

% Heater 1 duty
% Ix=IxHtr(1);
% IxQH1 = Ix(1);      QH1=Var(IxQH1);

% % Heater 2 duty
% Ix=IxHtr(2);
% IxQH2 = Ix(1);      QH2=Var(IxQH2);
% 
% % Heater 3 duty
% Ix=IxHtr(3);
% IxQH3 = Ix(1);      QH3=Var(IxQH3);
% 
% % Heater 4 duty
% Ix=IxHtr(4);
% IxQH4 = Ix(1);      QH4=Var(IxQH4);

% Compressor
Ix=IxCmp(1);                                
IxW = Ix(1);      W = Var(IxW);               % Compressor work
IxnC= Ix(2);      nC= Var(IxnC);              % Compressor efficiency

% 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,NVar,0);
    r=M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables reciding in this module
    Var(Ixyi) = [0 0 .39 .46 .15];     % Composition
    Var(IxFi) = .24;                   % Feed stream component flow                          [kmol/s]
    Var(IxTi) = 393.15;                % Feed stream temperature                             [C] 
    Var(IxPi) = 10;                    % Vapour stream pressure                              [bar a]

  
    r=Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(Neq,NVar,NC+2);
    i=1;
    J(i,IxCa10) = 1;    J( i,IxCa11) =-1; i=i+1;
    J(i,IxCa11) = 1;    J( i,IxCa12) =-1; i=i+1;
    J(i,IxCa12) = 1;    J( i,IxCa13) =-1; i=i+1;
    J(i,IxCa13) = 1;    J( i,IxCa14) =-1; i=i+1;
    J(i,IxCa14) = 1;    J( i,IxCa15) =-1; i=i+1;
    J(i,IxCa15) = 1;    J( i,IxCa16) =-1; i=i+1;
    J(i,IxCa16) = 1;    J( i,IxCa17) =-1; i=i+1;
    J(i,IxCa17) = 1;    J( i,IxCa18) =-1; i=i+1;
    J(i,IxCa18) = 1;    J( i,IxCa19) =-1; i=i+1;

    J(i,IxCa20) = 1;    J(i,IxCa21) =-1; i=i+1;
    J(i,IxCa21) = 1;    J(i,IxCa22) =-1; i=i+1;
    J(i,IxCa22) = 1;    J(i,IxCa23) =-1; i=i+1;
    J(i,IxCa23) = 1;    J(i,IxCa24) =-1; i=i+1;
    J(i,IxCa24) = 1;    J(i,IxCa25) =-1; i=i+1;
    J(i,IxCa25) = 1;    J(i,IxCa26) =-1; i=i+1;
    J(i,IxCa26) = 1;    J(i,IxCa27) =-1; i=i+1;
    J(i,IxCa27) = 1;    J(i,IxCa28) =-1; i=i+1;
    J(i,IxCa28) = 1;    J(i,IxCa29) =-1; i=i+1;

    J(i,IxCa30) = 1;    J(i,IxCa31) =-1; i=i+1;
    J(i,IxCa31) = 1;    J(i,IxCa32) =-1; i=i+1;
    J(i,IxCa32) = 1;    J(i,IxCa33) =-1; i=i+1;
    J(i,IxCa33) = 1;    J(i,IxCa34) =-1; i=i+1;
    J(i,IxCa34) = 1;    J(i,IxCa35) =-1; i=i+1;
    J(i,IxCa35) = 1;    J(i,IxCa36) =-1; i=i+1;
    J(i,IxCa36) = 1;    J(i,IxCa37) =-1; i=i+1;
    J(i,IxCa37) = 1;    J(i,IxCa38) =-1; i=i+1;
    J(i,IxCa38) = 1;    J(i,IxCa39) =-1; i=i+1;

    J(i,IxCa40) = 1;    J(i,IxCa41) =-1; i=i+1;
    J(i,IxCa41) = 1;    J(i,IxCa42) =-1; i=i+1;
    J(i,IxCa42) = 1;    J(i,IxCa43) =-1; i=i+1;
    J(i,IxCa43) = 1;    J(i,IxCa44) =-1; i=i+1;
    J(i,IxCa44) = 1;    J(i,IxCa45) =-1; i=i+1;
    J(i,IxCa45) = 1;    J(i,IxCa46) =-1; i=i+1;
    J(i,IxCa46) = 1;    J(i,IxCa47) =-1; i=i+1;
    J(i,IxCa47) = 1;    J(i,IxCa48) =-1; i=i+1;
    J(i,IxCa48) = 1;    J(i,IxCa49) =-1; i=i+1;

    J(i,IxPlo)       = 1;                i=i+1;
    J(i:NC+i-1,Ixyi) = speye(NC);        i=i+NC;
    J(i,IxFcw)       = 1;                i=i+1;
    J(i,IxTcwi)      = 1;                i=i+1;

    J(i,IxCa19)      = 1;                i=i+1;
    J(i,IxCa29)      = 1;                i=i+1;
    J(i,IxCa39)      = 1;                i=i+1;
    J(i,IxCa49)      = 1;                i=i+1;

    J(i,IxU2)        = 1;                i=i+1;
    J(i,IxU3)        = 1;                i=i+1;
    J(i,IxmFi)       = 1;                i=i+1;
    J(i,IxTi)        = 1;                i=i+1;
    J(i,IxnC)        = 1;                i=i+1;
    J(i,IxPvo)       = 1;                i=i+1;
%    J(i,IxmFC)       = 1;                i=i+1;
    J(i,IxH2HC)      = 1;                i=i+1;
    J(i,IxRON)       = 1;                i=i+1;
    J(i,IxTh1) = 1;  J(i,IxTh2) =-1;     i=i+1;
    J(i,IxTh2) = 1;  J(i,IxTh3) =-1;     i=i+1;
    J(i,IxTh3) = 1;  J(i,IxTh4) =-1;     i=i+1;
%    J(i,IxTh1)       = 1;                i=i+1;
%    J(i,IxTh2)       = 1;                i=i+1;
%    J(i,IxTh3)       = 1;                i=i+1;
%    J(i,IxTh4)       = 1;                i=i+1;

    if EvalScale
    disp('SS Specs')
    VarIx=[IxCa10 IxCa11 IxCa12 IxCa13 IxCa14 IxCa15 IxCa16 IxCa17 IxCa18 IxCa19 ...
           IxCa20 IxCa21 IxCa22 IxCa23 IxCa24 IxCa25 IxCa26 IxCa27 IxCa28 IxCa29 ...
           IxCa30 IxCa31 IxCa32 IxCa33 IxCa34 IxCa35 IxCa36 IxCa37 IxCa38 IxCa39 ...
           IxCa40 IxCa41 IxCa42 IxCa43 IxCa44 IxCa45 IxCa46 IxCa47 IxCa48 IxCa49 ...
           IxPlo Ixyi' IxFcw IxTcwi IxU2 IxU3 IxmFi IxTi IxnC IxPvo IxmFC IxTh1  ...
           IxTh2 IxTh3 IxTh4];
           
    Jus=full(J(:,VarIx));
    P=zeros(size(J));
    VarIxP=[IxCa10 IxCa11 IxCa12 IxCa13 IxCa14 IxCa15 IxCa16 IxCa17 IxCa18 ...
            IxCa20 IxCa21 IxCa22 IxCa23 IxCa24 IxCa25 IxCa26 IxCa27 IxCa28 ...
            IxCa30 IxCa31 IxCa32 IxCa33 IxCa34 IxCa35 IxCa36 IxCa37 IxCa38 ...
            IxCa40 IxCa41 IxCa42 IxCa43 IxCa44 IxCa45 IxCa46 IxCa47 IxCa48 ...
            IxPlo Ixyi' IxFcw IxTcwi IxCa19 IxCa29 IxCa39 IxCa49 IxU2 IxU3 ...
            IxmFi IxTi IxnC IxPvo IxmFC IxTh1 IxTh2 IxTh3 IxTh4];
    
    for i=1:Neq
     P(i,VarIxP(i))=1;
    end
 
    Snew =abs(inv(eye(Neq).*(J*Sv*P')));
    Js=full(Snew*J(:,VarIx)*Sv(VarIx,VarIx));
    disp(['cond reduction Unscaled/Scaled/NewScaling: ' num2str(cond(Jus)) ' / ' num2str(cond(full(S*J*Sv))) ' / ' num2str(cond(Js))])
    disp(['Max/Min values in Js                     : ' num2str(max(max(Js))) ' , ' num2str(min(min(Js)))])
    save Jus.txt Jus -ascii
    save Js.txt Js -ascii
    end

    r=S*J;
    
elseif strcmp(Mode,'b')
%    FeedComp=[0 0 .32 .56 .12]';             %          AI-015 vol%
    FeedComp=[0 0 3.0906e-1 5.5021e-1 1.3097e-1]';             %          AI-015 vol%
   r=  [ ...
       zeros(4*9,1)
       8         % Liquid outlet pressure
       FeedComp  % Feed composition               
       8         % E-1402 coolingwaterflow, 1 kmol/s = 64 T/h (15=960)
       288.15    % E-1402 coolingwater inlet temperature
       1         % Reactor 1 catalyst efficiency
       1         % Reactor 2 catalyst efficiency
       1         % Reactor 3 catalyst efficiency
       1         % Reactor 4 catalyst efficiency
       560       % E-1401 heat transfer coefficient
       200       % E-1402 heat transfer coefficient
       90        % Feed flow
       3.6237e+2 % Feed tempearture
       0.75      % Compressor efficiency
       10        % Vapour outlet presure
%       8        % Compressor recycle flow
       3.78450688493968       % H2HC ratio
       103       % RON
       0         % Heater 1 and 2 delta temperature
       0         % Heater 2 and 3 delta temperature
       0         % Heater 3 and 4 delta temperature
%       790      % Heater 1 outlet temperature
%       790      % Heater 2 outlet temperature
%       790      % Heater 3 outlet temperature
%       790      % Heater 4 outlet temperature
        ];
    r=S*r;
    
else
    FeedComp=[0 0 .32 .56 .12]';             %          AI-015 vol%
    % Compute dxdt
    r=[ ...
       zeros(4*9,1)
       Plo  - 8         % Liquid outlet pressure
       yi   - FeedComp  % Feed composition               
       Fcw  - 8         % E-1402 coolingwaterflow, 1 kmol/s = 64 T/h (15=960)
       Tcwi - 288.15    % E-1402 coolingwater inlet temperature
       Ca19 - 1         % Reactor 1 catalyst efficiency
       Ca29 - 1         % Reactor 2 catalyst efficiency
       Ca39 - 1         % Reactor 3 catalyst efficiency
       Ca49 - 1         % Reactor 4 catalyst efficiency
       U2   - 560       % E-1401 heat transfer coefficient
       U3   - 200       % E-1402 heat transfer coefficient
       mFi  - 85        % Feed flow
       Ti   - 358.5     % Feed tempearture
       nC   - 0.75      % Compressor efficiency
       Pvo  - 7.9      % Vapour outlet presure
       mFC  - 8        % Compressor recycle flow
       Th1  - 790      % Heater 1 outlet temperature
       Th2  - 790      % Heater 2 outlet temperature
       Th3  - 790      % Heater 3 outlet temperature
       Th4  - 790      % Heater 4 outlet temperature
      ];
    
    r=S*r;
end

