From db81da095e58ed32a5a299d3446dab3c31720f8c Mon Sep 17 00:00:00 2001
From: Bob Myhill <myhill.bob@gmail.com>
Date: Tue, 26 Nov 2024 10:51:17 +0000
Subject: [PATCH] tweak spock eos

---
 burnman/eos/spock.py | 42 +++++++++++++++++++-----------------------
 1 file changed, 19 insertions(+), 23 deletions(-)

diff --git a/burnman/eos/spock.py b/burnman/eos/spock.py
index 98c6b580..9492bbdb 100644
--- a/burnman/eos/spock.py
+++ b/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)