Repository: jobovy/galpy

Issue:
TwoPowerSphericalPotential returns NaN at R=0 for non-special cases
## 🐛 Bug

The TwoPowerSphericalPotential returns NaN outside of the built-in combinations of alpha and beta (e.g., alpha=1 and beta=3). This causes issues when initializing distribution functions like the Osipkov-Merritt DF. 


### Reproducible example

``` python
import numpy as np
from galpy.potential import TwoPowerSphericalPotential, PlummerPotential
from galpy.df import osipkovmerrittdf

pot = TwoPowerSphericalPotential(
    amp=1,
    a=1,
    alpha=0.0,
    beta=3,
)
dens_pot = PlummerPotential(amp=1, b=1)
print(pot(R=0, z=0))
DF = osipkovmerrittdf(pot=pot, denspot=dens_pot)  
# Fails here with ValueError: x must be increasing if s > 0
```

### Expected behavior

For alpha < 2, I believe the potential should be finite at R=0.


### System Details

macOS-12.4-x86_64-i386-64bit-Mach-O
Python 3.13.3 | packaged by conda-forge | (main, Apr 14 2025, 20:46:04) [Clang 18.1.8 ]
numpy 2.2.6
scipy 1.15.2
matplotlib 3.10.3
galpy 1.10.2
astropy 7.0.2


### Additional context

I'm primarily concerned with cases where 0 <= alpha < 2 and beta = 3 if addressing that regime is simpler than solving the issue for arbitrary alpha and beta.

Discussion:

Solution patch:
diff --git a/galpy/potential/TwoPowerSphericalPotential.py b/galpy/potential/TwoPowerSphericalPotential.py
index 6547ac29..af504dce 100644
--- a/galpy/potential/TwoPowerSphericalPotential.py
+++ b/galpy/potential/TwoPowerSphericalPotential.py
@@ -95,9 +95,8 @@ class TwoPowerSphericalPotential(Potential):
             return (
                 (1.0 / self.a)
                 * (
-                    r
-                    - self.a
-                    * (r / self.a) ** (3.0 - self.alpha)
+                    1
+                    - (r / self.a) ** (2.0 - self.alpha)
                     / (3.0 - self.alpha)
                     * special.hyp2f1(
                         3.0 - self.alpha,
@@ -107,10 +106,11 @@ class TwoPowerSphericalPotential(Potential):
                     )
                 )
                 / (self.alpha - 2.0)
-                / r
             )
         else:
-            r = numpy.sqrt(R**2.0 + z**2.0)
+            r = (
+                numpy.sqrt(R**2.0 + z**2.0) + 1e-11
+            )  # avoid division by zero and numerical instability of the hyp2f1 function
             return (
                 special.gamma(self.beta - 3.0)
                 * (


Test patch:
diff --git a/tests/test_potential.py b/tests/test_potential.py
index dbd997b3..1fcf4dd4 100644
--- a/tests/test_potential.py
+++ b/tests/test_potential.py
@@ -40,6 +40,7 @@ def test_normalize_potential():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialPowerSphericalPotential")
     pots.append("specialFlattenedPowerPotential")
@@ -139,6 +140,7 @@ def test_forceAsDeriv_potential():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -392,6 +394,7 @@ def test_2ndDeriv_potential():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -644,6 +647,8 @@ def test_2ndDeriv_potential():
                         continue  # Not implemented, or badly defined
                     if p == "NFWTwoPowerSphericalPotential":
                         continue  # Not implemented, or badly defined
+                    if p == "CoreNFWTwoPowerSphericalPotential":
+                        continue  # Not implemented, or badly defined
                     # Excluding KuzminDiskPotential at z = 0
                     if p == "KuzminDiskPotential" and Zs[jj] == 0:
                         continue
@@ -783,6 +788,7 @@ def test_poisson_potential():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -932,6 +938,7 @@ def test_poisson_surfdens_potential():
     pots.append('HernquistTwoPowerSphericalPotential')
     pots.append('JaffeTwoPowerSphericalPotential')
     pots.append('NFWTwoPowerSphericalPotential')
+    pots.append('CoreNFWTwoPowerSphericalPotential')
     pots.append('specialMiyamotoNagaiPotential')
     pots.append('specialMN3ExponentialDiskPotentialPD')
     pots.append('specialMN3ExponentialDiskPotentialSECH')
@@ -1090,6 +1097,7 @@ def test_evaluateAndDerivs_potential():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -1348,6 +1356,7 @@ def test_amp_mult_divide():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -1803,12 +1812,13 @@ def test_potential_at_zero():
             and not "toVertical" in p
         )
     ]
-    # pots.append('specialTwoPowerSphericalPotential')
-    # pots.append('DehnenTwoPowerSphericalPotential')
-    # pots.append('DehnenCoreTwoPowerSphericalPotential')
-    # pots.append('HernquistTwoPowerSphericalPotential')
-    # pots.append('JaffeTwoPowerSphericalPotential')
-    # pots.append('NFWTwoPowerSphericalPotential') # Difficult, and who cares?
+    pots.append("specialTwoPowerSphericalPotential")
+    pots.append("DehnenTwoPowerSphericalPotential")
+    pots.append("DehnenCoreTwoPowerSphericalPotential")
+    pots.append("HernquistTwoPowerSphericalPotential")
+    pots.append("JaffeTwoPowerSphericalPotential")
+    pots.append("NFWTwoPowerSphericalPotential")
+    pots.append("CoreNFWTwoPowerSphericalPotential")
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -1961,6 +1971,7 @@ def test_potential_at_infinity():
     pots.append("HernquistTwoPowerSphericalPotential")
     pots.append("JaffeTwoPowerSphericalPotential")
     # pots.append('NFWTwoPowerSphericalPotential') # Difficult, and who cares?
+    # pots.append('CoreNFWTwoPowerSphericalPotential')
     pots.append("specialMiyamotoNagaiPotential")
     pots.append("specialMN3ExponentialDiskPotentialPD")
     pots.append("specialMN3ExponentialDiskPotentialSECH")
@@ -8306,6 +8317,12 @@ class NFWTwoPowerSphericalPotential(TwoPowerSphericalPotential):
         return None
 
 
+class CoreNFWTwoPowerSphericalPotential(TwoPowerSphericalPotential):
+    def __init__(self):
+        TwoPowerSphericalPotential.__init__(self, amp=1.0, a=5.0, alpha=0.0, beta=3.0)
+        return None
+
+
 class specialPowerSphericalPotential(PowerSphericalPotential):
     def __init__(self):
         PowerSphericalPotential.__init__(self, amp=1.0, alpha=2.0)