Add Gaussian, TriaxialGaussian, and AxisymmetricGaussian potentials

[prashjet]

Adds three new potential classes:

• GaussianPotential: spherical Gaussian density profile, closed-form via erf.
• TriaxialGaussianPotential: Chandrasekhar ellipsoidal quadrature, mirroring the NFW/TriaxialNFW structure.
• AxisymmetricGaussianPotential: the q1=1 special case of TriaxialGaussianPotential as its own class (single shape parameter q2), for API clarity in the common oblate/prolate case.

All three are cross-checked against galpy in tests; AxisymmetricGaussianPotential is additionally cross-checked against TriaxialGaussianPotential(q1=1). Profiling showed no measurable speedup for the axisymmetric class over the triaxial one with q1 fixed — it earns its keep on API ergonomics alone.

Note: this overlaps with #756, which adds a spherical GaussianPotential. This PR additionally adds triaxial and axisymmetric variants under a gaussian/ package rather than a single module. Happy to rebase on #756 or coordinate however maintainers prefer.

[Copilot]

Pull request overview

Adds a new Gaussian-density potential family to galax.potential, including spherical, triaxial, and axisymmetric variants, and introduces unit tests that cross-check against galpy for validation.

Changes:

• Added GaussianPotential with closed-form potential/mass/density helpers.
• Added TriaxialGaussianPotential (ellipsoidal quadrature) and AxisymmetricGaussianPotential (q1=1 specialization).
• Added galpy cross-check helpers and new unit tests for all three potentials.

Reviewed changes

Copilot reviewed 10 out of 10 changed files in this pull request and generated 3 comments.

Show a summary per file

File	Description
tests/unit/potential/builtin/test_triaxialgaussian.py	Adds unit tests for TriaxialGaussianPotential, including galpy equivalence checks.
tests/unit/potential/builtin/test_gaussian.py	Adds unit tests for GaussianPotential, including galpy equivalence checks.
tests/unit/potential/builtin/test_axisymmetricgaussian.py	Adds unit tests for AxisymmetricGaussianPotential, including galpy equivalence checks.
tests/unit/potential/builtin/test_common.py	Adds a shared helper to compare Gaussian(-family) potentials against galpy.
src/galax/potential/_src/builtin/gaussian/base.py	Implements spherical GaussianPotential and closed-form helper functions.
src/galax/potential/_src/builtin/gaussian/triaxial.py	Implements TriaxialGaussianPotential via Gauss–Legendre quadrature.
src/galax/potential/_src/builtin/gaussian/axisymmetric.py	Implements AxisymmetricGaussianPotential as the q1=1 quadrature specialization.
src/galax/potential/_src/builtin/gaussian/init.py	Exposes the Gaussian potential classes from the gaussian package.
src/galax/potential/_src/builtin/init.py	Re-exports the new Gaussian potential classes from the builtin collection.
src/galax/potential/init.py	Re-exports the new Gaussian potential classes at the top-level galax.potential API.

Suppressed comments (2)

src/galax/potential/_src/builtin/gaussian/triaxial.py:66

• m is documented as the “Total mass”, but the implementation (via rho0 = m / ((2π)^{3/2} r_s^3) and the ellipsoidal coordinates) implies the true total mass scales with shape: M_tot = m * q1 * q2 (also reflected in the galpy cross-check helper). This should be clarified to avoid users passing a total mass and getting a different normalization.

        dimensions="mass", doc="Total mass of the potential."

src/galax/potential/_src/builtin/gaussian/axisymmetric.py:83

• m is documented as the “Total mass”, but (as in the triaxial class) the normalization implied by rho0 = m / ((2π)^{3/2} r_s^3) means the true total mass depends on flattening: M_tot = m * q2 for the axisymmetric case. Clarifying this will prevent accidental mis-normalization by callers.

        dimensions="mass", doc="Total mass of the potential."

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[prashjet]

Thanks for the thorough review! Pushed two follow-up commits addressing the feedback:

Mass parameter (@nstarman, Copilot):

• Renamed m → m_tot on all three classes.
• For TriaxialGaussianPotential/AxisymmetricGaussianPotential, fixed the density normalization to divide by q1*q2 (resp. q2), so m_tot is the true total mass for any flattening — not just in the spherical case. Previously the parameter silently meant "mass scale," matching the (infinite-mass) NFW convention, which is misleading for a profile with finite mass.

r=0 NaN (Copilot):

• Fixed the 0/0 in GaussianPotential.potential() by substituting the analytic limit Φ(0) = -G m_tot √(2/π) / r_s explicitly. Note: gradient/hessian at exactly r=0 are still NaN — that's a separate, pre-existing limitation shared by every radial potential in the codebase (e.g. HernquistPotential has the same issue), coming from vector_norm not being differentiable at the origin. Out of scope here.

Stray G in density docstrings (Copilot):

• Fixed in TriaxialGaussianPotential/AxisymmetricGaussianPotential.

Duplicated GaussLegendreIntegrator (@nstarman):

• Consolidated into potential/_src/utils.py alongside the other cross-potential helpers, with a GaussLegendreIntegrator.for_order() factory so each class's __post_init__/integration_order boilerplate collapses to one line. Used by TriaxialNFWPotential too now.

AbstractGaussianPotential (@nstarman): replied inline — holding off since there's no AbstractNFWPotential either despite that family having more variants; happy to do a follow-up covering both families if that's the preferred direction.

All existing + updated tests pass (1026 passed locally).

[nstarman]

@prashjet if you rebase on main it should fix the Format CI job.

[Copilot]

Pull request overview

Copilot reviewed 12 out of 12 changed files in this pull request and generated 1 comment.

Suppressed comments (1)

src/galax/potential/_src/builtin/gaussian/base.py:104

• The density expression uses an undefined symbol u; every other Gaussian density formula in this module (e.g. the standalone density function docstring and the class docstring) expresses the exponent as r^2 / (2 r_s^2). Using u here is inconsistent and, since u is also the imported unxt alias in this file, potentially confusing.

            \rho(r) = \rho_0 \exp\left(-\frac{u^2}{2}\right)

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.14%. Comparing base (c8eeeff) to head (c25e5af).

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #796      +/-   ##
==========================================
+ Coverage   96.00%   96.14%   +0.13%     
==========================================
  Files         153      157       +4     
  Lines        6079     6298     +219     
==========================================
+ Hits         5836     6055     +219     
  Misses        243      243              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Copilot]

Pull request overview

Copilot reviewed 13 out of 13 changed files in this pull request and generated no new comments.

Suppressed comments (2)

src/galax/potential/_src/builtin/gaussian/base.py:192

• This doctest is missing the blank line between the "A quick sanity check:" prose and the >>> prompt, unlike rho0_of_m just above (lines 169-171). In reStructuredText a doctest block must be preceded by a blank line, so as written the >>> will render as part of the preceding paragraph rather than as an interactive example. Add the blank line for consistency and correct rendering.

    A quick sanity check:
    >>> m_of_rho0({"rho0": 1.0, "r_s": 1.0}) - (2 * jnp.pi) ** 1.5

tests/unit/potential/builtin/test_axisymmetricgaussian.py:106

• The PR description states that AxisymmetricGaussianPotential is "additionally cross-checked against TriaxialGaussianPotential(q1=1)", but no such equivalence test is present in this file (or anywhere in the test suite). The committed tests only compare against hard-coded expected values and against galpy; and TestTriaxialGaussianPotential uses q1=1.1, so the two are never exercised at the equal configuration. Since this class is justified purely as the q1=1 special case of TriaxialGaussianPotential, an explicit equivalence test would be the most valuable guard against regressions. Either add that cross-check test or update the PR description to match what is committed.

    @pytest.mark.skipif(not OptDeps.GALPY.installed, reason="requires galpy")
    def test_method_galpy(
        self, pot: gp.AxisymmetricGaussianPotential, x: gt.QuSz3
    ) -> None:
        """Test the equivalence of potential/density between galpy and galax."""
        assert_gaussian_matches_galpy(pot, x)

[nstarman]

Thanks! @prashjet. LMK what other features would be useful for your project.