Gravitational potentials (gala.potential
)¶
Introduction¶
This subpackage provides a number of classes for working with parametric models
of gravitational potentials. There are a number of built-in potentials
implemented in C and Cython for speed and there are base classes that allow for
easy creation of new custom potential classes in
pure-Python. The potential objects have methods for computing, for example, the
potential energy, gradient, density, or mass profiles. These are particularly
useful in combination with the integrate
and dynamics
subpackages.
Also defined in this subpackage are a set of reference frames which can be used
for numerical integration of orbits in non-static reference frames. See the page
on Hamiltonian objects and reference frames for more information. Potential objects
can be combined with a reference frame and stored in a
Hamiltonian
object that provides an easy interface to
numerical orbit integration.
For code blocks below and any pages linked below, I assume the following imports have already been excuted:
>>> import astropy.units as u
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> import gala.potential as gp
>>> from gala.units import galactic, solarsystem, dimensionless
Getting started: built-in methods of potential classes¶
The built-in potentials are all initialized by passing in keyword argument
parameter values as Quantity
objects or as numeric
values in a specified unit system. To see what parameters are available for a
given potential, check the documentation for the individual classes below. You
must also specify a UnitSystem
when initializing a potential. A
unit system is a set of non-reducible units that define (at minimum) the length,
mass, time, and angle units. A few common unit systems are built in to the
package (e.g., galactic
, solarsystem
, dimensionless
).
All of the built-in potential objects have defined methods to evaluate the potential energy and the gradient/acceleration at a given position(s). For example, here we will create a potential object for a 2D point mass located at the origin with unit mass:
>>> ptmass = gp.KeplerPotential(m=1.*u.Msun, units=solarsystem)
>>> ptmass
<KeplerPotential: m=1.00 (AU,yr,solMass,rad)>
If you pass in parameters with different units, they will be converted to the specified unit system:
>>> gp.KeplerPotential(m=1047.6115*u.Mjup, units=solarsystem)
<KeplerPotential: m=1.00 (AU,yr,solMass,rad)>
If no units are specified for a parameter, it is assumed to already be
consistent with the UnitSystem
passed in:
>>> gp.KeplerPotential(m=1., units=solarsystem)
<KeplerPotential: m=1.00 (AU,yr,solMass,rad)>
The potential classes work well with the astropy.units
framework, but to
ignore units you can use the DimensionlessUnitSystem
or pass
None
as the unit system:
>>> gp.KeplerPotential(m=1., units=None)
<KeplerPotential: m=1.00 (dimensionless)>
With an instantiated potential object, we can evaluate the potential energy at some position or a set of positions. Here, we’ll evaluate a 3D point mass potential at the 3D position (1,-1,0):
>>> ptmass.energy([1.,-1.,0.]*u.au)
<Quantity [-27.91440236] AU2 / yr2>
These functions also accept both Quantity
objects or
plain ndarray
-like objects (in which case the position is
assumed to be in the unit system of the potential):
>>> ptmass.energy([1.,-1.,0.])
<Quantity [-27.91440236] AU2 / yr2>
This also works for multiple positions by passing in a 2D position (but see Conventions for a description of the interpretation of different axes):
>>> pos = np.array([[1.,-1.,0],
... [2.,3.,0]]).T
>>> ptmass.energy(pos*u.au)
<Quantity [-27.91440236, -10.94892941] AU2 / yr2>
We may also compute the gradient or acceleration:
>>> ptmass.gradient([1.,-1.,0]*u.au)
<Quantity [[ 13.95720118],
[-13.95720118],
[ 0. ]] AU / yr2>
>>> ptmass.acceleration([1.,-1.,0]*u.au)
<Quantity [[-13.95720118],
[ 13.95720118],
[ -0. ]] AU / yr2>
Some of the potential objects also have methods implemented for computing the corresponding mass density and the Hessian of the potential (matrix of 2nd derivatives) at given locations. For example:
>>> pot = gp.HernquistPotential(m=1E9*u.Msun, c=1.*u.kpc, units=galactic)
>>> pot.density([1.,-1.,0]*u.kpc)
<Quantity [7997938.82200887] solMass / kpc3>
>>> pot.hessian([1.,-1.,0]*u.kpc)
<Quantity [[[ -4.68318131e-05],
[ 5.92743432e-04],
[ 0.00000000e+00]],
[[ 5.92743432e-04],
[ -4.68318131e-05],
[ 0.00000000e+00]],
[[ 0.00000000e+00],
[ 0.00000000e+00],
[ 5.45911619e-04]]] 1 / Myr2>
Another useful method is mass_enclosed()
,
which numerically estimates the mass enclosed within a spherical shell defined
by the specified position. This numerically estimates
\(\frac{d \Phi}{d r}\) along the vector pointing at the specified position
and estimates the enclosed mass simply as
\(M(<r)\approx\frac{r^2}{G} \frac{d \Phi}{d r}\). This function can
be used to compute, for example, a mass profile:
>>> pot = gp.NFWPotential(m=1E11*u.Msun, r_s=20.*u.kpc, units=galactic)
>>> pos = np.zeros((3,100)) * u.kpc
>>> pos[0] = np.logspace(np.log10(20./100.), np.log10(20*100.), pos.shape[1]) * u.kpc
>>> m_profile = pot.mass_enclosed(pos)
>>> plt.loglog(pos[0], m_profile, marker='')
>>> plt.xlabel("$r$ [{}]".format(pos.unit.to_string(format='latex')))
>>> plt.ylabel("$M(<r)$ [{}]".format(m_profile.unit.to_string(format='latex')))
(Source code, png)
Plotting isopotentials¶
Potential objects also provide more specialized methods such as
plot_contours()
, which is a fast
way to plot either 1D slices or 2D contour plots of isopotentials. To plot a 1D
slice over the dimension of interest, pass in a grid of values for that
dimension and numerical values for the others. For example, to make a 1D plot of
the potential value as a function of \(x\) position at \(y=0, z=1\):
>>> p = gp.MiyamotoNagaiPotential(m=1E11, a=6.5, b=0.27, units=galactic)
>>> fig, ax = plt.subplots(1,1)
>>> p.plot_contours(grid=(np.linspace(-15,15,100), 0., 1.), marker='', ax=ax)
>>> E_unit = p.units['energy'] / p.units['mass']
>>> ax.set_xlabel("$x$ [{}]".format(p.units['length'].to_string(format='latex')))
>>> ax.set_ylabel("$\Phi(x,0,1)$ [{}]".format(E_unit.to_string(format='latex')))
(Source code, png)
To instead make a 2D contour plot over \(x\) and \(z\) along with
\(y=0\), pass in a 1D grid of values for \(x\) and a 1D grid of values
for \(z\) (the meshgridding is taken care of internally). Here, we choose
to draw on a pre-defined matplotlib
axes object so we can set the labels and
aspect ratio of the plot:
>>> fig,ax = plt.subplots(1, 1, figsize=(12,4))
>>> x = np.linspace(-15,15,100)
>>> z = np.linspace(-5,5,100)
>>> p.plot_contours(grid=(x, 1., z), ax=ax)
>>> ax.set_xlabel("$x$ [kpc]")
>>> ax.set_ylabel("$z$ [kpc]")
(Source code, png)
Saving / loading potential objects¶
Potential objects can be pickled
and can therefore be stored for later use. However, pickles are saved as binary
files. It may be useful to save to or load from text-based specifications of
Potential objects. This can be done with gala.potential.save()
and
gala.potential.load()
, or with the save()
and method:
>>> from gala.potential import load
>>> pot = gp.NFWPotential(m=6E11*u.Msun, r_s=20.*u.kpc,
... units=galactic)
>>> pot.save("potential.yml")
>>> load("potential.yml")
<NFWPotential: m=6.00e+11, r_s=20.00, a=1.00, b=1.00, c=1.00 (kpc,Myr,solMass,rad)>
Exporting potentials as sympy
expressions¶
Most of the potential classes can be exported to a sympy
expression that can
be used to manipulate or evaluate the form of the potential. To access this
functionality, the potential classes have a .to_sympy()
classmethod (note:
this requires that sympy
is installed):
>>> expr, vars_, pars = gp.LogarithmicPotential.to_sympy()
>>> str(expr)
'0.5*v_c**2*log(r_h**2 + z**2/q3**2 + y**2/q2**2 + x**2/q1**2)'
This method also returns a dictionary containing the coordinate variables used
in the expression as sympy
symbols, here defined as vars_
:
>>> vars_
{'x': x, 'y': y, 'z': z}
A second dictionary containing the potential parameters as sympy
symbols is
also returned, here defined as pars
:
>>> pars
{'v_c': v_c, 'r_h': r_h, 'q1': q1, 'q2': q2, 'q3': q3, 'phi': phi, 'G': G}
The expressions and variables returned can be used to perform operations on the
potential expression. For example, to create a sympy
expression for the
gradient of the potential:
>>> import sympy as sy
>>> grad = sy.derive_by_array(expr, list(vars_.values()))
>>> grad[0] # dPhi/dx
1.0*v_c**2*x/(q1**2*(r_h**2 + z**2/q3**2 + y**2/q2**2 + x**2/q1**2))
Using gala.potential¶
More details are provided in the linked pages below:
API¶
gala.potential.potential Package¶
Functions¶
|
Create a potential class from an expression for the potential. |
|
Read a potential specification file and return a |
|
Write a |
Classes¶
|
An implementation of the |
|
|
|
A baseclass for defining gravitational potentials implemented in C. |
|
A potential composed of several distinct components. |
|
Represents an N-dimensional harmonic oscillator. |
|
The Hénon-Heiles potential. |
|
Hernquist potential for a spheroid. |
|
The Isochrone potential. |
|
Jaffe potential for a spheroid. |
|
The Kepler potential for a point mass. |
|
The Kuzmin flattened disk potential. |
|
The Galactic potential used by Law and Majewski (2010) to represent the Milky Way as a three-component sum of disk, bulge, and halo. |
|
Approximation of a Triaxial NFW Potential with the flattening in the density, not the potential. |
|
Triaxial logarithmic potential. |
|
A simple, triaxial model for a galaxy bar. |
|
A simple mass-model for the Milky Way consisting of a spherical nucleus and bulge, a Miyamoto-Nagai disk, and a spherical NFW dark matter halo. |
|
Miyamoto-Nagai potential for a flattened mass distribution. |
|
General Navarro-Frenk-White potential. |
|
A null potential with 0 mass. |
|
Plummer potential for a spheroid. |
|
A baseclass for defining pure-Python gravitational potentials. |
|
A spherical power-law density profile with an exponential cutoff. |
|
Satoh potential for a flattened mass distribution. |
|
Stone potential from Stone & Ostriker (2015). |
Class Inheritance Diagram¶
gala.potential.frame.builtin Package¶
Classes¶
|
Represents a constantly rotating reference frame. |
|
Represents a static intertial reference frame. |
Class Inheritance Diagram¶
gala.potential.hamiltonian Package¶
Classes¶
|
Represents a composition of a gravitational potential and a reference frame. |