Add differential energy distribution for spherical DFs

HISTORY.txt

@@ -5,6 +5,9 @@ v1.9.0 (XXXX-XX-XX)
   add Orbit.SOS and Orbit.plotSOS methods to compute and plot the surface of section
   of an orbit.
 
+- Added a method, dMdE, to calculate the differential energy distribution of
+  spherical distribution functions.
+
 - Started using black for code formatting.
 
 - Allow potentials' density and surface density to be plotted on physical axes.

doc/source/reference/df.rst

@@ -78,6 +78,7 @@ General instance routines
    :maxdepth: 1
 
    __call__ <sphericaldfcall.rst>
+   dMdE <sphericaldfdmde.rst>
    beta <sphericaldfbeta.rst>
    sigmar <sphericaldfsigmar.rst>
    sigmar <sphericaldfsigmat.rst>

doc/source/reference/sphericaldfdmde.rst

@@ -0,0 +1,4 @@
+galpy.df.sphericaldf.dMdE
+=========================
+
+.. automethod:: galpy.df.sphericaldf.dMdE

galpy/df/constantbetadf.py

@@ -76,6 +76,33 @@ def _call_internal(self, *args):
         E, L, _ = args
         return L ** (-2 * self._beta) * self.fE(E)
 
+    def _dMdE(self, E):
+        if not hasattr(self, "_rphi"):
+            self._rphi = self._setup_rphi_interpolator()
+        fE = self.fE(E)
+        out = numpy.zeros_like(E)
+        out[fE > 0.0] = (
+            (2.0 * numpy.pi) ** 2.5
+            * special.gamma(1.0 - self._beta)
+            / 2.0 ** (self._beta - 1.0)
+            / special.gamma(1.5 - self._beta)
+            * fE[fE > 0.0]
+            * numpy.array(
+                [
+                    integrate.quad(
+                        lambda r: r ** (2.0 - 2.0 * self._beta)
+                        * (tE - _evaluatePotentials(self._pot, r, 0.0))
+                        ** (0.5 - self._beta),
+                        0.0,
+                        self._rphi(tE),
+                    )[0]
+                    for ii, tE in enumerate(E)
+                    if fE[ii] > 0.0
+                ]
+            )
+        )
+        return out
+
     def _sample_eta(self, r, n=1):
         """Sample the angle eta which defines radial vs tangential velocities"""
         if not hasattr(self, "_coseta_icmf_interp"):

galpy/df/isotropicHernquistdf.py

@@ -97,6 +97,23 @@ def fE(self, E):
             fE[Etilde_out] = 0
         return fE
 
+    def _dMdE(self, E):
+        # E already in internal units here
+        fE = self.fE(E)
+        A = -self._psi0 / E[fE > 0.0]
+        out = numpy.zeros_like(E)
+        out[fE > 0.0] = (
+            (4.0 * numpy.pi) ** 2.0
+            * fE[fE > 0.0]
+            * self._pot.a**3.0
+            * numpy.sqrt(-2.0 * E[fE > 0.0])
+            * (
+                numpy.sqrt(A - 1.0) * (0.125 * A**2.0 - 5.0 / 12.0 * A - 1.0 / 3.0)
+                + 0.125 * A * (A**2.0 - 4.0 * A + 8.0) * numpy.arccos(A**-0.5)
+            )
+        )
+        return out
+
     def _icmf(self, ms):
         """Analytic expression for the normalized inverse cumulative mass
         function. The argument ms is normalized mass fraction [0,1]"""

galpy/df/osipkovmerrittHernquistdf.py

@@ -85,15 +85,15 @@ def fQ(self, Q):
         """
         Qtilde = numpy.atleast_1d(conversion.parse_energy(Q, vo=self._vo) / self._psi0)
         # Handle potential Q outside of bounds
-        Qtilde_out = numpy.where(numpy.logical_or(Qtilde < 0, Qtilde > 1))[0]
+        Qtilde_out = numpy.where(numpy.logical_or(Qtilde < 0.0, Qtilde > 1.0))[0]
         if len(Qtilde_out) > 0:
             # Dummy variable now and 0 later, prevents numerical issues
             Qtilde[Qtilde_out] = 0.5
         sqrtQtilde = numpy.sqrt(Qtilde)
         # The 'ergodic' part
         fQ = (
             sqrtQtilde
-            / (1 - Qtilde) ** 2.0
+            / (1.0 - Qtilde) ** 2.0
             * (
                 (1.0 - 2.0 * Qtilde) * (8.0 * Qtilde**2.0 - 8.0 * Qtilde - 3.0)
                 + (

galpy/df/osipkovmerrittdf.py

@@ -85,6 +85,64 @@ def _call_internal(self, *args):
         E, L, _ = args
         return self.fQ(-E - 0.5 * L**2.0 / self._ra2)
 
+    def _dMdE(self, E):
+        if not hasattr(self, "_rphi"):
+            self._rphi = self._setup_rphi_interpolator()
+
+        def Lintegrand(t, L2lim, E):
+            return self((E, numpy.sqrt(L2lim - t**2.0)), use_physical=False)
+
+        # Integrate where Q > 0
+
+        out = (
+            16.0
+            * numpy.pi**2.0
+            * numpy.array(
+                [
+                    integrate.quad(
+                        lambda r: r
+                        * integrate.quad(
+                            Lintegrand,
+                            numpy.sqrt(
+                                numpy.amax(
+                                    [
+                                        (0.0),
+                                        (
+                                            2.0
+                                            * r**2.0
+                                            * (
+                                                tE
+                                                - _evaluatePotentials(self._pot, r, 0.0)
+                                            )
+                                            + 2.0 * tE * self._ra2
+                                        ),
+                                    ]
+                                )
+                            ),
+                            numpy.sqrt(
+                                2.0
+                                * r**2.0
+                                * (tE - _evaluatePotentials(self._pot, r, 0.0))
+                            ),
+                            args=(
+                                2.0
+                                * r**2.0
+                                * (tE - _evaluatePotentials(self._pot, r, 0.0)),
+                                tE,
+                            ),
+                        )[0],
+                        0.0,
+                        self._rphi(tE),
+                    )[0]
+                    for ii, tE in enumerate(E)
+                ]
+            )
+        )
+        # Numerical issues can make the integrand's sqrt argument negative, only
+        # happens at dMdE ~ 0, so just set to zero
+        out[numpy.isnan(out)] = 0.0
+        return out
+
     def _sample_eta(self, r, n=1):
         """Sample the angle eta which defines radial vs tangential velocities"""
         # cumulative distribution of x = cos eta satisfies

galpy/df/sphericaldf.py

@@ -168,6 +168,32 @@ def __call__(self, *args, **kwargs):
             L = numpy.atleast_1d(numpy.sqrt(vtotSq - vrad**2.0) * r)
         return self._call_internal(E, L, Lz)  # Some function for each sub-class
 
+    @physical_conversion("energydensity", pop=True)
+    def dMdE(self, E):
+        """
+        NAME:
+
+            dMdE
+
+        PURPOSE:
+
+            Compute the different energy distribution dM/dE: the amount of mass per unit energy
+
+        INPUT:
+
+            E - energy; can be a Quantity
+
+        OUTPUT:
+
+            dM/dE
+
+        HISTORY:
+
+            2023-05-23 - Written - Bovy (UofT)
+
+        """
+        return self._dMdE(numpy.atleast_1d(conversion.parse_energy(E, vo=self._vo)))
+
     def vmomentdensity(self, r, n, m, **kwargs):
         """
         NAME:
@@ -705,6 +731,34 @@ def _call_internal(self, *args):
         """
         return self.fE(args[0])
 
+    def _dMdE(self, E):
+        if not hasattr(self, "_rphi"):
+            self._rphi = self._setup_rphi_interpolator()
+        fE = self.fE(E)
+        out = numpy.zeros_like(E)
+        out[fE > 0.0] = (
+            16.0
+            * numpy.pi**2.0
+            * numpy.sqrt(2.0)
+            * fE[fE > 0.0]
+            * numpy.array(
+                [
+                    integrate.quad(
+                        lambda r: r**2.0
+                        * numpy.sqrt(tE - _evaluatePotentials(self._pot, r, 0.0)),
+                        0.0,
+                        self._rphi(tE),
+                    )[0]
+                    for ii, tE in enumerate(E)
+                    if fE[ii] > 0.0
+                ]
+            )
+        )
+        # Numerical issues can make the integrand's sqrt argument negative, only
+        # happens at dMdE ~ 0, so just set to zero
+        out[numpy.isnan(out)] = 0.0
+        return out
+
     def _vmomentdensity(self, r, n, m):
         if m % 2 == 1 or n % 2 == 1:
             return 0.0
@@ -776,3 +830,44 @@ def __init__(self, pot=None, denspot=None, rmax=None, scale=None, ro=None, vo=No
         sphericaldf.__init__(
             self, pot=pot, denspot=denspot, rmax=rmax, scale=scale, ro=ro, vo=vo
         )
+
+    def _dMdE(self, E):
+        if not hasattr(self, "_rphi"):
+            self._rphi = self._setup_rphi_interpolator()
+
+        def Lintegrand(t, L2lim, E):
+            return self((E, numpy.sqrt(L2lim - t**2.0)), use_physical=False)
+
+        out = (
+            16.0
+            * numpy.pi**2.0
+            * numpy.array(
+                [
+                    integrate.quad(
+                        lambda r: r
+                        * integrate.quad(
+                            Lintegrand,
+                            0.0,
+                            numpy.sqrt(
+                                2.0
+                                * r**2.0
+                                * (tE - _evaluatePotentials(self._pot, r, 0.0))
+                            ),
+                            args=(
+                                2.0
+                                * r**2.0
+                                * (tE - _evaluatePotentials(self._pot, r, 0.0)),
+                                tE,
+                            ),
+                        )[0],
+                        0.0,
+                        self._rphi(tE),
+                    )[0]
+                    for ii, tE in enumerate(E)
+                ]
+            )
+        )
+        # Numerical issues can make the integrand's sqrt argument negative, only
+        # happens at dMdE ~ 0, so just set to zero
+        out[numpy.isnan(out)] = 0.0
+        return out

galpy/potential/PowerSphericalPotentialwCutoff.py

@@ -365,15 +365,26 @@ def _mass(self, R, z=None, t=0.0):
         """
         if z is not None:
             raise AttributeError  # use general implementation
-        return (
+        R = numpy.array(R)
+        out = numpy.ones_like(R)
+        out[~numpy.isinf(R)] = (
             2.0
             * numpy.pi
-            * R ** (3.0 - self.alpha)
+            * R[~numpy.isinf(R)] ** (3.0 - self.alpha)
             / (1.5 - self.alpha / 2.0)
             * special.hyp1f1(
-                1.5 - self.alpha / 2.0, 2.5 - self.alpha / 2.0, -((R / self.rc) ** 2.0)
+                1.5 - self.alpha / 2.0,
+                2.5 - self.alpha / 2.0,
+                -((R[~numpy.isinf(R)] / self.rc) ** 2.0),
             )
         )
+        out[numpy.isinf(R)] = (
+            2.0
+            * numpy.pi
+            * self.rc ** (3.0 - self.alpha)
+            * special.gamma(1.5 - self.alpha / 2.0)
+        )
+        return out
 
     @kms_to_kpcGyrDecorator
     def _nemo_accpars(self, vo, ro):

galpy/util/conversion.py

@@ -840,6 +840,7 @@ def parse_numdens(x, ro=None, vo=None):
     "phasespacedensity2d": True,
     "phasespacedensityvelocity": True,
     "phasespacedensityvelocity2": True,
+    "energydensity": False,
     "dimensionless": False,
 }
 _voNecessary = copy.copy(_roNecessary)
@@ -851,6 +852,7 @@ def parse_numdens(x, ro=None, vo=None):
 _voNecessary["velocity2"] = True
 _voNecessary["velocity_kms"] = True
 _voNecessary["energy"] = True
+_voNecessary["energydensity"] = True
 
 
 # Determine whether or not outputs will be physical or not
@@ -1057,6 +1059,10 @@ def wrapped(*args, **kwargs):
                     fac = 1.0 / vo / ro**3.0
                     if _apy_units:
                         u = 1 / (units.km / units.s) / units.kpc**3
+                elif quantity.lower() == "energydensity":
+                    fac = 1.0 / vo**2.0
+                    if _apy_units:
+                        u = 1 / (units.km / units.s) ** 2
                 elif quantity.lower() == "dimensionless":
                     fac = 1.0
                     if _apy_units:

tests/test_quantity.py

@@ -16061,6 +16061,9 @@ def test_sphericaldf_method_returntype():
     assert isinstance(
         dfh(o.R(), o.vR(), o.vT(), o.z(), o.vz(), o.phi()), units.Quantity
     ), "sphericaldf method __call__ does not return Quantity when it should"
+    assert isinstance(
+        dfh.dMdE(o.E(pot=pot)), units.Quantity
+    ), "sphericaldf method dMdE does not return Quantity when it should"
     assert isinstance(
         dfh.vmomentdensity(1.1, 0, 0), units.Quantity
     ), "sphericaldf method vmomentdensity does not return Quantity when it should"
@@ -16121,6 +16124,12 @@ def test_sphericaldf_method_returnunit():
         raise AssertionError(
             "sphericaldf method __call__ does not return Quantity with the right units"
         )
+    try:
+        dfh.dMdE(o.E(pot=pot)).to(1 / (units.km / units.s) ** 2)
+    except units.UnitConversionError:
+        raise AssertionError(
+            "sphericaldf method dMdE does not return Quantity with the right units"
+        )
     try:
         dfh.vmomentdensity(1.1, 0, 0).to(1 / units.kpc**3)
     except units.UnitConversionError:
@@ -16213,6 +16222,13 @@ def test_sphericaldf_method_value():
         )
         < 10.0**-8.0
     ), "sphericaldf method __call__ does not return correct Quantity"
+    assert (
+        numpy.fabs(
+            dfh.dMdE(o.E(pot=pot)).to(1 / (units.km / units.s) ** 2).value
+            - dfh_nou.dMdE(o.E(pot=pot)) / vo**2
+        )
+        < 10.0**-8.0
+    ), "sphericaldf method dMdE does not return correct Quantity"
     assert (
         numpy.fabs(
             dfh.vmomentdensity(1.1, 0, 0).to(1 / units.kpc**3).value
@@ -16329,6 +16345,13 @@ def test_sphericaldf_method_inputAsQuantity():
         )
         < 10.0**-8.0
     ), "sphericaldf method __call__ does not return correct Quantity"
+    assert (
+        numpy.fabs(
+            dfh.dMdE(o.E(pot=pot)).to(1 / (units.km / units.s) ** 2).value
+            - dfh_nou.dMdE(o.E(pot=pot, use_physical=False)) / vo**2
+        )
+        < 10.0**-8.0
+    ), "sphericaldf method dMdE does not return correct Quantity"
     assert (
         numpy.fabs(
             dfh.vmomentdensity(1.1 * ro * units.kpc, 0, 0).to(1 / units.kpc**3).value