Issues with Sampling from Power-Law Potential DFs

[slizewsk]

🐛 Bug: Sampling Fails for PowerSphericalPotential

Attempting to sample from isotropic and constant anisotropy distribution functions (eddingtondf, constantbetadf) based on PowerSphericalPotential fails with a divide by zero error.

This issue appears both in the self-consistent case and when using a separate tracer density. While the potentials do diverge at very small r, the class should be designed to handle this appropriately.

Reproducible Example

For the self-consistent case:

from galpy.potential import PowerSphericalPotential
from galpy.df import eddingtondf

gamma = 0.5
sph_pot = PowerSphericalPotential(amp=1, alpha=gamma + 2)
iso_df = eddingtondf(pot=sph_pot, rmax=1e4)

samps = iso_df.sample(n=3000, return_orbit=True)  # Fails here

Full traceback: https://pastebin.com/5Z32qx0u

Expected Behavior

• The DF should allow bound orbit sampling for both isotropic and anisotropic forms.
• This issue does not occur for isotropicHernquistdf and constantbetaHernquistdf. See https://pastebin.com/XcfnDTtt for a working example using the Hernquist model to show the

System Details

• macOS-15.3.1-arm64-arm-64bit
• Python 3.11.8 | packaged by conda-forge | (main, Feb 16 2024, 20:49:36) [Clang 16.0.6 ]
• numpy 1.26.4
• scipy 1.13.1
• matplotlib 3.8.4
• galpy 1.10.2
• astropy 6.1.0

Additional Context

• I am testing PowerSphericalPotential for a Bayesian model fitting package using CmdStanPy.
• The potential is defined as: $\Phi(r) = -\frac{\Phi_0}{r^\gamma}$
• With a self-consistent density: $\rho(r) = \frac{\gamma (1-\gamma) \Phi_0}{4\pi G} r^{-\gamma-2}$
• In galpy, PowerSphericalPotential(alpha=γ+2) should correspond to this density.
• Sampling also fails for the KeplerPotential with the same error.
• For the NFWPotential, sampling works for the isotropic DF but fails for the constant-anisotropy DF. However, this produces a different error.
• I noticed the eddingtondf ._Emin returns -inf.. perhaps related?

Let me know if you'd like more details or a different test case.

[github-actions]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. If the issue has been resolved since the last activity, please close the issue. Thank you for your contributions.