function r=E1402_D1405(Var,Hx,Shiv,Shil,Shov,Shol,Mode)
% Model of E-1402 and D-1405
% Cold side feed is a mixture of vapour recycle gas and liquid feed. Cold side outlet is in vapour phase
% Hot side has vapour feed. Hot side outlet has a vapour and a liquid phase
% Heat transfer is a lumped model where e-Ntu model is used based on actual temperatures.
% Streams
% Shsiv    : Hot side inlet vapour phase
% Shsil    : Hot side inlet liquid phase
% Shsov    : Hot side outlet vapour phase
% Shsol    : Hot side outlet liquid phase

 GlobalVars; 
 I  = ones(1,NC);
 A  = 10;                           % Heat surface area

% Properties of water
% Cp(300K)=75.35 [J/mol*K]           NIST Chemistry WebBook
% Cp(400K)=76.74 [J/mol*K]           NIST Chemistry WebBook
% MW      =18.02 [g/mol]
Cpcw = 76; % [kJ/kmol*K]

Neq=2*NC+8;                        % Number of equations

% Scaling
S=sparse(diag([
  8.1473e-006
  2.1412e-005
  1.0803e+000
  1.0693e+000
  4.2462e+000
  4.2462e+000
  5.5117e+000
  1.0000e+000
  1.0000e+000
  1.0000e-004
  1.2338e-002
  1.6729e-001
  1.0000e+000
  1.0000e+000
  1.0000e+000
  2.0000e-003
  1.2500e-001
  1.2500e-001
  ]));
  
if Unscaled
    S=speye(size(S));
end

% Variables
% Liquid inlet stream (cold side)
Ix=IxStr(Shil);
Ixxhil = Ix(1:NC);   xhil = Var(Ixxhil);
IxFhil = Ix(NC+1);   Fhil = Var(IxFhil);
IxThil = Ix(NC+2);   Thil = Var(IxThil);
IxPhil = Ix(NC+3);   Phil = Var(IxPhil);

% Vapour inlet stream (cold side)
Ix=IxStr(Shiv);
Ixxhiv = Ix(1:NC);   xhiv = Var(Ixxhiv);
IxFhiv = Ix(NC+1);   Fhiv = Var(IxFhiv);
IxThiv = Ix(NC+2);   Thiv = Var(IxThiv);
IxPhiv = Ix(NC+3);   Phiv = Var(IxPhiv);

% Vapour outlet stream (cold side)
Ix=IxStr(Shol);
Ixxhol = Ix(1:NC);   xhol = Var(Ixxhol);
IxFhol = Ix(NC+1);   Fhol = Var(IxFhol);
IxThol = Ix(NC+2);   Thol = Var(IxThol);
IxPhol = Ix(NC+3);   Phol = Var(IxPhol);

% Vapour inlet stream (hot side)
Ix=IxStr(Shov);
Ixxhov = Ix(1:NC);   xhov = Var(Ixxhov);
IxFhov = Ix(NC+1);   Fhov = Var(IxFhov);
IxThov = Ix(NC+2);   Thov = Var(IxThov);
IxPhov = Ix(NC+3);   Phov = Var(IxPhov);

% Heatexchanger variables
Ix=IxHex3(Hx);         
IxQ    = Ix(1);       Q    = Var(IxQ);
IxFcw  = Ix(2);       Fcw  = Var(IxFcw);
IxTcwi = Ix(3);       Tcwi = Var(IxTcwi);
IxTcwo = Ix(4);       Tcwo = Var(IxTcwo);
IxU    = Ix(5);       U    = Var(IxU);

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(Ixxhov) = [.99 .01 .0 .0 .0]';                 % Vapour stream composition   
    Var(IxFhov) = 1;                                   % Vapour stream component flow     [kmol/s]
    Var(IxThov) = 303;                                 % Vapour stream temperature        [C] 
    Var(IxPhov) = 6;                                   % Vapour stream pressure           [bar a]

    Var(Ixxhol) = [.0 .0 .3 .3 .4]';                   % Liquid stream composition   
    Var(IxFhol) = 1;                                   % Liquid stream component flow     [kmol/s]
    Var(IxThol) = 303;                                 % Liquid stream temperature        [C] 
    Var(IxPhol) = 6;                                   % Liquid stream pressure           [bar a]

    Var(IxQ)    =  7000 ;                              % Condenser duty                   [kW]
    Var(IxFcw)  =  8;                                  % Coolingwater flow                [kmol/s]
    Var(IxTcwi) =  288;                                % Coolingwater inlet temperature   [C]
    Var(IxTcwo) =  300;                                % Coolingwater outlet temperaure   [C]
    Var(IxU)    =  200;                                % Heat transfer coefficient
    
    r = Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(Neq,NVar,2*NC^2+8*NC+9);
    II=speye(NC);

    hhiv =h(Thiv,1);
    hhil =h(Thil,0);
    hhov =h(Thov,1);
    hhol =h(Thol,0);

    dhhiv_dT =dh_dT(Thiv,1);
    dhhil_dT =dh_dT(Thil,0);
    dhhov_dT =dh_dT(Thov,1);
    dhhol_dT =dh_dT(Thol,0);

    % Heatexchanger calculations
    % Inlet and outlet enthalpies and derivatives
    hhi       = (Fhil*xhil'*hhil+Fhiv*xhiv'*hhiv)/(Fhil+Fhiv); 
    dhhi_dFhil=-Fhiv*(xhil'*hhil-xhiv'*hhiv)/(Fhil+Fhiv)^2; %???
    dhhi_dFhiv=-Fhil*(xhil'*hhil-xhiv'*hhiv)/(Fhil+Fhiv)^2;
    dhhi_dxhil= Fhil*hhil/(Fhil+Fhiv);
    dhhi_dxhiv= Fhiv*hhiv/(Fhil+Fhiv);
    dhhi_dThil=-Fhil*xhil'*dhhil_dT/(Fhil+Fhiv); %???
    dhhi_dThiv= Fhiv*xhiv'*dhhiv_dT/(Fhil+Fhiv);

    hho       = (Fhol*xhol'*hhol+Fhov*xhov'*hhov)/(Fhol+Fhov);
    dhho_dFhol= Fhov*(xhol'*hhol-xhov'*hhov)/(Fhol+Fhov)^2;
    dhho_dFhov=-Fhol*(xhol'*hhol-xhov'*hhov)/(Fhol+Fhov)^2;
    dhho_dxhol= Fhol*hhol/(Fhol+Fhov);
    dhho_dxhov= Fhov*hhov/(Fhol+Fhov);
    dhho_dThol= Fhol*xhol'*dhhol_dT/(Fhol+Fhov);
    dhho_dThov= Fhov*xhov'*dhhov_dT/(Fhol+Fhov);

    
    % e-Ntu model
    Cc  = Fcw*Cpcw;
    Ch  = Fhiv*(hho-hhi)/(Thov-Thiv);
    
    Cmin= min(Cc,Ch);
    Cmax= max(Cc,Ch);
    
    Rc  = Cmin/Cmax;
    Ntu = U*A/Cmin;
    
    % Pure countercurrent flow
    e  = (1-exp(-Ntu*(1-Rc)))/(1-Rc*exp(-Ntu*(1-Rc)));

    if Cmin==Cc
       %Cmin = Cc and Cmax=Ch
       dCmin_dFcw= Cpcw;
       dCmin_dTcwi= 0;
       dCmin_dTcwo= 0;
       
       dCmin_dFhiv= 0; 
       dCmin_dFhil= 0; 
       dCmin_dFhov= 0; 
       dCmin_dFhol= 0; 
       dCmin_dxhiv= I*0;                   
       dCmin_dxhil= I*0;                   
       dCmin_dxhov= I*0;  
       dCmin_dxhol= I*0;  
       dCmin_dThiv = 0;
       dCmin_dThil = 0;
       dCmin_dThov = 0;
       dCmin_dThol = 0;
       
       dCmax_dFcw = 0;
       dCmax_dTcwi= 0;
       dCmax_dTcwo= 0;
    
       dCmax_dFhiv= ((hho-hhi)-Fhiv*dhhi_dFhiv)/(Thov-Thiv);
       dCmax_dFhil= Fhiv*dhhi_dFhil/(Thov-Thiv);
       dCmax_dFhov= Fhiv*dhho_dFhov/(Thov-Thiv);
       dCmax_dFhol= Fhiv*dhho_dFhol/(Thov-Thiv);
       
       dCmax_dxhiv=-Fhiv*dhhi_dxhiv'/(Thov-Thiv);                   
       dCmax_dxhil=-Fhiv*dhhi_dxhil'/(Thov-Thiv);                   
       dCmax_dxhov= Fhiv*dhho_dxhov'/(Thov-Thiv);                   
       dCmax_dxhol= Fhiv*dhho_dxhol'/(Thov-Thiv);                   
       
       dCmax_dThiv = Fhiv*(-dhhi_dThiv/(Thov-Thiv)+(hho-hhi)/(Thov-Thiv)^2);
       dCmax_dThil = Fhiv*  dhhi_dThil/(Thov-Thiv);
       dCmax_dThov = Fhiv*( dhho_dThov/(Thov-Thiv)-(hho-hhi)/(Thov-Thiv)^2);
       dCmax_dThol = Fhiv*  dhho_dThol/(Thov-Thiv);
       
    else
%        %Cmin = Ch and Cmax=Cc
       dCmax_dFcw= Cpcw;
       dCmax_dTcwi= 0;
       dCmax_dTcwo= 0;
            
       dCmax_dFhiv= 0; 
       dCmax_dFhil= 0; 
       dCmax_dFhov= 0; 
       dCmax_dFhol= 0; 
       dCmax_dxhiv= I*0;                   
       dCmax_dxhil= I*0;                   
       dCmax_dxhov= I*0;  
       dCmax_dxhol= I*0;  
       dCmax_dThiv = 0;
       dCmax_dThil = 0;
       dCmax_dThov = 0;
       dCmax_dThol = 0;
       
       dCmin_dFcw = 0;
       dCmin_dTcwi= 0;
       dCmin_dTcwo= 0;
    
       dCmin_dFhiv= ((hho-hhi)-Fhiv*dhhi_dFhiv)/(Thov-Thiv);
       dCmin_dFhil= Fhiv*dhhi_dFhil/(Thov-Thiv);
       dCmin_dFhov= Fhiv*dhho_dFhov/(Thov-Thiv);
       dCmin_dFhol= Fhiv*dhho_dFhol/(Thov-Thiv);
       
       dCmin_dxhiv=-Fhiv*dhhi_dxhiv'/(Thov-Thiv);                   
       dCmin_dxhil=-Fhiv*dhhi_dxhil'/(Thov-Thiv);                   
       dCmin_dxhov= Fhiv*dhho_dxhov'/(Thov-Thiv);                   
       dCmin_dxhol= Fhiv*dhho_dxhol'/(Thov-Thiv);                   
       
       dCmin_dThiv = Fhiv*(-dhhi_dThiv/(Thov-Thiv)+(hho-hhi)/(Thov-Thiv)^2);
       dCmin_dThil = Fhiv*  dhhi_dThil/(Thov-Thiv);
       dCmin_dThov = Fhiv*( dhho_dThov/(Thov-Thiv)-(hho-hhi)/(Thov-Thiv)^2);
       dCmin_dThol = Fhiv*  dhho_dThol/(Thov-Thiv);

    end
    dRc_dCmin = 1/Cmax;           dRc_dCmax=-Cmin/Cmax^2;
    dNtu_dCmin=-U*A/Cmin^2;
    dNtu_dU   = A/Cmin;
    de_dNtu   =(Rc-1)^2*exp(-Ntu*(1-Rc))/(Rc*exp(-Ntu*(1-Rc))-1)^2;   
    de_dRc    =exp(-Ntu*(1-Rc))*(1-Ntu+Rc*Ntu-exp(-Ntu*(1-Rc)))/(Rc*exp(-Ntu*(1-Rc))-1)^2;


    % Coolingwater side
    J(1,IxFcw) = Cpcw*(Tcwo-Tcwi);
    J(1,IxTcwi)=-Fcw*Cpcw;        J(1,IxTcwo)= Fcw*Cpcw;
    J(1,IxQ)   = -1;
    
    %Heat transfeer
    J(2,IxQ)   =- 1;
    J(2,IxFcw)  =  (de_dNtu*dNtu_dCmin*dCmin_dFcw  + de_dRc*(dRc_dCmin*dCmin_dFcw  + dRc_dCmax*dCmax_dFcw ))*Cmin*(Thiv-Tcwi)+e*dCmin_dFcw *(Thiv-Tcwi); 
    J(2,IxTcwi) =  (de_dNtu*dNtu_dCmin*dCmin_dTcwi + de_dRc*(dRc_dCmin*dCmin_dTcwi + dRc_dCmax*dCmax_dTcwi))*Cmin*(Thiv-Tcwi)+e*dCmin_dTcwi*(Thiv-Tcwi)-e*Cmin;;
    J(2,IxTcwo) =  (de_dNtu*dNtu_dCmin*dCmin_dTcwo + de_dRc*(dRc_dCmin*dCmin_dTcwo + dRc_dCmax*dCmax_dTcwo))*Cmin*(Thiv-Tcwi)+e*dCmin_dTcwo*(Thiv-Tcwi);
    J(2,Ixxhiv) = ((de_dNtu*dNtu_dCmin*dCmin_dxhiv + de_dRc*(dRc_dCmin*dCmin_dxhiv + dRc_dCmax*dCmax_dxhiv))*Cmin*(Thiv-Tcwi)+e*dCmin_dxhiv*(Thiv-Tcwi)); 
    J(2,Ixxhil) = ((de_dNtu*dNtu_dCmin*dCmin_dxhil + de_dRc*(dRc_dCmin*dCmin_dxhil + dRc_dCmax*dCmax_dxhil))*Cmin*(Thiv-Tcwi)+e*dCmin_dxhil*(Thiv-Tcwi)); 
    J(2,Ixxhov) = ((de_dNtu*dNtu_dCmin*dCmin_dxhov + de_dRc*(dRc_dCmin*dCmin_dxhov + dRc_dCmax*dCmax_dxhov))*Cmin*(Thiv-Tcwi)+e*dCmin_dxhov*(Thiv-Tcwi));  
    J(2,Ixxhol) = ((de_dNtu*dNtu_dCmin*dCmin_dxhol + de_dRc*(dRc_dCmin*dCmin_dxhol + dRc_dCmax*dCmax_dxhol))*Cmin*(Thiv-Tcwi)+e*dCmin_dxhol*(Thiv-Tcwi));  
    J(2,IxFhiv) =  (de_dNtu*dNtu_dCmin*dCmin_dFhiv + de_dRc*(dRc_dCmin*dCmin_dFhiv + dRc_dCmax*dCmax_dFhiv))*Cmin*(Thiv-Tcwi)+e*dCmin_dFhiv*(Thiv-Tcwi); 
    J(2,IxFhil) =  (de_dNtu*dNtu_dCmin*dCmin_dFhil + de_dRc*(dRc_dCmin*dCmin_dFhil + dRc_dCmax*dCmax_dFhil))*Cmin*(Thiv-Tcwi)+e*dCmin_dFhil*(Thiv-Tcwi); 
    J(2,IxFhov) =  (de_dNtu*dNtu_dCmin*dCmin_dFhov + de_dRc*(dRc_dCmin*dCmin_dFhov + dRc_dCmax*dCmax_dFhov))*Cmin*(Thiv-Tcwi)+e*dCmin_dFhov*(Thiv-Tcwi); 
    J(2,IxFhol) =  (de_dNtu*dNtu_dCmin*dCmin_dFhol + de_dRc*(dRc_dCmin*dCmin_dFhol + dRc_dCmax*dCmax_dFhol))*Cmin*(Thiv-Tcwi)+e*dCmin_dFhol*(Thiv-Tcwi); 
    J(2,IxThiv) =  (de_dNtu*dNtu_dCmin*dCmin_dThiv + de_dRc*(dRc_dCmin*dCmin_dThiv + dRc_dCmax*dCmax_dThiv))*Cmin*(Thiv-Tcwi)+e*dCmin_dThiv*(Thiv-Tcwi) + e*Cmin;
    J(2,IxThil) =  (de_dNtu*dNtu_dCmin*dCmin_dThil + de_dRc*(dRc_dCmin*dCmin_dThil + dRc_dCmax*dCmax_dThil))*Cmin*(Thiv-Tcwi)+e*dCmin_dThil*(Thiv-Tcwi);
    J(2,IxThol) =  (de_dNtu*dNtu_dCmin*dCmin_dThol + de_dRc*(dRc_dCmin*dCmin_dThol + dRc_dCmax*dCmax_dThol))*Cmin*(Thiv-Tcwi)+e*dCmin_dThol*(Thiv-Tcwi);
    J(2,IxThov) =  (de_dNtu*dNtu_dCmin*dCmin_dThov + de_dRc*(dRc_dCmin*dCmin_dThov + dRc_dCmax*dCmax_dThov))*Cmin*(Thiv-Tcwi)+e*dCmin_dThov*(Thiv-Tcwi);
    J(2,IxU)    =   de_dNtu*dNtu_dU*Cmin*(Thiv-Tcwi);
    
% -----------Flash-----------------------
    % Component balance
    J(3:NC+2,Ixxhiv) = Fhiv*II;        J(3:NC+2,Ixxhil) = Fhil*II;        J(3:NC+2,Ixxhov) =-Fhov*II;          J(3:NC+2,Ixxhol) =-Fhol*II; 
    J(3:NC+2,IxFhiv) = xhiv;           J(3:NC+2,IxFhil) = xhil;           J(3:NC+2,IxFhov) =-xhov;             J(3:NC+2,IxFhol) =-xhol;  
    % Sum of components equal to one
    J(NC+3,Ixxhov)= I;                  
    J(NC+4,Ixxhol)= I;
    % Energy balance
    J(NC+5,IxFhiv) = xhiv'*hhiv;          J(NC+5,IxFhil) = xhil'*hhil;            J(NC+5,IxFhov) =-xhov'*hhov;            J(NC+5,IxFhol) =-xhol'*hhol;
    J(NC+5,Ixxhiv) = Fhiv*hhiv';          J(NC+5,Ixxhil) = Fhil*hhil';            J(NC+5,Ixxhov) =-Fhov*hhov';            J(NC+5,Ixxhol) =-Fhol*hhol';
    J(NC+5,IxThiv) = Fhiv*xhiv'*dhhiv_dT; J(NC+5,IxThil) = Fhil*xhil'*dhhil_dT;   J(NC+5,IxThov) =-Fhov*xhov'*dhhov_dT;   J(NC+5,IxThol) =-Fhol*xhol'*dhhol_dT;
    J(NC+5,IxQ) =-1;
    % Equlibrium 
    J(NC+6:2*NC+5,Ixxhov) = II;                    J(NC+6:2*NC+5,Ixxhol) =-K(Thov,Phov); 
    J(NC+6:2*NC+5,IxThov) =-dK_dT(Thov,Phov)*xhol; J(NC+6:2*NC+5,IxPhov) =-dK_dP(Thov,Phov)*xhol;
    % Outlet stream temperature
    J(2*NC+6,IxThov) = 1;   J(2*NC+6,IxThol) =-1;
    % Feed stream pressure
    J(2*NC+7,IxPhiv) = 1;  J(2*NC+7,IxPhov) =-1;
    J(2*NC+8,IxPhil) = 1;  J(2*NC+8,IxPhov) =-1;

    r = S*J;
     
else
    % Compute dxdt
    r=spalloc(Neq,1,Neq);
    % Compute enthalpies
    hhiv =h(Thiv,1);
    hhil =h(Thil,0);
    hhov =h(Thov,1);
    hhol =h(Thol,0);

    dhhiv_dT =dh_dT(Thiv,1);
    dhhil_dT =dh_dT(Thil,0);
    dhhov_dT =dh_dT(Thov,1);
    dhhol_dT =dh_dT(Thol,0);

    hhi       = (Fhil*xhil'*hhil+Fhiv*xhiv'*hhiv)/(Fhil+Fhiv); 
    hho       = (Fhol*xhol'*hhol+Fhov*xhov'*hhov)/(Fhol+Fhov);

    % e-Ntu model
    Cc  = Fcw*Cpcw;
    Ch  = Fhiv*(hho-hhi)/(Thov-Thiv);
    
    Cmin= min(Cc,Ch);
    Cmax= max(Cc,Ch);
    
    Rc  = Cmin/Cmax;
    Ntu = U*A/Cmin;
    
    % Pure countercurrent flow
    e  = (1-exp(-Ntu*(1-Rc)))/(1-Rc*exp(-Ntu*(1-Rc)));

    r =[...
        % Hx cold side and heat transfeer equations
        Fcw*Cpcw*(Tcwo-Tcwi)-Q                              % [kJ/s]    Energy balance
        e*Cmin*(Thiv-Tcwi)-Q                                % [kJ/s]    Heat transfeer
        % Flash equatiions
        Fhiv*xhiv+Fhil*xhil-Fhov*xhov-Fhol*xhol             % [kmol/s]  Component balance
        I*xhov-1                                            % [molfrac] Sum of components
        I*xhol-1                                            % [molfrac] Sum of components
        Fhiv*xhiv'*hhiv+Fhil*xhil'*hhil-Fhov*xhov'*hhov-Fhol*xhol'*hhol-Q   
                                                            % [kJ/s]    Energy balance
        xhov-K(Thov,Phov)*xhol                              % [molfrac] Equilibrium 
        Thov-Thol                                           % [C]       Outlet stream temperature
        Phiv-Phov                                           % [bar]     Feed stream pressure
        Phil-Phov                                           % [bar]     Feed stream pressure
       ];               

    r = S*r;
    
end

return
