MilkyWayPotential2022#

class gala.potential.potential.MilkyWayPotential2022(units=<UnitSystem (kpc, Myr, solMass, rad)>, disk=None, halo=None, bulge=None, nucleus=None)[source]#

Bases: CCompositePotential

A mass-model for the Milky Way consisting of a spherical nucleus and bulge, a 3-component sum of Miyamoto-Nagai disks to represent an exponential disk, and a spherical NFW dark matter halo.

The disk model is fit to the Eilers et al. 2019 rotation curve for the radial dependence, and the shape of the phase-space spiral in the solar neighborhood is used to set the vertical structure in Darragh-Ford et al. 2023.

Other parameters are fixed by fitting to a compilation of recent mass measurements of the Milky Way, from 10 pc to ~150 kpc.

Parameters:

unitsUnitSystem (optional): Set of non-reducable units that specify (at minimum) the length, mass, time, and angle units.
diskdict (optional): Parameters to be passed to the MN3ExponentialDiskPotential.
bulgedict (optional): Parameters to be passed to the HernquistPotential.
halodict (optional): Parameters to be passed to the NFWPotential.
nucleusdict (optional): Parameters to be passed to the HernquistPotential.
Note: in subclassing, order of arguments must match order of potential
components added at bottom of init.

Attributes Summary

`Wrapper`
`ndim`
`parameters`
`units`

Methods Summary

`__call__`(q)	Call self as a function.
`acceleration`([q, t])	Compute the gravitational acceleration at the given position(s).
`as_interop`(package, **kwargs)	Interoperability with other Galactic dynamics packages
`circular_velocity`([q, t])	Estimate the circular velocity at given position(s) assuming spherical symmetry.
`clear`()
`copy`()
`density`([q, t])	Compute the mass density at the given position(s).
`energy`([q, t])	Compute the gravitational potential energy at the given position(s).
`fromkeys`(/, iterable[, value])	Create a new ordered dictionary with keys from iterable and values set to value.
`get`(key[, default])	Return the value for key if key is in the dictionary, else default.
`gradient`([q, t])	Compute the gradient of the gravitational potential.
`hessian`([q, t])	Compute the Hessian matrix of the gravitational potential.
`integrate_orbit`(args, *kwargs)	Integrate an orbit in the current potential using the integrator class provided.
`items`()
`keys`()
`mass_enclosed`(q, t)	Estimate the mass enclosed within the given position by assuming the potential is spherical.
`move_to_end`(/, key[, last])	Move an existing element to the end (or beginning if last is false).
`plot_contours`(grid[, t, filled, ax, labels, ...])	Plot equipotentials contours.
`plot_density_contours`(grid[, t, filled, ax, ...])	Plot density contours.
`plot_rotation_curve`(R_grid[, t, ax, labels])	Plot the rotation curve or circular velocity curve for this potential on the input grid of cylindrical radii.
`pop`(/, key[, default])	If the key is not found, return the default if given; otherwise, raise a KeyError.
`popitem`(/[, last])	Remove and return a (key, value) pair from the dictionary.
`replace_units`(units)	Change the unit system of this potential.
`replicate`(**kwargs)	Return a copy of the potential instance with some parameter values changed.
`save`(f)	Save the potential to a text file.
`setdefault`(/, key[, default])	Insert key with a value of default if key is not in the dictionary.
`to_latex`()	Return a string LaTeX representation of this potential
`to_sympy`()	Return a representation of this potential class as a sympy expression
`update`([E, ]**F)	If E is present and has a .keys() method, then does: for k in E.keys(): D[k] = E[k] If E is present and lacks a .keys() method, then does: for k, v in E: D[k] = v In either case, this is followed by: for k in F: D[k] = F[k]
`values`()

Attributes Documentation

Wrapper = None#

ndim = 3#

parameters#

units#

Methods Documentation

__call__(q)#: Call self as a function.

acceleration(q=None, t=0.0, **coord_kwargs)#

Compute the gravitational acceleration at the given position(s).

The gravitational acceleration is the negative gradient of the potential: \(\vec{a} = -\nabla \phi\). This is the acceleration experienced by a test particle in the gravitational field.

Parameters:

qPhaseSpacePosition, Quantity, array_like, optional: Position(s) at which to compute the gravitational acceleration. If the input has no units (i.e., is an ndarray), it is assumed to be in the same unit system as the potential. If using symmetry coordinates, pass those as keyword arguments instead and leave q as None.
tnumeric, Quantity, optional: Time at which to evaluate the acceleration. Default is 0.
**coord_kwargs: For potentials with spherical or cylindrical symmetry, you can optionally provide coordinates in the natural coordinate system. For spherical potentials, use r=.... For cylindrical potentials, use R=... and optionally z=... (defaults to 0).

Returns:

accQuantity: The gravitational acceleration vector(s). Has the same shape as the input position array q. Units are acceleration (e.g., m/s² in SI units).

See also

gradient: Compute the potential gradient (negative acceleration).

Notes

This method is equivalent to -self.gradient(q, t) and is provided for convenience in orbital integration and dynamics calculations.

as_interop(package, **kwargs)#

Interoperability with other Galactic dynamics packages

Parameters:

packagestr: The package to export the potential to. Currently supported packages are "galpy" and "agama".
kwargs: Any additional keyword arguments are passed to the interop function.

circular_velocity(q=None, t=0.0, **coord_kwargs)#

Estimate the circular velocity at given position(s) assuming spherical symmetry.

The circular velocity is the speed required for a circular orbit at the given radius in a spherically symmetric potential. It is computed using \(v_{\rm circ}(r) = \sqrt{r |dΦ/dr|}\), where the radial derivative is evaluated from the potential gradient.

Parameters:

qPhaseSpacePosition, Quantity, array_like, optional: Position(s) at which to estimate the circular velocity. The calculation uses the spherical radius from the origin. If the input has no units, it is assumed to be in the same unit system as the potential. If using symmetry coordinates, pass those as keyword arguments instead and leave q as None.
tnumeric, Quantity, optional: Time at which to evaluate the circular velocity. Default is 0.
**coord_kwargs: For potentials with spherical or cylindrical symmetry, you can optionally provide coordinates in the natural coordinate system. For spherical potentials, use r=.... For cylindrical potentials, use R=... and optionally z=... (defaults to 0).

Returns:

vcircQuantity: Circular velocity at the spherical radius corresponding to each position. For input shape (n_dim, n_positions), returns shape (n_positions,). Units are velocity (e.g., m/s in SI units).

Notes

This method assumes the potential is approximately spherically symmetric. The circular velocity is computed using the relation:

\[v_{\rm circ}(r) = \sqrt{r \left| \frac{d\Phi}{dr} \right|}\]

where the radial derivative is computed from the Cartesian gradient via \(dΦ/dr = \vec{\nabla}Φ \cdot \hat{r}\).

For exactly spherical potentials, this gives the speed of circular orbits. For non-spherical potentials, this provides an approximation useful for initial orbit estimates.

clear() → None. Remove all items from od.#

copy() → a shallow copy of od#

density(q=None, t=0.0, **coord_kwargs)#

Compute the mass density at the given position(s).

For potentials that have an associated mass distribution, this method computes the mass density rho(q, t) at the specified positions and time. The density is related to the potential via Poisson’s equation: \(\nabla^2 \phi = 4\pi G \rho\).

Parameters:

qPhaseSpacePosition, Quantity, array_like, optional: Position(s) at which to evaluate the mass density. If the input has no units (i.e., is an ndarray), it is assumed to be in the same unit system as the potential. Shape should be (n_dim,) for a single position or (n_dim, n_positions) for multiple positions. If using symmetry coordinates, pass those as keyword arguments instead and leave q as None.
tnumeric, Quantity, optional: Time at which to evaluate the mass density. Default is 0.
**coord_kwargs: For potentials with spherical or cylindrical symmetry, you can optionally provide coordinates in the natural coordinate system. For spherical potentials, use r=.... For cylindrical potentials, use R=... and optionally z=... (defaults to 0).

Returns:

densQuantity: The mass density at the specified position(s). For input shape (n_dim, n_positions), returns shape (n_positions,). Units are mass density (e.g., kg/m³ in SI units).

Raises:

NotImplementedError: If the potential does not have an implemented density function.

Notes

Not all potential models have an implemented density function. For potentials without a density implementation, this method will raise a NotImplementedError.

The density is computed using the relationship with the potential’s Laplacian (when available) or from the underlying mass model.

energy(q=None, t=0.0, **coord_kwargs)#

Compute the gravitational potential energy at the given position(s).

The potential energy per unit mass is evaluated at the specified position(s) and time.

Parameters:

qPhaseSpacePosition, Quantity, array_like, optional: Position(s) at which to evaluate the potential. If the input has no units (i.e., is an ndarray), it is assumed to be in the same unit system as the potential. Shape should be (n_dim,) for a single position or (n_dim, n_positions) for multiple positions. If using symmetry coordinates, pass those as keyword arguments instead and leave q as None.
tnumeric, Quantity, optional: Time at which to evaluate the potential. Default is 0.
**coord_kwargs: For potentials with spherical or cylindrical symmetry, you can optionally provide coordinates in the natural coordinate system. For spherical potentials, use r=.... For cylindrical potentials, use R=... and optionally z=... (defaults to 0).

Returns:

EQuantity: The gravitational potential energy per unit mass. For input shape (n_dim, n_positions), returns shape (n_positions,). Units are specific energy (e.g., m²/s² in SI units).

Notes

The potential energy is related to the gravitational acceleration by \(\vec{a} = -\nabla \phi\), where φ is the potential energy per unit mass.

Examples

Using Cartesian coordinates (works for all potentials):

>>> import astropy.units as u
>>> import numpy as np
>>> pot = SomePotential(...)
>>> xyz = np.array([[1., 0., 0.]]).T * u.kpc
>>> pot.energy(xyz)

For spherical potentials, you can use spherical radius:

>>> pot = HernquistPotential(m=1e10*u.Msun, c=5*u.kpc)
>>> r = np.linspace(0.1, 10, 100) * u.kpc
>>> pot.energy(r=r)

For cylindrical potentials, you can use R and z:

>>> pot = MiyamotoNagaiPotential(m=1e11*u.Msun, a=3*u.kpc, b=0.3*u.kpc)
>>> R = np.linspace(1, 15, 100) * u.kpc
>>> pot.energy(R=R, z=0*u.kpc)
>>> pot.energy(R=R)  # z defaults to 0

classmethod fromkeys(/, iterable, value=None)#: Create a new ordered dictionary with keys from iterable and values set to value.

get(key, default=None, /)#: Return the value for key if key is in the dictionary, else default.

gradient(q=None, t=0.0, **coord_kwargs)#

Compute the gradient of the gravitational potential.

Parameters:

qPhaseSpacePosition, Quantity, array_like, optional: Position(s) at which to evaluate the potential gradient. If the input has no units (i.e., is an ndarray), it is assumed to be in the same unit system as the potential. Shape should be (n_dim,) for a single position or (n_dim, n_positions) for multiple positions. If using symmetry coordinates, pass those as keyword arguments instead and leave q as None.
tnumeric, Quantity, optional: Time at which to evaluate the potential gradient. Default is 0.
**coord_kwargs: For potentials with spherical or cylindrical symmetry, you can optionally provide coordinates in the natural coordinate system. For spherical potentials, use r=.... For cylindrical potentials, use R=... and optionally z=... (defaults to 0).

Returns:

gradQuantity: The gradient of the gravitational potential. Has the same shape as the input position array q. Units are acceleration (e.g., m/s² in SI units). To get gravitational acceleration, use acceleration() or compute -gradient().