Complete three more jacobians of temperature dependent fit equations

chemicals/dippr.py

@@ -50,15 +50,19 @@
 -----------------------
 .. autofunction:: chemicals.dippr.EQ101_fitting_jacobian
 .. autofunction:: chemicals.dippr.EQ102_fitting_jacobian
-
+.. autofunction:: chemicals.dippr.EQ105_fitting_jacobian
+.. autofunction:: chemicals.dippr.EQ106_fitting_jacobian
+.. autofunction:: chemicals.dippr.EQ107_fitting_jacobian
 
 """
 
 from __future__ import division
 
 __all__ = ['EQ100', 'EQ101', 'EQ102', 'EQ104', 'EQ105', 'EQ106', 'EQ107',
            'EQ114', 'EQ115', 'EQ116', 'EQ127', 
-           'EQ101_fitting_jacobian', 'EQ102_fitting_jacobian']
+           'EQ101_fitting_jacobian', 'EQ102_fitting_jacobian',
+           'EQ106_fitting_jacobian', 'EQ105_fitting_jacobian',
+           'EQ107_fitting_jacobian']
 
 from chemicals.utils import log, exp, sinh, cosh, atan, atanh, sqrt, tanh
 from cmath import log as clog
@@ -404,6 +408,125 @@ def EQ102_fitting_jacobian(Ts, A, B, C, D):
         out[i][3] = -x2*x5
     return out
 
+def EQ105_fitting_jacobian(Ts, A, B, C, D):
+    r'''Compute and return the Jacobian of the property predicted by 
+    DIPPR Equation # 105 with respect to all the coefficients. This is used in
+    fitting parameters for chemicals.
+
+    Parameters
+    ----------
+    Ts : list[float]
+        Temperatures of the experimental data points, [K]
+    A-D : float
+        Parameter for the equation; chemical and property specific [-]
+
+    Returns
+    -------
+    jac : list[list[float, 4], len(Ts)]
+        Matrix of derivatives of the equation with respect to the fitting
+        parameters, [various]
+
+    '''
+    N = len(Ts)
+#    out = np.zeros((N, 4)) # numba: uncomment
+    out = [[0.0]*4 for _ in range(N)] # numba: delete
+    for i in range(N):
+        r = out[i]
+        x0 = 1.0 - Ts[i]/C
+
+        if D < 1.0 and x0 < 0.0:
+            r[0] = 1.0/B
+            r[1] = -A/(B*B)
+        else:
+            x1 = x0**D
+            x2 = x1 + 1.0
+            x3 = A*B**(-x1 - 1.0)
+            x4 = x1*x3*log(B)
+            r[0] = B**(-x2)
+            r[1] = -x2*x3/B
+            r[2] = -D*Ts[i]*x4/(C*C*x0)
+            r[3] = -x4*log(x0)
+    return out
+
+def EQ106_fitting_jacobian(Ts, Tc, A, B, C, D, E):
+    r'''Compute and return the Jacobian of the property predicted by 
+    DIPPR Equation # 106 with respect to all the coefficients. This is used in
+    fitting parameters for chemicals.
+
+    Parameters
+    ----------
+    Ts : list[float]
+        Temperatures of the experimental data points, [K]
+    Tc : float
+        Critical temperature, [K]
+    A-E : float
+        Parameter for the equation; chemical and property specific [-]
+
+    Returns
+    -------
+    jac : list[list[float, 5], len(Ts)]
+        Matrix of derivatives of the equation with respect to the fitting
+        parameters, [various]
+
+    '''
+    N = len(Ts)
+#    out = np.zeros((N, 5)) # numba: uncomment
+    out = [[0.0]*5 for _ in range(N)] # numba: delete
+    for i in range(N):
+        x0 = Ts[i]/Tc
+        if x0 != 1.0:
+            x1 = 1.0 - x0
+            x2 = x1**(B + x0*(C + x0*(D + E*x0)))
+            x3 = A*x2*log(x1)
+            r = out[i]
+            r[0] = x2
+            r[1] = x3
+            r[2] = x0*x3
+            r[3] = x0*x0*x3
+            r[4] = x0*x0*x0*x3
+    return out
+
+def EQ107_fitting_jacobian(Ts, A, B, C, D, E):
+    r'''Compute and return the Jacobian of the property predicted by 
+    DIPPR Equation # 107 with respect to all the coefficients. This is used in
+    fitting parameters for chemicals.
+
+    Parameters
+    ----------
+    Ts : list[float]
+        Temperatures of the experimental data points, [K]
+    A-E : float
+        Parameter for the equation; chemical and property specific [-]
+
+    Returns
+    -------
+    jac : list[list[float, 5], len(Ts)]
+        Matrix of derivatives of the equation with respect to the fitting
+        parameters, [various]
+
+    '''
+    N = len(Ts)
+#    out = np.zeros((N, 5)) # numba: uncomment
+    out = [[0.0]*5 for _ in range(N)] # numba: delete
+    for i in range(N):
+        r = out[i]
+        x1 = 1.0/Ts[i]
+        x0 = x1*x1
+        x2 = C*x1
+        x3 = sinh(x2)
+        x3_inv = 1.0/x3
+        x4 = x0*x3_inv*x3_inv
+        x5 = E*x1
+        x6 = cosh(x5)
+        x6_inv = 1.0/x6
+        x7 = x0*x6_inv*x6_inv
+        r[0] = 1.0
+        r[1] = C*C*x4
+        r[2] = 2.0*B*C*x4*(-x2*cosh(x2)*x3_inv + 1.0)
+        r[3] = E*E*x7
+        r[4] = 2.0*D*E*x7*(-x5*sinh(x5)*x6_inv + 1.0)
+    return out
+
 def EQ104(T, A, B, C=0.0, D=0.0, E=0.0, order=0):
     r'''DIPPR Equation #104. Often used in calculating second virial
     coefficients of gases. All 5 parameters are required.
@@ -788,7 +911,11 @@ def EQ107(T, A=0, B=0, C=0, D=0, E=0, order=0):
        doi:10.1016/0378-3812(81)85002-9.
     '''
     if order == 0:
-        return A + B*((C/T)/sinh(C/T))**2 + D*((E/T)/cosh(E/T))**2
+        C_T = C/T
+        t0 = 2.0*C_T/(trunc_exp(C_T) - trunc_exp(-C_T))
+        E_T = E/T
+        t1 = 2.0*E_T/(trunc_exp(-E_T) + trunc_exp(E_T))
+        return A + B*t0*t0 + D*t1*t1
     elif order == 1:
         return (2*B*C**3*cosh(C/T)/(T**4*sinh(C/T)**3)
                 - 2*B*C**2/(T**3*sinh(C/T)**2)

chemicals/numba.py

@@ -75,9 +75,13 @@ def transform_complete_chemicals(replaced, __funcs, __all__, normal, vec):
 
              'vapor_pressure.Wagner_fitting_jacobian',
              'vapor_pressure.Wagner_original_fitting_jacobian',
+             'vapor_pressure.Antoine_fitting_jacobian',
 
              'dippr.EQ102_fitting_jacobian',
              'dippr.EQ101_fitting_jacobian',
+             'dippr.EQ106_fitting_jacobian',
+             'dippr.EQ105_fitting_jacobian',
+             'dippr.EQ107_fitting_jacobian',
 
              'vapor_pressure.Yaws_Psat_fitting_jacobian',
              'viscosity.mu_Yaws_fitting_jacobian',

chemicals/vapor_pressure.py

@@ -57,6 +57,7 @@
 -----------------------
 .. autofunction:: chemicals.vapor_pressure.Wagner_fitting_jacobian
 .. autofunction:: chemicals.vapor_pressure.Wagner_original_fitting_jacobian
+.. autofunction:: chemicals.vapor_pressure.Antoine_fitting_jacobian
 .. autofunction:: chemicals.vapor_pressure.Yaws_Psat_fitting_jacobian
 
 Vapor Pressure Estimation Correlations
@@ -165,6 +166,7 @@
            'Antoine_coeffs_from_point', 'Antoine_AB_coeffs_from_point',
            'DIPPR101_ABC_coeffs_from_point', 'Wagner_original_fitting_jacobian',
            'Wagner_fitting_jacobian', 'Yaws_Psat_fitting_jacobian',
+           'Antoine_fitting_jacobian',
            'TDE_PVExpansion']
 
 import os
@@ -1327,8 +1329,8 @@ def Wagner_original_fitting_jacobian(Ts, Tc, Pc, a, b, c, d):
 
     Parameters
     ----------
-    T : float
-        Temperature of fluid, [K]
+    Ts : list[float]
+        Temperatures of fluid data points, [K]
     Tc : float
         Critical temperature, [K]
     Pc : float
@@ -1369,8 +1371,8 @@ def Wagner_fitting_jacobian(Ts, Tc, Pc, a, b, c, d):
 
     Parameters
     ----------
-    T : float
-        Temperature of fluid, [K]
+    Ts : list[float]
+        Temperatures of fluid data points, [K]
     Tc : float
         Critical temperature, [K]
     Pc : float
@@ -1402,6 +1404,50 @@ def Wagner_fitting_jacobian(Ts, Tc, Pc, a, b, c, d):
         row[3] = x4*x5
     return out
 
+def Antoine_fitting_jacobian(Ts, A, B, C, base=10.0):
+    r'''Calculates the jacobian of the Antoine vapor pressure equation
+    for use in fitting these parameters when experimental values are known.
+    
+    Requires three coefficients specific to each chemical.
+
+    Parameters
+    ----------
+    Ts : list[float]
+        Temperatures of fluid data points, [K]
+    A : float
+        Antoine `A` parameter, [-]
+    B : float
+        Antoine `B` parameter, [K]
+    C : float
+        Antoine `C` parameter, [K]
+    base : float, optional
+        Optional base of logarithm; 10 by default, [-]
+
+    Returns
+    -------
+    jac : list[list[float, 3], len(Ts)]
+        Matrix of derivatives of the equation with respect to the fitting
+        parameters, [various]
+    '''
+    N = len(Ts)
+#    out = np.zeros((N, 3)) # numba: uncomment
+    out = [[0.0]*3 for _ in range(N)] # numba: delete
+    ln_base = log(base)
+    for i in range(N):
+        row = out[i]
+        x0 = C + Ts[i]
+        if x0 <= 0.0:
+            row[0] = 0.0
+            row[1] = 0.0
+            row[2] = 0.0
+        else:
+            x1 = 1.0/x0
+            x2 = base**(A - B*x1)*ln_base
+            row[0] = x2
+            row[1] = -x1*x2
+            row[2] =  B*x2*x1*x1
+    return out
+
 def Wagner(T, Tc, Pc, a, b, c, d):
     r'''Calculates vapor pressure using the Wagner equation (2.5, 5 form).
 

tests/test_dippr.py

@@ -227,6 +227,10 @@ def test_EQ107_more():
 
     with pytest.raises(Exception):
         EQ107(20., *coeffs, order=1E100)
+
+    # Case that requires overflow handling
+    EQ107(**{'T': 377.77777777777777, 'A': 1539249.2020718465, 'B': -46807441.804555826, 'C': -409401.9169728528, 'D': -2164118.45731599, 'E': 339.5030595758336, 'order': 0})
+
 
 def test_EQ114_more():
     # T derivative
@@ -367,3 +371,56 @@ def test_EQ102_fitting():
     der_analytical = EQ102_fitting_jacobian([T], A, B, C, D)
     assert_close1d(der_analytical, [der_num])
     assert_close1d(der_analytical, der_expect, rtol=1e-13)
+
+def test_EQ105_fitting():
+    T, A, B, C, D = 300., 0.70824, 0.26411, 507.6, 0.27537
+    der_num = [derivative(lambda A: EQ105(T, A, B, C, D), A, dx=A*1e-5),
+                 derivative(lambda B: EQ105(T, A, B, C, D), B, dx=B*1e-5),
+                 derivative(lambda C: EQ105(T, A, B, C, D), C, dx=C*1e-5),
+                 derivative(lambda D: EQ105(T, A, B, C, D), D, dx=D*1e-5)]
+    der_expect = [[10.721182221195127, -51.225752844842916, 0.006195731841673147, -7.0661094492142285]]
+    der_analytical = EQ105_fitting_jacobian([T], A, B, C, D)
+    assert_close1d(der_analytical, [der_num])
+    assert_close1d(der_analytical, der_expect, rtol=1e-13)
+
+
+    T, A, B, C, D = 304, 3733.087734888731, 0.30552803535622014, 301.6993863907116, 0.3415888512743092
+    der_num = [derivative(lambda A: EQ105(T, A, B, C, D), A, dx=A*1e-5),
+                 derivative(lambda B: EQ105(T, A, B, C, D), B, dx=B*1e-5),
+                 derivative(lambda C: EQ105(T, A, B, C, D), C, dx=C*1e-5),
+                 derivative(lambda D: EQ105(T, A, B, C, D), D, dx=D*1e-5)]
+    der_analytical = EQ105_fitting_jacobian([T], A, B, C, D)
+    assert_close1d(der_analytical, [der_num])
+
+
+def test_EQ106_fitting():
+    T, Tc, A, B, C, D, E = 300, 647.096, 0.17766, 2.567, -3.3377, 1.9699, 0.25
+    der_num = [derivative(lambda A: EQ106(T, Tc, A, B, C, D, E), A, dx=A*1e-5),
+                 derivative(lambda B: EQ106(T, Tc, A, B, C, D, E), B, dx=B*1e-5),
+                 derivative(lambda C: EQ106(T, Tc, A, B, C, D, E), C, dx=C*1e-5),
+                 derivative(lambda D: EQ106(T, Tc, A, B, C, D, E), D, dx=D*1e-5),
+                 derivative(lambda E: EQ106(T, Tc, A, B, C, D, E), E, dx=E*1e-5),
+              ]
+
+    der_expect = [[0.4007741423076755, -0.04435095583995359, -0.020561534535812425, -0.009532527415937865, -0.004419372434354963]]
+    der_analytical = EQ106_fitting_jacobian([T], Tc, A, B, C, D, E)
+    assert_close1d(der_analytical, [der_num])
+    assert_close1d(der_analytical, der_expect, rtol=1e-13)
+
+    # Test case with errors
+    EQ106_fitting_jacobian(Ts=[466.],Tc = 466.0, A = 47700.0, B = 0.37, C = 0.0, D = 0.0, E = 0.0)
+
+def test_EQ107_fitting():
+    T, A, B, C, D, E = 250.0, 33363., 26790., 2610.5, 8896., 1169
+    der_num = [derivative(lambda A: EQ107(T, A, B, C, D, E), A, dx=A*1e-5),
+                 derivative(lambda B: EQ107(T, A, B, C, D, E), B, dx=B*3e-3, order=3),
+                 derivative(lambda C: EQ107(T, A, B, C, D, E), C, dx=C*1e-3, order=9),
+                 derivative(lambda D: EQ107(T, A, B, C, D, E), D, dx=D*1e-5, order=3),
+                 derivative(lambda E: EQ107(T, A, B, C, D, E), E, dx=E*1e-5),
+              ]
+
+    der_expect = [[1.0, 3.7138247806865474e-07, -7.197214038362036e-05, 0.00758947296962729, -0.42452325330497365]]
+    der_analytical = EQ107_fitting_jacobian([T], A, B, C, D, E)
+    assert_close1d(der_analytical, [der_num])
+    assert_close1d(der_analytical, der_expect, rtol=1e-13)
+

tests/test_vapor_pressure.py

@@ -378,6 +378,27 @@ def test_Wagner_fitting_jacobian():
     assert_close1d(der_expect, der_analytical, rtol=1e-13)
     assert_close1d(der_analytical, [der_num], rtol=1e-7)
 
+def test_Antoine_fitting_jacobian():
+
+    T, A, B, C = 100.0, 8.7687, 395.744, -6.469
+    der_num = [derivative(lambda A: Antoine(T, A, B, C), A, dx=A*1e-5),
+     derivative(lambda B: Antoine(T, A, B, C), B, dx=B*1e-5),
+     derivative(lambda C: Antoine(T, A, B, C), C, dx=C*1e-5)]
+    der_expect =  [[79389.37469005348, -848.802800034785, 3591.414774748115]]
+    der_analytical = Antoine_fitting_jacobian([T], A, B, C)
+    assert_close1d(der_analytical, [der_num], rtol=1e-7)
+    assert_close1d(der_expect, der_analytical, rtol=1e-13)
+
+    # Zero point
+    T, A, B, C = 30, 3.45604+5, 1044.038, -53.893
+    der_num = [derivative(lambda A: Antoine(T, A, B, C), A, dx=A*1e-5),
+     derivative(lambda B: Antoine(T, A, B, C), B, dx=B*1e-5),
+     derivative(lambda C: Antoine(T, A, B, C), C, dx=C*1e-5)]
+    der_expect =  [[0.0, 0.0, 0.0]]
+    der_analytical = Antoine_fitting_jacobian([T], A, B, C)
+    assert_close1d(der_analytical, [der_num], rtol=1e-7)
+    assert_close1d(der_expect, der_analytical, rtol=1e-13)
+
 def test_Yaws_Psat():
     assert_close(Yaws_Psat(T=400.0, A=28.588+ log10(101325/760), B=-2469, C=-7.351, D=2.8025E-10, E=2.7361E-6), 708657.0891069275, rtol=1e-14)