Skip to content

Commit

Permalink
Add two functions for DIPPR fitting
Browse files Browse the repository at this point in the history
  • Loading branch information
@CalebBell
CalebBell committed Nov 28, 2021
1 parent bc55f00 commit 9f93262
Show file tree
Hide file tree
Showing 2 changed files with 165 additions and 1 deletion.
151 changes: 150 additions & 1 deletion chemicals/dippr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@
'EQ114', 'EQ115', 'EQ116', 'EQ127',
'EQ101_fitting_jacobian', 'EQ102_fitting_jacobian',
'EQ106_fitting_jacobian', 'EQ105_fitting_jacobian',
'EQ107_fitting_jacobian']
'EQ107_fitting_jacobian',
'EQ106_AB', 'EQ106_ABC']

from chemicals.utils import log, exp, sinh, cosh, atan, atanh, sqrt, tanh
from cmath import log as clog
Expand Down Expand Up @@ -837,6 +838,154 @@ def EQ106(T, Tc, A, B, C=0.0, D=0.0, E=0.0, order=0):
else:
raise ValueError(order_not_found_msg)

def EQ106_AB(T, Tc, val, der):
r'''Calculate the coefficients `A` and `B` of the DIPPR Equation #106 using
the value of the function and its derivative at a specific point.
.. math::
A = val \left(\frac{1}{Tc} \left(- T + Tc\right)\right)^{- \frac{der}{val} \left(T - Tc\right)}
.. math::
B = \frac{der}{val} \left(T - Tc\right)
Parameters
----------
T : float
Temperature, [K]
Tc : float
Critical temperature, [K]
val : float
Property value [constant-specific]
der : float
First temperature derivative of property value [constant-specific/K]
Returns
-------
A : float
Parameter for the equation [constant-specific]
B : float
Parameter for the equation [-]
Notes
-----
Examples
--------
>>> val = EQ106(300, 647.096, A=0.17766, B=2.567)
>>> der = EQ106(300, 647.096, A=0.17766, B=2.567, order=1)
>>> EQ106_AB(300, 647.096, val, der)
(0.17766, 2.567)
'''

'''# Derived with:
from sympy import *
T, Tc, A, B, val, der = symbols('T, Tc, A, B, val, der')
Tr = T/Tc
expr = A*(1 - Tr)**B
Eq0 = Eq(expr, val)
Eq1 = Eq(diff(expr, T), der)
s = solve([Eq0, Eq1], [A, B])
'''
x0 = T - Tc
x1 = der*x0/val
A, B = val*(-x0/Tc)**(-x1), x1
return (A, B)

def EQ106_ABC(T, Tc, val, der, der2):
r'''Calculate the coefficients `A`, `B`, and `C` of the DIPPR Equation #106
using, the value of the function and its first and second derivative at a
specific point.
.. math::
A = val \left(\frac{1}{Tc} \left(- T + Tc\right)\right)^{\frac{1}{val^{2}
\left(\log{\left (\frac{1}{Tc} \left(- T + Tc\right) \right )} + 2\right)}
\left(T \left(\log{\left (\frac{1}{Tc} \left(- T + Tc\right) \right )}
+ 1\right) \left(- T der^{2} + T der_{2} val + Tc der^{2} - Tc der_{2}
val + der val\right) - T \left(- T der^{2} + T der_{2} val + Tc der^{2}
- Tc der_{2} val + der val\right) - Tc \left(- T der^{2} + T der_{2} val
+ Tc der^{2} - Tc der_{2} val + der val\right) \log{\left (\frac{1}{Tc}
\left(- T + Tc\right) \right )} - der val \left(T - Tc\right)
\left(\log{\left (\frac{1}{Tc} \left(- T + Tc\right) \right )}
+ 2\right)\right)}
.. math::
B = \frac{1}{val^{2} \left(\log{\left (\frac{1}{Tc} \left(- T + Tc\right)
\right )} + 2\right)} \left(- T \left(\log{\left (\frac{1}{Tc}
\left(- T + Tc\right) \right )} + 1\right) \left(- T der^{2} + T der_{2}
val + Tc der^{2} - Tc der_{2} val + der val\right) + Tc \left(- T der^{2}
+ T der_{2} val + Tc der^{2} - Tc der_{2} val + der val\right)
\log{\left (\frac{1}{Tc} \left(- T + Tc\right) \right )} + der val
\left(T - Tc\right) \left(\log{\left (\frac{1}{Tc} \left(- T + Tc\right) \right )} + 2\right)\right)
.. math::
C = \frac{Tc \left(- T der^{2} + T der_{2} val + Tc der^{2} - Tc der_{2}
val + der val\right)}{val^{2} \left(\log{\left (\frac{1}{Tc}
\left(- T + Tc\right) \right )} + 2\right)}
Parameters
----------
T : float
Temperature, [K]
Tc : float
Critical temperature, [K]
val : float
Property value [constant-specific]
der : float
First temperature derivative of property value [constant-specific/K]
der2 : float
Second temperature derivative of property value [constant-specific/K^2]
Returns
-------
A : float
Parameter for the equation [constant-specific]
B : float
Parameter for the equation [-]
C : float
Parameter for the equation [-]
Notes
-----
Examples
--------
>>> val = EQ106(300, 647.096, A=0.17766, B=2.567, C=-0.01)
>>> der = EQ106(300, 647.096, A=0.17766, B=2.567, C=-0.01, order=1)
>>> der2 = EQ106(300, 647.096, A=0.17766, B=2.567, C=-0.01, order=2)
>>> EQ106_ABC(300, 647.096, val, der, der2)
(0.17766, 2.567, -0.01)
'''
'''# Broken in recent versions of SymPy, SymPy 1.1 is good
from sympy import *
T, Tc, A, B, C, val, der, der2 = symbols('T, Tc, A, B, C, val, der, der2')
Tr = T/Tc
expr = A*(1 - Tr)**(B + C*Tr)
Eq0 = Eq(expr, val)
Eq1 = Eq(diff(expr, T), der)
Eq2 = Eq(diff(expr, T, 2), der2)
s = solve([Eq0, Eq1, Eq2], [A, B, C])
'''
x0 = T - Tc
x1 = -x0/Tc
x2 = log(x1)
x3 = x2 + 2
x4 = 1/(val*val*x3)
x5 = der*val
x6 = der2*val
x7 = der*der
x8 = T*x6 - T*x7 - Tc*x6 + Tc*x7 + x5
x9 = T*x8
x10 = Tc*x8
x11 = x0*x3*x5 + x10*x2 - x9*(x2 + 1)
A, B, C = val*x1**(-x4*(x11 + x9)), x11*x4, x10*x4
return (A, B, C)


def EQ107(T, A=0, B=0, C=0, D=0, E=0, order=0):
r'''DIPPR Equation #107. Often used in calculating ideal-gas heat capacity.
Expand Down
15 changes: 15 additions & 0 deletions tests/test_dippr.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,21 @@ def test_EQ106_more():
assert 0.0 == EQ106(647.097, 647.096, 0.17766, 2.567, -3.3377, 1.9699)


def test_EQ106_AB_and_ABC():
Tmin, Tc, A, B = 194.0, 592.5, 0.056, 1.32
C = -0.01
Tmax = 590.0

kwargs = {'T': 590.0, 'Tc': 592.5, 'val': 4.106957515154657e-05, 'der': -2.1684735680016856e-05}

res = EQ106_AB(**kwargs)
assert_close1d(res, (0.056, 1.32), rtol=1e-13)

kwargs = {'T': 590.0, 'Tc': 592.5, 'val': 4.33678101801987e-05, 'der': -2.2721462154946353e-05, 'der2': 2.814732501661309e-06}
res = EQ106_ABC(**kwargs)
assert_close1d(res, (0.056, 1.32, -0.01), rtol=1e-9)


def test_EQ101_more():
assert_close(EQ101(300.0, 73.649, -7258.2, -7.3037, 4.1653E-6, 2, order=1), 208.00259945348506, rtol=1e-13)
assert_close(derivative(lambda T: EQ101(T, 73.649, -7258.2, -7.3037, 4.1653E-6, 2), 300, dx=1e-4),
Expand Down

0 comments on commit 9f93262

Please sign in to comment.