Merge pull request #614 from bobmyhill/improve_spock

tweak spock eos

burnman/eos/spock.py

@@ -24,24 +24,35 @@ def jit(fn):
         return fn
 
 
-def gammaincc(a, x):
+def generalised_gammainc(a, x1, x2):
     """
-    An implementation of the non-regularised upper incomplete gamma
+    An implementation of the generalised incomplete gamma
     function. Computed using the relationship with the regularised
     lower incomplete gamma function (scipy.special.gammainc).
     Uses the recurrence relation wherever a<0.
+
+    We could have used mpmath.gammainc(a, x1, x2) directly,
+    but it is significantly slower than this implementation.
     """
     n = int(-np.floor(a))
     if n > 0:
         a = a + n
-        u_gamma = exp1(x) if a == 0 else (1.0 - gammainc(a, x)) * gamma(a)
+        u_gamma = (
+            exp1(x1) - exp1(x2)
+            if np.abs(a) < 1.0e-12
+            else (gammainc(a, x2) - gammainc(a, x1)) * gamma(a)
+        )
 
+        esubx1 = np.exp(-x1)
+        esubx2 = np.exp(-x2)
         for _ in range(n):
             a = a - 1.0
-            u_gamma = (u_gamma - np.power(x, a) * np.exp(-x)) / a
+            u_gamma = (
+                u_gamma + np.power(x2, a) * esubx2 - np.power(x1, a) * esubx1
+            ) / a
         return u_gamma
     else:
-        return (1.0 - gammainc(a, x)) * gamma(a)
+        return (gammainc(a, x2) - gammainc(a, x1)) * gamma(a)
 
 
 @jit(nopython=True)
@@ -98,13 +109,7 @@ def pressure(self, temperature, volume, params):
             * np.exp(bi)
             / ai
             * np.power(bi, ci / ai)
-            * (
-                gammaincc(
-                    -ci / ai,
-                    bi * np.exp(ai * lnVrel),
-                )
-                - gammaincc(-ci / ai, bi)
-            )
+            * (generalised_gammainc(-ci / ai, bi * np.exp(ai * lnVrel), bi))
         )
 
     def molar_internal_energy(self, pressure, temperature, volume, params):
@@ -127,18 +132,9 @@ def molar_internal_energy(self, pressure, temperature, volume, params):
         I1 = (
             np.power(bi, 1.0 / ai)
             * Vrel
-            * (
-                gammaincc(
-                    -ci / ai,
-                    bi * np.exp(ai * lnVrel),
-                )
-                - gammaincc(-ci / ai, bi)
-            )
+            * (generalised_gammainc(-ci / ai, bi * np.exp(ai * lnVrel), bi))
         )
-        I2 = gammaincc(
-            (1.0 - ci) / ai,
-            bi * np.exp(ai * lnVrel),
-        ) - gammaincc((1.0 - ci) / ai, bi)
+        I2 = generalised_gammainc((1.0 - ci) / ai, bi * np.exp(ai * lnVrel), bi)
 
         return params["E_0"] + params["P_0"] * (volume - params["V_0"]) + f * (I1 - I2)
 
@@ -221,9 +217,22 @@ def validate_parameters(self, params):
         if params["Kprime_0"] < 0.0 or params["Kprime_0"] > 10.0:
             warnings.warn("Unusual value for Kprime_0", stacklevel=2)
         if params["Kdprime_0"] > 0.0:
-            warnings.warn("Unusual value for Kdprime_0", stacklevel=2)
+            warnings.warn("Kdprime_0 should be negative", stacklevel=2)
+        if (
+            -params["Kdprime_0"] * params["K_0"]
+            < params["Kprime_0"] - params["Kprime_inf"]
+        ):
+            warnings.warn(
+                "Kdprime_0*K_0 is expected to be more "
+                "negative than (Kprime_0 - Kprime_inf)",
+                stacklevel=2,
+            )
         if (
             params["Kprime_inf"] < 5.0 / 3.0
             or params["Kprime_inf"] > params["Kprime_0"]
         ):
-            warnings.warn("Unusual value for Kprime_inf", stacklevel=2)
+            warnings.warn(
+                "Kprime_inf is expected to be greater "
+                "than the Thomas-Fermi limit (5/3)",
+                stacklevel=2,
+            )