Function hmez(T As Double, P As Double, Taby As Range, _
                 Propiedades As Range) As Double

    Dim cp, Ccp(), Ccpm(5), y(), h0(), hideal, tau, Mcp, tem, href, T0, R As Double
    Dim n, i, j As Integer
    
    n = Propiedades.Rows.Count
    
    ReDim Ccp(n, 5), y(n), h0(n)
    
    T0 = 298.15 'K
    
    For i = 1 To n
        y(i) = Taby(i).Value
        h0(i) = Propiedades(i, 11) 'kJ/mol
        For j = 1 To 5
            Ccp(i, j) = Propiedades(i, j).Value
        Next j
    Next i
    
    For j = 1 To 5
        Ccpm(j) = 0
        For i = 1 To n
            Ccpm(j) = Ccpm(j) + y(i) * Ccp(i, j)
        Next i
    Next j
    
    href = 0
    For i = 1 To n
        href = href + y(i) * h0(i)
    Next i
    

    tau = T / T0
    
    tem = Ccpm(1) _
         + Ccpm(2) * 10 ^ (-3) / 2 * T0 * (tau + 1) _
         + Ccpm(3) * 10 ^ (-5) / 3 * T0 * T0 * (tau * tau + tau + 1) _
         + Ccpm(4) * 10 ^ (-8) / 4 * T0 * T0 * T0 * (tau * tau * tau + tau * tau + tau + 1) _
         + Ccpm(5) * 10 ^ (-11) / 5 * T0 * T0 * T0 * T0 * (tau * tau * tau * tau + tau * tau * tau + tau * tau + tau + 1)
          
    Mcp = 8.3145 / 1000 * tem 'kJ/mol/K
    hideal = href + Mcp * (T - T0) 'kJ/mol

    
    'Entalpia residual
    
    Dim Tc(), Pc(), w(), vc(), zc()  As Double
    Dim Tcij, Pcij, wij, vcij, zcij, Trij As Double
    

    
    Dim Bo, B1, Bij, B As Double
    Dim dBodT, dB1dT, dBdTij, dBdT, hR As Double
           
                
    ReDim Tc(n), Pc(n), w(n), vc(n), zc(n)

    
    For i = 1 To n
        Tc(i) = Propiedades(i, 6).Value
        Pc(i) = Propiedades(i, 7).Value
        vc(i) = Propiedades(i, 8).Value
        zc(i) = Propiedades(i, 9).Value
        w(i) = Propiedades(i, 10).Value
    Next i
    
    '8.3145 J mol^-1 K^-1 = 8.3145 Pa m^3 mol^-1 K^-1
    B = 0
    dBdT = 0
    For i = 1 To n
        For j = 1 To n
            Tcij = Sqr(Tc(i) * Tc(j)) 'K
            vcij = ((vc(i) ^ (1 / 3) + vc(j) ^ (1 / 3)) / 2) ^ 3 * 10 ^ (-6) 'm^3/mol
            zcij = (zc(i) + zc(j)) / 2
            wij = (w(i) + w(j)) / 2
            Pcij = 8.3145 * Tcij * zcij / vcij ' Pa -> Pa m^3 mol^-1 K^-1 K mol m^-3
            Trij = T / Tcij
            Bo = 0.083 - 0.422 / Trij ^ 1.6
            B1 = 0.139 - 0.172 / Trij ^ 4.2
            Bij = 8.314 * Tcij / Pcij * (Bo + wij * B1) ' m^3/mol -> Pa m^3 mol^-1 K^-1 K Pa^-1
            B = B + y(i) * y(j) * Bij 'm^3/mol
            dBodT = 0.6752 / Trij ^ 2.6
            dB1dT = 0.7224 / Trij ^ 5.2
            dBdTij = 8.314 / Pcij * (dBodT + wij * dB1dT) 'm^3/mol/K
            dBdT = dBdT + y(i) * y(j) * dBdTij 'm^3/mol/K
        Next j
    Next i
          
    hR = (P * 10 ^ 5) * (B - T * dBdT) * (1 / 1000) ' kJ/mol -> Pa m^3 mol^-1 /1000 ->kJ mol^-1
    
    ' 1 bar cm^3/mol -> 100 kJ/mol
    'hR = hR * 100 'kJ/mol
     
    
    hmez = hideal + hR 'kJ/mol
    
    
End Function