Address MultipoleExpansionPotential review comments

- Remove k parameter from public API (BPoly sets degree, old docstring misleading)
- Skip sin integral computation for m=0 (never used, zero bytes reach C)
- Enhance SCF test to verify only l=0,1 coefficients are non-zero
- Add C coverage tests for below-grid log branch (line 255), d²R (lines 286-288), and non-axi potential evaluation

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: jobovy

galpy/potential/MultipoleExpansionPotential.py

@@ -45,16 +45,17 @@ class MultipoleExpansionPotential(Potential, SphericalHarmonicPotentialMixin):
     def __init__(
         self,
         amp=1.0,
-        dens=lambda R, z: 1.0
-        / (2.0 * numpy.pi)
-        / numpy.sqrt(R**2 + z**2)
-        / (1 + numpy.sqrt(R**2 + z**2)) ** 3,
+        dens=lambda R, z: (
+            1.0
+            / (2.0 * numpy.pi)
+            / numpy.sqrt(R**2 + z**2)
+            / (1 + numpy.sqrt(R**2 + z**2)) ** 3
+        ),
         L=6,
         rgrid=numpy.geomspace(1e-3, 30, 1_001),
         symmetry=None,
         costheta_order=None,
         phi_order=None,
-        k=3,
         normalize=False,
         ro=None,
         vo=None,
@@ -78,8 +79,6 @@ def __init__(
             Gauss-Legendre quadrature order for theta. Default: ``max(20, L+1)``.
         phi_order : int, optional
             Number of uniform phi points for trapezoidal rule. Default: ``max(20, 2*L+1)``.
-        k : int, optional
-            Spline interpolation degree for radial functions (default: 3). Use k=5 for smoother second derivatives.
         normalize : bool or float, optional
             If True, normalize such that vc(1.,0.)=1., or, if given as a number, such that the force is this fraction of the force necessary to make vc(1.,0.)=1.
         ro : float or Quantity, optional
@@ -95,7 +94,7 @@ def __init__(
         # Parse density function
         dens_func = self._parse_density(dens)
         self._rgrid = rgrid
-        self._k = k
+        self._k = 3
         self._L = L
         # Set M based on symmetry
         if symmetry is not None and symmetry.startswith("spher"):
@@ -136,7 +135,7 @@ def __init__(
         self._rho_cos_splines = [
             [
                 InterpolatedUnivariateSpline(
-                    rgrid, beta_lm[l, m] * self._rho_cos[:, l, m], k=k
+                    rgrid, beta_lm[l, m] * self._rho_cos[:, l, m], k=self._k
                 )
                 for m in range(M)
             ]
@@ -145,7 +144,7 @@ def __init__(
         self._rho_sin_splines = [
             [
                 InterpolatedUnivariateSpline(
-                    rgrid, beta_lm[l, m] * self._rho_sin[:, l, m], k=k
+                    rgrid, beta_lm[l, m] * self._rho_sin[:, l, m], k=self._k
                 )
                 for m in range(M)
             ]
@@ -374,18 +373,14 @@ def _precompute_radial_integrals(self, beta_lm):
             pref = -4.0 * numpy.pi / (2 * l + 1)
             for m in range(min(l + 1, M)):
                 pref_blm = pref * beta_lm[l, m]
-                for rho_arr, I_inner_store, I_outer_store in [
-                    (
-                        self._rho_cos[:, l, m],
-                        self._I_inner_cos,
-                        self._I_outer_cos,
-                    ),
-                    (
-                        self._rho_sin[:, l, m],
-                        self._I_inner_sin,
-                        self._I_outer_sin,
-                    ),
-                ]:
+                # sin(m*phi) = 0 for m=0, so sin integrals are identically zero
+                # and never accessed; skip their computation.
+                pairs = [(self._rho_cos[:, l, m], self._I_inner_cos, self._I_outer_cos)]
+                if m > 0:
+                    pairs.append(
+                        (self._rho_sin[:, l, m], self._I_inner_sin, self._I_outer_sin)
+                    )
+                for rho_arr, I_inner_store, I_outer_store in pairs:
                     # Inner integral: integrand = r^{l+2} * rho_lm(r)
                     f_inner = rgrid ** (l + 2) * rho_arr
                     # Use spline integration for higher accuracy than trapezoid

tests/test_MultipoleExpansionPotential.py

@@ -120,6 +120,48 @@ def test_scf_potential_cross_validation():
         )
 
 
+def test_scf_only_l0_l1_nonzero():
+    """For a density with only l=0 and l=1 harmonics, verify the multipole expansion
+    has negligible l>=2 coefficients and correct l=0, l=1 radial functions."""
+    Acos = numpy.zeros((3, 3, 1))
+    Acos[0, 0, 0] = 1.0
+    Acos[1, 0, 0] = 0.1
+    Acos[0, 1, 0] = 0.05
+    scf = SCFPotential(Acos=Acos, a=1.0)
+    mp = MultipoleExpansionPotential(
+        dens=scf, L=6, symmetry="axisymmetric", rgrid=_FINE_RGRID
+    )
+    # Only l=0 and l=1 should have non-negligible raw coefficients
+    max_l0 = numpy.max(numpy.abs(mp._rho_cos[:, 0, 0]))
+    max_l1 = numpy.max(numpy.abs(mp._rho_cos[:, 1, 0]))
+    assert max_l0 > 0, "l=0 coefficient must be non-zero"
+    assert max_l1 > 0, "l=1 coefficient must be non-zero"
+    for l in range(2, mp._L):
+        max_lx = numpy.max(numpy.abs(mp._rho_cos[:, l, 0]))
+        assert max_lx < 1e-6 * max_l0, (
+            f"l={l} coefficient should be negligible: {max_lx:.2e} vs l=0 max {max_l0:.2e}"
+        )
+    # Verify l=0 radial function: at z=0, P_1^0(0) = 0 so only l=0 contributes to density.
+    # Therefore mp.dens(R, 0) should equal scf.dens(R, 0) to high precision.
+    for R in [0.5, 1.0, 2.0]:
+        d_mp = mp.dens(R, 0.0, use_physical=False)
+        d_scf = scf.dens(R, 0.0, use_physical=False)
+        assert abs(d_mp - d_scf) / abs(d_scf) < 1e-5, (
+            f"l=0 radial function mismatch at midplane R={R}: mp={d_mp}, scf={d_scf}"
+        )
+    # Verify l=1 radial function: the l=1 term breaks z-symmetry (cos θ changes sign).
+    # The difference dens(R,z) - dens(R,-z) isolates the odd-l (here l=1) contribution.
+    for R, z in [(1.0, 1.0), (0.5, 1.0), (2.0, 0.5)]:
+        diff_mp = mp.dens(R, z, use_physical=False) - mp.dens(R, -z, use_physical=False)
+        diff_scf = scf.dens(R, z, use_physical=False) - scf.dens(
+            R, -z, use_physical=False
+        )
+        assert abs(diff_scf) > 0, "SCF l=1 term should give z-asymmetric density"
+        assert abs(diff_mp - diff_scf) / abs(diff_scf) < 1e-4, (
+            f"l=1 radial function mismatch at R={R}, z={z}: diff_mp={diff_mp}, diff_scf={diff_scf}"
+        )
+
+
 def test_scf_density_cross_validation():
     Acos = numpy.zeros((3, 3, 1))
     Acos[0, 0, 0] = 1.0
@@ -227,9 +269,9 @@ def test_2arg_lambda_input():
     # rho = amp/(4*pi) * a / (r * (r+a)^3) for HernquistPotential(amp=2, a=1)
     coeff = 1.0 / (2.0 * numpy.pi)
     mp = MultipoleExpansionPotential(
-        dens=lambda R, z: coeff
-        / numpy.sqrt(R**2 + z**2)
-        / (1 + numpy.sqrt(R**2 + z**2)) ** 3,
+        dens=lambda R, z: (
+            coeff / numpy.sqrt(R**2 + z**2) / (1 + numpy.sqrt(R**2 + z**2)) ** 3
+        ),
         L=2,
         symmetry="spherical",
         rgrid=_FINE_RGRID,
@@ -358,9 +400,9 @@ def test_3arg_callable_density_input():
     """Test that a 3-argument callable density (R, z, phi) without units works."""
     coeff = 1.0 / (2.0 * numpy.pi)
     mp = MultipoleExpansionPotential(
-        dens=lambda R, z, phi: coeff
-        / numpy.sqrt(R**2 + z**2)
-        / (1 + numpy.sqrt(R**2 + z**2)) ** 3,
+        dens=lambda R, z, phi: (
+            coeff / numpy.sqrt(R**2 + z**2) / (1 + numpy.sqrt(R**2 + z**2)) ** 3
+        ),
         L=4,
         symmetry=None,
         rgrid=_DEFAULT_RGRID,
@@ -422,10 +464,12 @@ def test_2nd_derivs_on_z_axis():
     where dP/d(costheta) diverges for m>0, triggering the pole clamping."""
     coeff = 1.0 / (2.0 * numpy.pi)
     mp = MultipoleExpansionPotential(
-        dens=lambda R, z, phi: coeff
-        / numpy.sqrt(R**2 + z**2)
-        / (1 + numpy.sqrt(R**2 + z**2)) ** 3
-        * (1.0 + 0.1 * numpy.cos(2 * phi)),
+        dens=lambda R, z, phi: (
+            coeff
+            / numpy.sqrt(R**2 + z**2)
+            / (1 + numpy.sqrt(R**2 + z**2)) ** 3
+            * (1.0 + 0.1 * numpy.cos(2 * phi))
+        ),
         L=6,
         symmetry=None,
         rgrid=_FINE_RGRID,
@@ -478,21 +522,14 @@ def test_spherical_2nd_derivs_match_hernquist():
             )
 
 
-def test_spline_degree_k_parameter():
-    """Test that the k parameter is passed through to splines."""
+def test_internal_spline_degree():
+    """Test that the internal spline degree is set to 3."""
     hp = HernquistPotential(amp=2.0, a=1.0)
-    mp3 = MultipoleExpansionPotential(
-        dens=hp, L=2, symmetry="spherical", rgrid=_FINE_RGRID, k=3
-    )
-    mp5 = MultipoleExpansionPotential(
-        dens=hp, L=2, symmetry="spherical", rgrid=_FINE_RGRID, k=5
+    mp = MultipoleExpansionPotential(
+        dens=hp, L=2, symmetry="spherical", rgrid=_FINE_RGRID
     )
-    assert mp3._k == 3
-    assert mp5._k == 5
-    # Both should give reasonable results
-    for mp in [mp3, mp5]:
-        val = mp.R2deriv(1.0, 0.5, use_physical=False)
-        assert numpy.isfinite(val)
+    assert mp._k == 3
+    assert numpy.isfinite(mp.R2deriv(1.0, 0.5, use_physical=False))
 
 
 # --- Below/above grid extrapolation tests ---
@@ -581,3 +618,50 @@ def test_below_grid_density_clamped():
     assert d_below == d_at_rmin, (
         f"Density below grid should be clamped to rmin value: {d_below} != {d_at_rmin}"
     )
+
+
+# --- C code coverage ---
+
+
+def test_c_orbit_below_grid_l2():
+    """Cover C code line 255 (l=2 below-grid log formula) via orbit integration in C.
+    Orbit integration uses Rforce/zforce in C, which calls below_grid_integrals
+    for r < rmin.  With L=4, l=2 is included and triggers the log-branch."""
+    from galpy.orbit import Orbit
+
+    hp = HernquistPotential(amp=2.0, a=1.0)
+    mp = MultipoleExpansionPotential(
+        dens=hp,
+        L=4,  # includes l=0,1,2,3; l=2 hits the log branch in below_grid_integrals
+        symmetry="axisymmetric",
+        rgrid=numpy.geomspace(2.0, 20.0, 201),  # rmin=2 > orbit's minimum r
+    )
+    # Orbit starting at R=1.0 < rmin=2.0 so forces are evaluated below the grid
+    o = Orbit([1.0, 0.1, 1.0, 0.0, 0.1, 0.0])
+    ts = numpy.linspace(0, 10, 101)
+    o.integrate(ts, mp, method="leapfrog_c")
+    assert numpy.all(numpy.isfinite(o.R(ts))), "Orbit R should be finite"
+    assert numpy.all(numpy.isfinite(o.z(ts))), "Orbit z should be finite"
+
+
+def test_c_planar_liouville_below_grid_d2R():
+    """Cover C code lines 286-288 (d²R below-grid via EVAL_DERIV2 mode).
+    A planar orbit.integrate_dxdv with a C integrator calls
+    integratePlanarOrbit_dxdv in C, which calls MultipoleExpansionPotential-
+    PlanarR2deriv → compute_multipole_spher_2nd_derivs → eval_radial_lm with
+    EVAL_DERIV2. With rmin=2 > orbit radius the below-grid branch is hit."""
+    from galpy.orbit import Orbit
+
+    hp = HernquistPotential(amp=2.0, a=1.0)
+    mp = MultipoleExpansionPotential(
+        dens=hp,
+        L=4,
+        symmetry="axisymmetric",
+        rgrid=numpy.geomspace(2.0, 20.0, 201),  # rmin=2 > orbit's R=1.0
+    )
+    # Planar orbit [R, vR, vT, phi] at R=1.0 < rmin=2.0
+    o = Orbit([1.0, 0.1, 1.0, 0.5])
+    ts = numpy.linspace(0, 5, 51)
+    o.integrate_dxdv([1.0, 0.0, 0.0, 0.0], ts, mp, method="dopr54_c")
+    result = o.getOrbit_dxdv()
+    assert numpy.all(numpy.isfinite(result)), "integrate_dxdv result should be finite"

tests/test_actionAngle.py

@@ -3552,6 +3552,15 @@ def test_actionAngleStaeckel_wSpherical_conserved_actions_c():
         symmetry="spherical",
         normalize=1.0,
     )
+    mep_nonaxi = potential.MultipoleExpansionPotential(
+        dens=lambda R, z, phi: potential.HernquistPotential(normalize=1.0).dens(
+            R, z, phi
+        )
+        * (1.0 + 1e-9 * numpy.cos(phi)),
+        L=2,
+        symmetry=None,
+        normalize=1.0,
+    )
     pots = [
         lp,
         lpb,
@@ -3579,6 +3588,7 @@ def test_actionAngleStaeckel_wSpherical_conserved_actions_c():
         tpsp,
         tpsp_beta3,
         mep,
+        mep_nonaxi,
     ]
     for pot in pots:
         aAS = actionAngleStaeckel(pot=pot, c=True, delta=0.01)

tests/test_dynamfric.py

@@ -290,10 +290,10 @@ def test_dynamfric_c():
             normalize=1.0,
         ),
         potential.MultipoleExpansionPotential(
-            dens=lambda R, z, phi: potential.HernquistPotential(
-                normalize=1.0, a=3.5
-            ).dens(R, z)
-            * (1 + 0.01 * (numpy.cos(phi) + numpy.sin(phi))),
+            dens=lambda R, z, phi: (
+                potential.HernquistPotential(normalize=1.0, a=3.5).dens(R, z)
+                * (1 + 0.01 * (numpy.cos(phi) + numpy.sin(phi)))
+            ),
             L=2,
             normalize=1.0,
         ),

tests/test_quantity.py

@@ -10324,12 +10324,14 @@ def test_potential_paramunits():
     )
     # MultipoleExpansionPotential with density in units (2-arg)
     pot = potential.MultipoleExpansionPotential(
-        dens=lambda R, z: coeff
-        / numpy.sqrt(R**2 + z**2)
-        / (1 + numpy.sqrt(R**2 + z**2)) ** 3
-        * dens_scale
-        * units.Msun
-        / units.pc**3,
+        dens=lambda R, z: (
+            coeff
+            / numpy.sqrt(R**2 + z**2)
+            / (1 + numpy.sqrt(R**2 + z**2)) ** 3
+            * dens_scale
+            * units.Msun
+            / units.pc**3
+        ),
         symmetry="spherical",
         rgrid=numpy.geomspace(1e-2, 20, 201),
         ro=ro,
@@ -10346,12 +10348,14 @@ def test_potential_paramunits():
     )
     # MultipoleExpansionPotential with density in units (3-arg)
     pot = potential.MultipoleExpansionPotential(
-        dens=lambda R, z, phi: coeff
-        / numpy.sqrt(R**2 + z**2)
-        / (1 + numpy.sqrt(R**2 + z**2)) ** 3
-        * dens_scale
-        * units.Msun
-        / units.pc**3,
+        dens=lambda R, z, phi: (
+            coeff
+            / numpy.sqrt(R**2 + z**2)
+            / (1 + numpy.sqrt(R**2 + z**2)) ** 3
+            * dens_scale
+            * units.Msun
+            / units.pc**3
+        ),
         L=4,
         symmetry=None,
         rgrid=numpy.geomspace(1e-2, 20, 201),