Merge pull request #611 from bobmyhill/p_eos

fix BM4, MT pressure calcs

Author: bobmyhill

burnman/calibrants/Decker_1971.py

@@ -10,7 +10,7 @@
 
 from burnman.eos.mie_grueneisen_debye import MGDBase
 from burnman.classes.calibrant import Calibrant
-from burnman.eos.birch_murnaghan_4th import birch_murnaghan_fourth
+from burnman.eos.birch_murnaghan_4th import pressure_birch_murnaghan_fourth
 
 
 class NaCl_B1(Calibrant):
@@ -27,7 +27,7 @@ class NaCl_B1(Calibrant):
     def __init__(self):
         def _pressure_Decker_NaCl(volume, temperature, params):
             # Isothermal pressure (GPa)
-            P0 = birch_murnaghan_fourth(params["V_0"] / volume, params)
+            P0 = pressure_birch_murnaghan_fourth(params["V_0"] / volume, params)
 
             # Thermal pressure
             thermal_model = MGDBase()

burnman/eos/birch_murnaghan_4th.py

@@ -19,8 +19,8 @@ def bulk_modulus_fourth(volume, params):
     modulus.  Pressure must be in :math:`[Pa]`.
     """
 
-    x = params["V_0"] / volume
-    f = 0.5 * (pow(x, 2.0 / 3.0) - 1.0)
+    invVrel = params["V_0"] / volume
+    f = 0.5 * (pow(invVrel, 2.0 / 3.0) - 1.0)
 
     Xi = (3.0 / 4.0) * (4.0 - params["Kprime_0"])
     Zeta = (3.0 / 8.0) * (
@@ -45,24 +45,46 @@ def bulk_modulus_fourth(volume, params):
 
 
 def volume_fourth_order(pressure, params):
-    func = lambda x: birch_murnaghan_fourth(params["V_0"] / x, params) - pressure
+    """
+    Volume according to the fourth-order
+    Birch-Murnaghan equation of state.
+
+    :param pressure: Pressure in the same units that are supplied for the reference bulk
+    modulus (params['K_0'])
+    :type pressure: float
+    :param params: parameter dictionary
+    :type params: dictionary
+    :return: V/V_0
+    :rtype: float
+    """
+
+    def delta_pressure(x):
+        return pressure_birch_murnaghan_fourth(params["V_0"] / x, params) - pressure
+
     try:
-        sol = bracket(func, params["V_0"], 1.0e-2 * params["V_0"])
-    except:
+        sol = bracket(delta_pressure, params["V_0"], 1.0e-2 * params["V_0"])
+    except ValueError:
         raise ValueError(
             "Cannot find a volume, perhaps you are outside of the range of validity for the equation of state?"
         )
-    return opt.brentq(func, sol[0], sol[1])
+    return opt.brentq(delta_pressure, sol[0], sol[1])
 
 
-def birch_murnaghan_fourth(x, params):
+def pressure_birch_murnaghan_fourth(invVrel, params):
     """
-    equation for the fourth order birch-murnaghan equation of state, returns
-    pressure in the same units that are supplied for the reference bulk
+    Pressure according to the fourth-order
+    Birch-Murnaghan equation of state.
+
+    :param invVrel: V/V_0
+    :type invVrel: float or numpy array
+    :param params: parameter dictionary
+    :type params: dictionary
+    :return: Pressure in the same units that are supplied for the reference bulk
     modulus (params['K_0'])
+    :rtype: float or numpy array
     """
 
-    f = 0.5 * (pow(x, 2.0 / 3.0) - 1.0)
+    f = 0.5 * (pow(invVrel, 2.0 / 3.0) - 1.0)
     Xi = (3.0 / 4.0) * (4.0 - params["Kprime_0"])
     Zeta = (3.0 / 8.0) * (
         (params["K_0"] * params["Kprime_prime_0"])
@@ -93,7 +115,7 @@ def volume(self, pressure, temperature, params):
         return volume_fourth_order(pressure, params)
 
     def pressure(self, temperature, volume, params):
-        return birch_murnaghan_fourth(volume / params["V_0"], params)
+        return pressure_birch_murnaghan_fourth(params["V_0"] / volume, params)
 
     def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
         """

burnman/eos/modified_tait.py

@@ -30,15 +30,21 @@ def tait_constants(params):
     return a, b, c
 
 
-def modified_tait(x, params):
+def pressure_modified_tait(Vrel, params):
     """
-    equation for the modified Tait equation of state, returns
-    pressure in the same units that are supplied for the reference bulk
-    modulus (params['K_0'])
+    Pressure according to the modified Tait equation of state.
     EQ 2 from Holland and Powell, 2011
+
+    :param Vrel: V/V_0
+    :type Vrel: float or numpy array
+    :param params: Parameter dictionary
+    :type params: dictionary
+    :return: pressure in the same units that are supplied for the reference bulk
+    modulus (params['K_0'])
+    :rtype: float or numpy array
     """
     a, b, c = tait_constants(params)
-    return (np.power((x + a - 1.0) / a, -1.0 / c) - 1.0) / b + params["P_0"]
+    return (np.power((Vrel + a - 1.0) / a, -1.0 / c) - 1.0) / b + params["P_0"]
 
 
 def volume(pressure, params):
@@ -47,8 +53,8 @@ def volume(pressure, params):
     EQ 12
     """
     a, b, c = tait_constants(params)
-    x = 1 - a * (1.0 - np.power((1.0 + b * (pressure - params["P_0"])), -1.0 * c))
-    return x * params["V_0"]
+    Vrel = 1.0 - a * (1.0 - np.power((1.0 + b * (pressure - params["P_0"])), -1.0 * c))
+    return Vrel * params["V_0"]
 
 
 def bulk_modulus(pressure, params):
@@ -113,7 +119,7 @@ def pressure(self, temperature, volume, params):
         """
         Returns pressure [Pa] as a function of temperature [K] and volume[m^3]
         """
-        return modified_tait(params["V_0"] / volume, params)
+        return pressure_modified_tait(volume / params["V_0"], params)
 
     def isothermal_bulk_modulus_reuss(self, pressure, temperature, volume, params):
         """