Fixed default method in example. Reformatted BM4 to f notation

CaymanUnterborn · Feb 24, 2016 · b47c22f · b47c22f
1 parent cb8e698
commit b47c22f
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 12 deletions.
diff --git a/burnman/eos/birch_murnaghan_4th.py b/burnman/eos/birch_murnaghan_4th.py
@@ -15,31 +15,39 @@ def bulk_modulus_fourth(volume, params):
     modulus in the same units as the reference bulk
     modulus.  Pressure must be in :math:`[Pa]`.
     """
-    B1 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-5.))+(59./9.)
-    B2 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-13.))+(129./9.)
-    B3 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-11.))+(105./9.)
-    B4 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-3.))+(35./9.)
 
     x = params['V_0']/volume
+    f = 0.5*(pow(x, 2./3.) - 1.0)
+
+    Xi = (3./4.)*(4.-params['Kprime_0'])
+    Zeta = (3./8.)*((params['K_0']*params['Kprime_prime_0'])+params['Kprime_0']*(params['Kprime_0']-7.)+143./9.)
+
+    K = (5.*f*pow((1.+2.*f),5./2.)*params['K_0']*(1.-(2.*Xi*f)+(4.*Zeta*pow(f,2.)))) + \
+        (pow(1.+(2.*f),7./2.)*params['K_0']*(1.-(4.*Xi*f)+(12.*Zeta*pow(f,2.))))
+
 
-    K = (9./16.)*params['K_0']*( (-5.*B1/3.)*pow(x,-5./3.)+(7.*B2/3.)*pow(x,-7./3.)-3*B3*pow(x,-3.)+(11.*B4/3.)*pow(x,-11./3.))
     return K
 
 def volume_fourth_order(pressure,params):
-    func = lambda x: birch_murnaghan_fourth(x/params['V_0'], params) - pressure
+    func = lambda x: birch_murnaghan_fourth(params['V_0']/x, params) - pressure
     try:
         sol = bracket(func, params['V_0'], 1.e-2*params['V_0'])
     except:
         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])
 
 def birch_murnaghan_fourth(x, params):
-    B1 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-5.))+(59./9.)
-    B2 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-13.))+(129./9.)
-    B3 = (3.*params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*((3.*params['Kprime_0'])-11.))+(105./9.)
-    B4 = (params['K_0']*params['Kprime_prime_0'])+((params['Kprime_0']-4.)*(params['Kprime_0']-3.))+(35./9.)
+    """
+    equation for the fourth order birch-murnaghan equation of state, returns
+    pressure in the same units that are supplied for the reference bulk
+    modulus (params['K_0'])
+    """
+
+    f = 0.5*(pow(x, 2./3.) - 1.0)
+    Xi = (3./4.)*(4.-params['Kprime_0'])
+    Zeta = (3./8.)*((params['K_0']*params['Kprime_prime_0'])+params['Kprime_0']*(params['Kprime_0']-7.)+143./9.)
 
-    return (9.*params['K_0']/16.)*((-B1*pow(x, -5./3.)+(B2*pow(x, -7./3.))) - (B3*pow(x, -3.))+(B4*pow(x, -11./3.)))
+    return 3.*f*pow(1.+2.*f,5./2.)*params['K_0']*(1.-(2.*Xi*f)+(4.*Zeta*pow(f,2.)))
 
 
 class BM4(eos.EquationOfState):

diff --git a/examples/example_user_input_material.py b/examples/example_user_input_material.py
@@ -59,7 +59,7 @@
         (your choice in geotherm will not matter in this case)
     or 'vinet' (vinet equation of state, if you choose to ignore temperature
         (your choice in geotherm will not matter in this case)))"""
-    method = 'vinet'
+    method = 'slb3'
 
     #in form name_of_mineral (burnman.mineral <- creates list with parameters)
     class own_material (burnman.Mineral):