function r=E1401(Var,Hx,Sciv,Scil,Scov,Shiv,Shov,Shol,Mode)
% Model of E-1401
% 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
% Scsiv    : Cold side inlet vapour phase
% Scsil    : Cold side inlet liquid phase
% Scsov    : Cold side outlet vapour phase
% Shsiv    : Hot side inlet vapour phase
% Shsov    : Hot side outlet vapour phase
% Shsol    : Hot side outlet liquid phase

GlobalVars;

A  =  5;                           % Heat surface area
Neq=3*NC+10;                       % Number of equations

% Scaling

S=sparse(diag([
 1.6330e+000
  1.3583e+000
  1.3583e+000
  1.3583e+000
  1.3583e+000
  1.0000e+000
  1.3533e-005
  1.0000e-001
  1.0000e-001
  5.8729e-005
  1.2346e+000
  9.3367e-001
  3.1122e+000
  3.1122e+000
  1.3069e+000
  1.0000e+000
  1.0000e+000
  3.3333e-005
  1.0923e-002
  3.7609e-002
  2.9796e+000
  5.3086e+000
  3.3333e+000
  2.0000e-003
  1.0000e-001
  ]));

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

I=ones(1,NC);

% Variables
% Liquid inlet stream (cold side)
Ix=IxStr(Scil);
Ixxcil = Ix(1:NC);   xcil = Var(Ixxcil);
IxFcil = Ix(NC+1);   Fcil = Var(IxFcil);
IxTcil = Ix(NC+2);   Tcil = Var(IxTcil);
IxPcil = Ix(NC+3);   Pcil = Var(IxPcil);

% Vapour inlet stream (cold side)
Ix=IxStr(Sciv);
Ixxciv = Ix(1:NC);   xciv = Var(Ixxciv);
IxFciv = Ix(NC+1);   Fciv = Var(IxFciv);
IxTciv = Ix(NC+2);   Tciv = Var(IxTciv);
IxPciv = Ix(NC+3);   Pciv = Var(IxPciv);

% Vapour outlet stream (cold side)
Ix=IxStr(Scov);
Ixxcov = Ix(1:NC);   xcov = Var(Ixxcov);
IxFcov = Ix(NC+1);   Fcov = Var(IxFcov);
IxTcov = Ix(NC+2);   Tcov = Var(IxTcov);
IxPcov = Ix(NC+3);   Pcov = Var(IxPcov);

% Vapour inlet stream (hot 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 (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);

% Vapour outlet stream (hot 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);

% Heatexchanger variables
Ix=IxHex2(Hx);         
IxQ   = Ix(1);       Q   = Var(IxQ);
%IxTxi = Ix(2);       Txi = Var(IxTxi)+Tk;
IxU   = Ix(2);       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(Ixxcov) = [.76 .01 .07 .13 .03]';              % Vapour stream composition   
    Var(IxFcov) = 1;                                   % Vapour stream component flow     [kmol/s]
    Var(IxTcov) = 680;                                 % Vapour stream temperature        [C] 
    Var(IxPcov) = 11;                                  % Vapour stream pressure           [bar a]

    Var(Ixxhov) = [.9 .01 .01 .01 .07]';               % Vapour stream composition   
    Var(IxFhov) = .2;                                  % Vapour stream component flow     [kmol/s]
    Var(IxThov) = 350;                                 % Vapour stream temperature        [C] 
    Var(IxPhov) = 8;                                   % Vapour stream pressure           [bar a]

    Var(Ixxhol) = [.0 .0 .1 .2 .7]';                   % Vapour stream composition   
    Var(IxFhol) = 1.4;                                 % Vapour stream component flow     [kmol/s]
    Var(IxThol) = 350;                                 % Vapour stream temperature        [C] 
    Var(IxPhol) = 8;                                   % Vapour stream pressure           [bar a]

    Var(IxQ)    =  37000;                              % Condenser duty                   [kW]
    Var(IxU)    =  560;                                % 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);

    hciv =h(Tciv,1);
    hcil =h(Tcil,0);
    hcov =h(Tcov,1);
    hhiv =h(Thiv,1);
    hhov =h(Thov,1);
    hhol =h(Thol,0);

    dhciv_dT =dh_dT(Tciv,1);
    dhcil_dT =dh_dT(Tcil,0);
    dhcov_dT =dh_dT(Tcov,1);
    dhhiv_dT =dh_dT(Thiv,1);
    dhhov_dT =dh_dT(Thov,1);
    dhhol_dT =dh_dT(Thol,0);

    % Heatexchanger calculations
    % Inlet and outlet enthalpies and derivatives
    hci       = (Fcil*xcil'*hcil+Fciv*xciv'*hciv)/(Fcil+Fciv);
    dhci_dFcil= Fciv*(xcil'*hcil-xciv'*hciv)/(Fcil+Fciv)^2;
    dhci_dFciv=-Fcil*(xcil'*hcil-xciv'*hciv)/(Fcil+Fciv)^2;
    dhci_dxcil= Fcil*hcil/(Fcil+Fciv);
    dhci_dxciv= Fciv*hciv/(Fcil+Fciv);
    dhci_dTcil= Fcil*xcil'*dhcil_dT/(Fcil+Fciv);
    dhci_dTciv= Fciv*xciv'*dhciv_dT/(Fcil+Fciv);

    hco       = xcov'*hcov; 
    dhco_dxcov= hcov;   
    dhco_dTcov= xcov'*dhcov_dT;
    
    hhi       = xhiv'*hhiv; 
    dhhi_dxhiv= hhiv;   
    dhhi_dThiv= xhiv'*dhhiv_dT;

    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  = (Fcil+Fciv)*(hco-hci)/(Tcov-Tciv);
    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_dFcil= ((hco-hci)-(Fcil+Fciv)*dhci_dFcil)/(Tcov-Tciv);
       dCmin_dFciv= ((hco-hci)-(Fcil+Fciv)*dhci_dFciv)/(Tcov-Tciv);
       
       dCmin_dxcil=-(Fcil+Fciv)*dhci_dxcil'/(Tcov-Tciv);
       dCmin_dxciv=-(Fcil+Fciv)*dhci_dxciv'/(Tcov-Tciv);
       dCmin_dxcov= (Fcil+Fciv)*dhco_dxcov'/(Tcov-Tciv);

       dCmin_dTcil=-(Fcil+Fciv)*dhci_dTcil/(Tcov-Tciv);
%      dCmin_dTciv=-(Fcil+Fciv)*dhci_dTciv/(Tcov-Tciv);
       dCmin_dTciv=-(Fcil+Fciv)*(dhci_dTciv/(Tcov-Tciv)-(hco-hci)/(Tcov-Tciv )^2);
       dCmin_dTcov= (Fcil+Fciv)*(dhco_dTcov/(Tcov-Tciv)-(hco-hci)/(Tcov-Tciv )^2);
            
       dCmin_dFhiv= 0; 
       dCmin_dFhov= 0; 
       dCmin_dFhol= 0; 
       dCmin_dxhiv= I*0;                   
       dCmin_dxhov= I*0;  
       dCmin_dxhol= I*0;  
       dCmin_dThiv = 0;
       dCmin_dThov = 0;
       dCmin_dThol = 0;
       
       
       dCmax_dFcil= 0;
       dCmax_dFciv= 0;
       dCmax_dxcil= I*0;
       dCmax_dxciv= I*0;
       dCmax_dxcov= I*0;
       dCmax_dTcil= 0;
       dCmax_dTciv= 0;
       dCmax_dTcov= 0;

    
       dCmax_dFhiv= (hho-hhi)/(Thov-Thiv);
       dCmax_dFhov= Fhiv*dhho_dFhov/(Thov-Thiv);
       dCmax_dFhol= Fhiv*dhho_dFhol/(Thov-Thiv);
       
       dCmax_dxhiv=-Fhiv*dhhi_dxhiv'/(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_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_dFcil= ((hco-hci)-(Fcil+Fciv)*dhci_dFcil)/(Tcov-Tciv);
       dCmax_dFciv= ((hco-hci)-(Fcil+Fciv)*dhci_dFciv)/(Tcov-Tciv);
       
       dCmax_dxcil=-(Fcil+Fciv)*dhci_dxcil'/(Tcov-Tciv);
       dCmax_dxciv=-(Fcil+Fciv)*dhci_dxciv'/(Tcov-Tciv);
       dCmax_dxcov= (Fcil+Fciv)*dhco_dxcov'/(Tcov-Tciv);

       dCmax_dTcil=-(Fcil+Fciv)*dhci_dTcil/(Tcov-Tciv);
%      dCmax_dTciv=-(Fcil+Fciv)*dhci_dTciv/(Tcov-Tciv);
       dCmax_dTciv=-(Fcil+Fciv)*(dhci_dTciv/(Tcov-Tciv)-(hco-hci)/(Tcov-Tciv )^2);
       dCmax_dTcov= (Fcil+Fciv)*(dhco_dTcov/(Tcov-Tciv)-(hco-hci)/(Tcov-Tciv )^2);

       dCmax_dFhiv= 0; 
       dCmax_dFhov= 0; 
       dCmax_dFhol= 0; 
       dCmax_dxhiv= I*0;                   
       dCmax_dxhov= I*0;  
       dCmax_dxhol= I*0;  
       dCmax_dThiv = 0;
       dCmax_dThov = 0;
       dCmax_dThol = 0;
       
       dCmin_dFcil= 0;
       dCmin_dFciv= 0;
       dCmin_dxcil= I*0;
       dCmin_dxciv= I*0;
       dCmin_dxcov= I*0;
       dCmin_dTcil= 0;
       dCmin_dTciv= 0;
       dCmin_dTcov= 0;

       dCmin_dFhiv= (hho-hhi)/(Thov-Thiv);
       dCmin_dFhov= Fhiv*dhho_dFhov/(Thov-Thiv);
       dCmin_dFhol= Fhiv*dhho_dFhol/(Thov-Thiv);
       
       dCmin_dxhiv=-Fhiv*dhhi_dxhiv'/(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_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;

    % Component balance cold side
    J(1:NC,Ixxcil)= Fcil*II;         J(1:NC,Ixxciv)= Fciv*II;           J(1:NC,Ixxcov)=-Fcov*II; 
    J(1:NC,IxFcil)= xcil;            J(1:NC,IxFciv)= xciv;              J(1:NC,IxFcov)=-xcov;  
    % Sum of components equal to one
    J(NC+1,Ixxcov)= I;
    % Energy balance cold side
    J(NC+2,IxFcil) = xcil'*hcil;           J(NC+2,IxFciv) = xciv'*hciv;    J(NC+2,IxFcov) =-xcov'*hcov;  
    J(NC+2,Ixxcil) = Fcil*hcil';           J(NC+2,Ixxciv) = Fciv*hciv';    J(NC+2,Ixxcov) =-Fcov*hcov';  
    J(NC+2,IxTcil) = Fcil*xcil'*dhcil_dT;  J(NC+2,IxTciv) = Fciv*xciv'*dhciv_dT; 
    J(NC+2,IxTcov) =-Fcov*xcov'*dhcov_dT;
    J(NC+2,IxQ)    = 1;
    %Pressure drop
    J(NC+3,IxPciv) = 1;                 J(NC+3,IxPcov) =-1;
    J(NC+4,IxPcil) = 1;                 J(NC+4,IxPcov) =-1;

    %Heat transfeer

    J(NC+5,IxQ)   =- 1;

    J(NC+5,Ixxcil) = ((de_dNtu*dNtu_dCmin*dCmin_dxcil + de_dRc*(dRc_dCmin*dCmin_dxcil + dRc_dCmax*dCmax_dxcil))*Cmin*(Thiv-Tciv)+e*dCmin_dxcil*(Thiv-Tciv)); 
    J(NC+5,Ixxciv) = ((de_dNtu*dNtu_dCmin*dCmin_dxciv + de_dRc*(dRc_dCmin*dCmin_dxciv + dRc_dCmax*dCmax_dxciv))*Cmin*(Thiv-Tciv)+e*dCmin_dxciv*(Thiv-Tciv)); 
    J(NC+5,Ixxcov) = ((de_dNtu*dNtu_dCmin*dCmin_dxcov + de_dRc*(dRc_dCmin*dCmin_dxcov + dRc_dCmax*dCmax_dxcov))*Cmin*(Thiv-Tciv)+e*dCmin_dxcov*(Thiv-Tciv));  
    J(NC+5,IxFcil) =  (de_dNtu*dNtu_dCmin*dCmin_dFcil + de_dRc*(dRc_dCmin*dCmin_dFcil + dRc_dCmax*dCmax_dFcil))*Cmin*(Thiv-Tciv)+e*dCmin_dFcil*(Thiv-Tciv); 
    J(NC+5,IxFciv) =  (de_dNtu*dNtu_dCmin*dCmin_dFciv + de_dRc*(dRc_dCmin*dCmin_dFciv + dRc_dCmax*dCmax_dFciv))*Cmin*(Thiv-Tciv)+e*dCmin_dFciv*(Thiv-Tciv); 
    J(NC+5,IxTcil) =  (de_dNtu*dNtu_dCmin*dCmin_dTcil + de_dRc*(dRc_dCmin*dCmin_dTcil + dRc_dCmax*dCmax_dTcil))*Cmin*(Thiv-Tciv)+e*dCmin_dTcil*(Thiv-Tciv);
    J(NC+5,IxTciv) =  (de_dNtu*dNtu_dCmin*dCmin_dTciv + de_dRc*(dRc_dCmin*dCmin_dTciv + dRc_dCmax*dCmax_dTciv))*Cmin*(Thiv-Tciv)+e*dCmin_dTciv*(Thiv-Tciv)- e*Cmin;
    J(NC+5,IxTcov) =  (de_dNtu*dNtu_dCmin*dCmin_dTcov + de_dRc*(dRc_dCmin*dCmin_dTcov + dRc_dCmax*dCmax_dTcov))*Cmin*(Thiv-Tciv)+e*dCmin_dTcov*(Thiv-Tciv);

    J(NC+5,Ixxhiv) = ((de_dNtu*dNtu_dCmin*dCmin_dxhiv + de_dRc*(dRc_dCmin*dCmin_dxhiv + dRc_dCmax*dCmax_dxhiv))*Cmin*(Thiv-Tciv)+e*dCmin_dxhiv*(Thiv-Tciv)); 
    J(NC+5,Ixxhov) = ((de_dNtu*dNtu_dCmin*dCmin_dxhov + de_dRc*(dRc_dCmin*dCmin_dxhov + dRc_dCmax*dCmax_dxhov))*Cmin*(Thiv-Tciv)+e*dCmin_dxhov*(Thiv-Tciv));  
    J(NC+5,Ixxhol) = ((de_dNtu*dNtu_dCmin*dCmin_dxhol + de_dRc*(dRc_dCmin*dCmin_dxhol + dRc_dCmax*dCmax_dxhol))*Cmin*(Thiv-Tciv)+e*dCmin_dxhol*(Thiv-Tciv));  
    J(NC+5,IxFhiv) =  (de_dNtu*dNtu_dCmin*dCmin_dFhiv + de_dRc*(dRc_dCmin*dCmin_dFhiv + dRc_dCmax*dCmax_dFhiv))*Cmin*(Thiv-Tciv)+e*dCmin_dFhiv*(Thiv-Tciv); 
    J(NC+5,IxFhov) =  (de_dNtu*dNtu_dCmin*dCmin_dFhov + de_dRc*(dRc_dCmin*dCmin_dFhov + dRc_dCmax*dCmax_dFhov))*Cmin*(Thiv-Tciv)+e*dCmin_dFhov*(Thiv-Tciv); 
    J(NC+5,IxFhol) =  (de_dNtu*dNtu_dCmin*dCmin_dFhol + de_dRc*(dRc_dCmin*dCmin_dFhol + dRc_dCmax*dCmax_dFhol))*Cmin*(Thiv-Tciv)+e*dCmin_dFhol*(Thiv-Tciv); 
    J(NC+5,IxThiv) =  (de_dNtu*dNtu_dCmin*dCmin_dThiv + de_dRc*(dRc_dCmin*dCmin_dThiv + dRc_dCmax*dCmax_dThiv))*Cmin*(Thiv-Tciv)+e*dCmin_dThiv*(Thiv-Tciv) + e*Cmin;
    J(NC+5,IxThol) =  (de_dNtu*dNtu_dCmin*dCmin_dThol + de_dRc*(dRc_dCmin*dCmin_dThol + dRc_dCmax*dCmax_dThol))*Cmin*(Thiv-Tciv)+e*dCmin_dThol*(Thiv-Tciv);
    J(NC+5,IxThov) =  (de_dNtu*dNtu_dCmin*dCmin_dThov + de_dRc*(dRc_dCmin*dCmin_dThov + dRc_dCmax*dCmax_dThov))*Cmin*(Thiv-Tciv)+e*dCmin_dThov*(Thiv-Tciv);
    J(NC+5,IxU)    =   de_dNtu*dNtu_dU*Cmin*(Thiv-Tciv);

    % -----------Flash-----------------------
     % Component balance
    J(NC+6:2*NC+5,Ixxhiv) = Fhiv*II;        J(NC+6:2*NC+5,Ixxhov) =-Fhov*II;          J(NC+6:2*NC+5,Ixxhol) =-Fhol*II; 
    J(NC+6:2*NC+5,IxFhiv) = xhiv;           J(NC+6:2*NC+5,IxFhov) =-xhov;             J(NC+6:2*NC+5,IxFhol) =-xhol;  
    % Sum of components equal to one
    J(2*NC+6,Ixxhov)= I;                  
    J(2*NC+7,Ixxhol)= I;
    % Energy balance
    J(2*NC+8,IxFhiv) = xhiv'*hhiv;          J(2*NC+8,IxFhov) =-xhov'*hhov;            J(2*NC+8,IxFhol) =-xhol'*hhol;
    J(2*NC+8,Ixxhiv) = Fhiv*hhiv';          J(2*NC+8,Ixxhov) =-Fhov*hhov';            J(2*NC+8,Ixxhol) =-Fhol*hhol';
    J(2*NC+8,IxThiv) = Fhiv*xhiv'*dhhiv_dT; J(2*NC+8,IxThov) =-Fhov*xhov'*dhhov_dT;   J(2*NC+8,IxThol) =-Fhol*xhol'*dhhol_dT;
    J(2*NC+8,IxQ) =-1;
    % Equlibrium 
    J(2*NC+9:3*NC+8,Ixxhov) = II;                    J(2*NC+9:3*NC+8,Ixxhol) =-K(Thov,Phov); 
    J(2*NC+9:3*NC+8,IxThov) =-dK_dT(Thov,Phov)*xhol; J(2*NC+9:3*NC+8,IxPhov) =-dK_dP(Thov,Phov)*xhol;
    % Outlet stream temperature
    J(3*NC+9,IxThov) = 1;   J(3*NC+9,IxThol) =-1;
    % Outlet stream pressure
    J(3*NC+10,IxPhiv) = 1;  J(3*NC+10,IxPhov) =-1;

    r = S*J;

else
    % Compute dxdt
    r=spalloc(Neq,1,Neq);
    % Compute enthalpies
    hciv =h(Tciv,1);
    hcil =h(Tcil,0);
    hcov =h(Tcov,1);
    hhiv =h(Thiv,1);
    hhov =h(Thov,1);
    hhol =h(Thol,0);
    
    hci       = (Fcil*xcil'*hcil+Fciv*xciv'*hciv)/(Fcil+Fciv);
    hco       = xcov'*hcov; 
    hhi       = xhiv'*hhiv; 
    hho       = (Fhol*xhol'*hhol+Fhov*xhov'*hhov)/(Fhol+Fhov);
    
    % e-Ntu model
    Cc  = (Fcil+Fciv)*(hco-hci)/(Tcov-Tciv );
    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
        Fciv*xciv+Fcil*xcil-Fcov*xcov                       % [kmol/s]  Component balance
        I*xcov-1                                            % [molfrac] Sum of components
        Fciv*xciv'*hciv+Fcil*xcil'*hcil-Fcov*xcov'*hcov+Q   % [kJ/s]    Energy balance
        Pciv-Pcov                                           % [bar]     Pressure drop vapour
        Pcil-Pcov                                           % [bar]     Pressure drop liquid
        e*Cmin*(Thiv-Tciv)-Q                                % [kJ/s]    Heat transfeer
        % Flash equations
        Fhiv*xhiv-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-Fhov*xhov'*hhov-Fhol*xhol'*hhol-Q   % [kmol/s]  Energy balance
        xhov-K(Thov,Phov)*xhol                              % [molfrac] Equilibrium 
        Thov-Thol                                           % [C]       Outlet stream temperature
        Phiv-Phov                                           % [kmol/s]  Feed stream pressure
       ];               

    r = S*r;
    
end