Interpolated spherical potential

The interpSphericalPotential class provides a general interface to generate interpolated instances of spherical potentials or lists of such potentials. This interpolated potential can be used in any function where other three-dimensional galpy potentials can be used. This includes functions that use C to speed up calculations.

The interpSphericalPotential interpolates the radial force of a spherical potential and determines the potential and its second derivative from the base radial-force interpolation object. To set up an interpSphericalPotential instance, either provide it with a function that returns the radial force or with a galpy potential or list of potentials, and also provide the radial interpolation grid in each case.

For example, to use a function that gives the radial force, do

>>> from galpy import potential
>>> ip= potential.interpSphericalPotential(rforce=lambda r: -1./r,
                        rgrid=numpy.geomspace(0.01,20,101),Phi0=0.)

which sets up an interpSphericalPotential instance that has the same radial force as the spherical LogarithmicHaloPotential. If you have a function that gives the enclosed mass within a given radius, simply pass it divided by \(-r^2\) to set up a interpSphericalPotential instance for this enclosed-mass profile. Note that the force function has to return the force in galpy’s internal units and it has to take the radius in internal units. For example, if you have the enclosed mass in solar masses, divide it by the following mass_conversion factor

>>> from galpy.util import conversion
>>> mass_conversion= conversion.mass_in_msol(vo,ro)

where vo and ro are the usual unit-conversion parameters (they cannot be Quantities in this case, they need to be floats in km/s and kpc). To convert the radius in kpc to/from internal units, simply divide/multiply by ro. The radial interpolation grid also specifies radii in internal units.

Alternatively, you can specify a galpy potential or list of potentials and (again) the radial interpolation grid, as for example,

>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> ip= potential.interpSphericalPotential(rforce=lp,
                      rgrid=numpy.geomspace(0.01,20,101))

Note that, because the potential is defined through integration of the (negative) radial force, we need to specify the potential at the smallest grid point, which is done through the Phi0= keyword in the first example. When using a galpy potential (or list), this value is automatically determined.

Also note that the density of the potential is assumed to be zero outside of the final radial grid point. That is, the potential outside of the final grid point is \(-GM/r\) where \(M\) is the mass within the final grid point. If during an orbit integration, the orbit strays outside of the interpolation grid, a warning is issued.

Warning

The density of a interpSphericalPotential instance is assumed to be zero outside of the largest radial grid point.

class galpy.potential.interpSphericalPotential(self, rforce=None, rgrid=numpy.geomspace(0.01, 20, 101), Phi0=None, ro=None, vo=None)[source]

Class that interpolates a spherical potential on a grid

__init__(rforce=None, rgrid=array([1.00000000e-02, 1.07897231e-02, 1.16418125e-02, 1.25611933e-02, 1.35531798e-02, 1.46235057e-02, 1.57783578e-02, 1.70244112e-02, 1.83688683e-02, 1.98195003e-02, 2.13846920e-02, 2.30734906e-02, 2.48956574e-02, 2.68617250e-02, 2.89830576e-02, 3.12719166e-02, 3.37415321e-02, 3.64061789e-02, 3.92812590e-02, 4.23833909e-02, 4.57305052e-02, 4.93419489e-02, 5.32385966e-02, 5.74429717e-02, 6.19793759e-02, 6.68740305e-02, 7.21552273e-02, 7.78534923e-02, 8.40017626e-02, 9.06355759e-02, 9.77932769e-02, 1.05516238e-01, 1.13849099e-01, 1.22840026e-01, 1.32540986e-01, 1.43008054e-01, 1.54301731e-01, 1.66487295e-01, 1.79635182e-01, 1.93821388e-01, 2.09127911e-01, 2.25643225e-01, 2.43462792e-01, 2.62689611e-01, 2.83434817e-01, 3.05818320e-01, 3.29969499e-01, 3.56027954e-01, 3.84144304e-01, 4.14481068e-01, 4.47213595e-01, 4.82531087e-01, 5.20637682e-01, 5.61753643e-01, 6.06116627e-01, 6.53983058e-01, 7.05629612e-01, 7.61354813e-01, 8.21480762e-01, 8.86354997e-01, 9.56352500e-01, 1.03187787e+00, 1.11336765e+00, 1.20129286e+00, 1.29616174e+00, 1.39852263e+00, 1.50896719e+00, 1.62813382e+00, 1.75671131e+00, 1.89544286e+00, 2.04513037e+00, 2.20663904e+00, 2.38090242e+00, 2.56892779e+00, 2.77180196e+00, 2.99069756e+00, 3.22687986e+00, 3.48171402e+00, 3.75667303e+00, 4.05334618e+00, 4.37344830e+00, 4.71882962e+00, 5.09148650e+00, 5.49357296e+00, 5.92741311e+00, 6.39551462e+00, 6.90058320e+00, 7.44553820e+00, 8.03352956e+00, 8.66795596e+00, 9.35248448e+00, 1.00910718e+01, 1.08879871e+01, 1.17478366e+01, 1.26755904e+01, 1.36766110e+01, 1.47566846e+01, 1.59220541e+01, 1.71794555e+01, 1.85361568e+01, 2.00000000e+01]), Phi0=None, ro=None, vo=None)[source]

Initialize an interpolated, spherical potential.

Parameters:

  • rforce (function or galpy Potential instance or list thereof, optional) – Either a function that gives the radial force (in internal units) as a function of r (in internal units) or a galpy Potential instance or list thereof. The default is None.

  • rgrid (numpy.ndarray, optional) – Radial grid in internal units on which to evaluate the potential for interpolation (note that beyond rgrid[-1], the potential is extrapolated as -GM(<rgrid[-1])/r). The default is numpy.geomspace(0.01,20,101).

  • Phi0 (float, optional) – Value of the potential at rgrid[0] in internal units (only necessary when rforce is a function, for galpy potentials automatically determined). The default is None.

  • ro (float or Quantity, optional) – Distance scale for translation into internal units (default from configuration file).

  • vo (float or Quantity, optional) – Velocity scale for translation into internal units (default from configuration file).

Notes

  • 2020-07-13 - Written - Bovy (UofT)