diff --git a/burnman/calibrants/Decker_1971.py b/burnman/calibrants/Decker_1971.py index 5a2dddef..1a64ccc7 100644 --- a/burnman/calibrants/Decker_1971.py +++ b/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() diff --git a/burnman/eos/birch_murnaghan_4th.py b/burnman/eos/birch_murnaghan_4th.py index 93ee5b4a..8ba753e2 100644 --- a/burnman/eos/birch_murnaghan_4th.py +++ b/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): """ diff --git a/burnman/eos/modified_tait.py b/burnman/eos/modified_tait.py index 90d176b2..f2bc5e2a 100644 --- a/burnman/eos/modified_tait.py +++ b/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): """