TimeInterpolatedPotential#

class gala.potential.potential.TimeInterpolatedPotential(potential_cls, time_knots, interpolation_method='cspline', *, units=None, origin=None, R=None, **kwargs)[source]#

Bases: CPotentialBase

A time-interpolated wrapper for any potential class.

This class allows any PotentialBase subclass to have time-varying parameters, origin, and rotation by interpolating between values specified at discrete time knots using GSL splines.

Parameters:

potential_clsPotentialBase subclass

The potential class to wrap.

time_knotsarray_like

Array of time values for interpolation knots. Must be monotonically increasing.

interpolation_methodstr, optional

Interpolation type. Any GSL interpolation type is supported: https://www.gnu.org/software/gsl/doc/html/interp.html Common options are: - ‘linear’: Linear interpolation - ‘cspline’: Cubic spline interpolation (default) - ‘akima’: Akima spline interpolation. This avoids unphysical wiggles in

regions where the second derivative in the underlying curve is rapidly changing, however it does not have a continuous second derivative.

‘steffen’: Steffen spline interpolation. This guarantees monotonicity of the interpolating function between the given data points. Therefore, minima and maxima can only occur exactly at the data points, and there can never be spurious oscillations between data points.

unitsUnitSystem, optional

Unit system for the potential

originarray_like, optional

Either a constant origin vector, or an array of origin vectors with shape (n_knots, n_dim).

Rarray_like, optional

Either a constant rotation matrix, or an array of rotation matrices with shape (n_knots, n_dim, n_dim).

**kwargs

Potential parameters. Each parameter can be either a constant value, or an array with shape (n_knots, *parameter_shape) for a time-varying parameter.

Examples

Create a Kepler potential with time-varying mass:

>>> import astropy.units as u
>>> from gala.potential import KeplerPotential
>>> from gala.units import galactic
>>>
>>> # Time knots in Myr
>>> times = np.linspace(0, 100, 11) * u.Myr
>>> # Mass growing linearly with time
>>> masses = np.linspace(1e10, 2e10, 11) * u.Msun
>>>
>>> pot = TimeInterpolatedPotential(
...     KeplerPotential, times, m=masses, units=galactic
... )
>>> pot.energy([1., 0, 0] * u.pc, t=0*u.Myr)
<Quantity [-44.98502151] kpc2 / Myr2>
>>> pot.energy([1., 0, 0] * u.pc, t=50*u.Myr)
<Quantity [-67.47753227] kpc2 / Myr2>

Create a potential with a time-varying rotation:

>>> # Rotation matrices for 90 degree rotation over 1 Gyr
>>> R_times = np.linspace(0, 1, 11) * u.Gyr
>>> angles = np.linspace(0, np.pi / 2, 11)
>>> Rs = np.array([R.from_rotvec([0, 0, angle]).as_matrix() for angle in angles])
>>> pot = gp.TimeInterpolatedPotential(
...     gp.LongMuraliBarPotential,
...     R_times,
...     m=1e10 * u.Msun,
...     a=3 * u.kpc,
...     b=1 * u.kpc,
...     c=0.5 * u.kpc,
...     R=Rs,
...     units=galactic,
... )
>>> pot.gradient([5., 0, 0] * u.kpc, t=0.*u.Gyr)[0, 0]
<Quantity 0.00207787 kpc / Myr2>
>>> pot.gradient([5., 0, 0] * u.kpc, t=0.5*u.Gyr)[0, 0]
<Quantity 0.0015879 kpc / Myr2>

Attributes Summary

`Wrapper`
`interpolation_method`
`ndim`
`potential_cls`
`time_knots`
`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.
`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).
`gradient`([q, t])	Compute the gradient of the gravitational potential.
`hessian`([q, t])	Compute the Hessian matrix of the gravitational potential.
`integrate_orbit`(w0[, Integrator, ...])	Integrate an orbit in the current potential using the integrator class provided.
`mass_enclosed`(q, t)	Estimate the mass enclosed within the given position by assuming the potential is spherical.
`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.
`replace_units`(units)	Change the unit system of this potential.
`replicate`(**kwargs)	Create a copy of this potential with possibly different parameters.
`save`(f)	Save the potential to a text file.
`to_latex`()	Return a string LaTeX representation of this potential
`to_sympy`()	Return a representation of this potential class as a sympy expression

Attributes Documentation

Wrapper = None#

interpolation_method = <PotentialParameter: interpolation_method>#

ndim = 3#

potential_cls = <PotentialParameter: potential_cls>#

time_knots = <PotentialParameter: time_knots [time]>#

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.

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

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().