galpy: A python LIBRARY FOR GALACTIC DYNAMICS

The overall structure of the galpy package is illustrated in Figure 1, where the package and its main subpackages are imported in a typical python session. These subpackages have fairly obvious names: galpy.potential contains classes and functions related to gravitational potentials, galpy.orbit consists of the Orbit class and all of the functionality related to it, galpy.actionAngle has specialized classes and functions for dealing with action–angle coordinates, and galpy.df is composed of classes and functions related to DFs. The galpy.util subpackage collects a large variety of utility functions related to plotting, coordinate transformations, and unit conversions.

Zoom In Zoom Out Reset image size

As described in more detail in Section 3, specific potentials are included as sub-classes of a general galpy.potential.Potential class under galpy.potential. For example, galpy.potential contains a Navarro–Frenk–White potential (NFW; Navarro et al. 1997), which can be imported and used as

This initialization is explained further below. Potential objects form the basis of almost all of the galpy functionality.

The Orbit object is another central part of galpy's infrastructure. It can be set up in a variety of flexible ways that are described in Section 4. For the general overview, we note that Orbit instances can be initialized as

which sets up an Orbit instance o that represents the initial conditions of an orbit. This orbit can then be integrated and its orbital characteristics can be evaluated using member functions.

The galpy.actionAngle module contains classes that represent different ways of calculating action–angle coordinates. These are always set up for a specific potential and can then be used to calculate actions, frequencies, and angles for specific orbits. A quick example that illustrates this using the NFW potential and Orbit instance defined above is

where the output are the radial, azimuthal, and vertical action of the orbit. The galpy.actionAngle module is described further in Section 5.

The galpy.df subpackage contains a few classes that represent different DFs. These typically depend on Potential and Orbit instances and in some cases they use actionAngle instances when the DF depends on the actions. For example, one action-based DF that is included in galpy is the quasi-isothermal DF (Binney 2010; Binney & McMillan 2011), which can be used as follows:

where the output is the value of the DF for this orbit. The various DF classes currently included in galpy are discussed in Section 6.

Finally, the utility subpackage galpy.util contains plotting utilities in galpy.util.bovy_plot, coordinate-transformation routines in galpy.util.bovy_coords, and unit-conversion functions in galpy.util.bovy_conversion. The latter are presented in the next subsection, as they relate to the system of units used preferentially by galpy. The coordinate transformations are those between equatorial and Galactic coordinates, Galactic coordinates and Galactocentric coordinates, and uncertainty propagation between some of these coordinate systems. These are discussed in some detail in the Appendix. The plotting routines are simple wrappers of matplotlib plotting functions; these are documented at http://galpy.readthedocs.org/en/latest/reference/bovyplot.html.

While generally galpy does not care about the units of its inputs (as long as they are consistent), some of the functionality will work best when using natural units. These are normalized units where the circular velocity is one at a cylindrical radius of one and height zero. Positions are therefore scaled to the physical radius that is defined to be "1" and velocities are scaled to the physical velocity at that radius. For example, for the Milky Way a reasonable choice is to normalize positions by 8 kpc and velocities by 220 km s⁻¹ in a model where the Sun is assumed to be at 8 kpc from the Galactic center and the circular velocity at the Sun is 220 km s⁻¹ (e.g., Bovy et al. 2012). Potential instances can easily be set up in this system of units by using the normalize= keyword when initializing (see below).

Functions to translate between physical units and the natural units defined above are contained in the galpy.util.bovy_conversion module. This module only contains non-trivial conversions; conversions of positions, velocities, energies, and actions are straightforward and therefore not included. Some of the conversion functions are illustrated in Figure 2. Further conversions to alternative units are also included: For example, force_in_10m13kms2 to convert forces to 10⁻¹³ km s⁻²; dens_in_gevcc to convert densities to GeV cm⁻³ (useful for dark-matter detection work); and dens_in_meanmatterdens to express densities in units of the mean-matter density in the universe (useful for dealing with virial quantities).

Zoom In Zoom Out Reset image size

As mentioned before, galpy will work in different units than natural units as well, as long as the system of units is consistent (e.g., the forces have to be in units consistent with the positions and velocities for the orbit integration to work). Natural units are only really depended on in situations where numerical determinations of zeros of a function or of integrals are difficult and approximations are made. An example is the evaluation of actions and frequencies, which often involve integrals that can be difficult to evaluate for close-to-circular or very eccentric orbits. In the code, approximations are made that depend in some cases on quantities being expressed in natural units. Even if most galpy functions work well in any system of units, because galpy is being developed and tested primarily in natural units, it is recommended to use these units when using galpy.

All three-dimensional Potential classes inherit from the general Potential class. (Two-dimensional Potential classes inherit from a general planarPotential class and one-dimensional classes inherit from a general linearPotential class). Typically, library-use of any Potential instance is through methods that are defined in the general Potential class. For basic functionality such as evaluating the potential or its forces at a point, this general Potential routine multiplies the subclass routines with an amplitude parameter; for more complicated routines such as calculating the circular velocity or the location of Lindblad frequencies, more work is done by the general class. The philosophy behind this implementation is that it allows users to implement additional potentials with a minimum of effort, only requiring the basic properties of a potential to be implemented. Allowing the general Potential class to multiply in an amplitude is helpful for using the natural units described above.

A basic implementation of a new specific potential in galpy requires the user to write a class that inherits from the general galpy.potential.Potential class and that has methods _evaluate, _Rforce, and _zforce. Any new class should take an amp parameter upon initialization and pass this to the general Potential class' initializer, by calling galpy.potential.Potential.__init__(self,amp=amp) within the new class' __init__ function. More advanced functionality or non-axisymmetric potentials require additional functions to be implemented specifying the second derivatives and derivatives with respect to azimuth. All of the pure-python functionality of galpy can then be applied to the new potential. While the forces and second derivatives could be obtained by taking numerical derivatives of the potential, such that a new potential could be defined using only the _evaluate function, this is not currently supported. The test suite described in Section 8.2 automatically checks for new potentials defined at the top-level of galpy that the forces are the negative derivative of the potential and that the second derivatives are correct, by using numerical derivatives. The density can be explicitly implemented using the _dens method; if it is not included then the density is computed by using the Poisson equation if all of the relevant second derivatives are implemented.

Many of the methods available for any Potential instance are illustrated in Figure 3, using a double-exponential disk potential as an example. In this example, we use the normalize=1. keyword to normalize the potential such that it has a circular velocity of 1 at R = 1 and z = 0, the preferred set of units in galpy (see Section 2.2). The normalize=X keyword in general normalizes the potential such that minus the radial force at R = 1 is equal to X. This can be used to normalize a sum of potentials to have a circular velocity of 1 at R = 1, by making sure that the individual potentials i's X_i sum to 1.

Zoom In Zoom Out Reset image size

I give an example of a new Potential class in Figure 4. This class implements a potential that smoothly interpolates between two given potentials as a function of time. We will use it in Section 5 to look at changes to orbits when the potential is changed adiabatically. Once this new Potential class is defined, it can be used like any other potential (note that use of some methods requires the second derivatives of the potential to be implemented as well).

Zoom In Zoom Out Reset image size

Almost all functions in galpy that take a Potential instance as an argument can also take lists of these. That is, new potentials can be specified using arbitrary combinations of basic potentials. When the potential is used, the contribution from each constituent potential is added to the total potential, force, etc. The galpy.potential module contains functions evaluatePotentials, evaluateRforces, evaluatephiforces, evaluatezforces, as well as second derivatives evaluateR2derivs, evaluatez2derivs, evaluateRzderivs, and evaluateDensities that can be used to evaluate Potential functions for lists of potentials. Similar functions exist for the other Potential methods used in Figure 3 (e.g., vcirc). New functionality that makes use of Potential objects should be written using these functions to evaluate the potential and its properties to sustain the support for lists of Potential objects. That is, one should make use of functions that can act on lists of Potential instances, such as evaluateRforces, rather than instance methods, such as Rforce.

In certain parts of the codebase, galpy uses C to speed up computations. Currently, this is limited to orbit integration and the calculation of certain types of action–angle coordinates. To make use of this functionality, Potential classes implemented in python need to also be explicitly implemented in C by the user and this needs to be registered by setting the hasC attribute of the python Potential subclass in question to True. This procedure is illustrated in Figure 5 in the case of LogarithmicHaloPotential, a simple logarithmic potential contained in galpy. This minimal implementation only implements the potential itself and the force components; glue is written to make it possible to call this function from within python. By following the procedure given in Figure 5, new, user-contributed potentials will be automatically incorporated into the galpy framework: wherever the code automatically switches to C routines to speed up the code, this will be done for the new Potential subclass as well, whether it is used on its own or whether it is contained in a list of Potential instances, all of which have their own C implementation.

Zoom In Zoom Out Reset image size

Full details of the procedure for adding C implementations can be found in the online documentation at http://galpy.readthedocs.org/en/latest/potential.html#adding-potentials-to-the-galpy-framework.

All potentials in galpy should have python implementations as many functions only use the python Potential methods; the C implementations are only an additional feature for potentials to allow certain computations to be sped up.

galpy contains a large number of axisymmetric, three-dimensional potentials that can be combined to form a realistic model for a galaxy. Many of the methods that are defined for each potential are illustrated in Figure 3. The rotation curves for a subset of the axisymmetric potentials are shown in Figure 6. Some of the potentials shown in this figure are implemented as special cases of more general potential classes. This is the case in particular for the KeplerPotential, which is a special case of general PowerSphericalPotentials with densities ρ∝r^−α; the $\rho \propto r^{-1}\,e^{-r^2}$, which is a special case of the PowerSphericalPotentialwCutoff potential with density ρ∝r^−α exp (− (r/rc)²); and the triplet of Hernquist, Jaffe, and NFW potentials, which are subclasses of a general TwoPowerSphericalPotential with density ρ∝(r/a)^−α (1 + r/a)^{α − β}.

Zoom In Zoom Out Reset image size

Because some of the axisymmetric potentials' evaluation of the potential itself or of the derivatives is computationally expensive, a general class interpRZPotential is provided that tabulates the potential and its forces on a grid for a given axisymmetric potential and uses spline interpolation to evaluate the potential. This is implemented simply by interpolating separately the potential itself, the forces and second derivatives, and the circular velocity curve, its derivative, and the various frequencies (epicycle and vertical). That is, no effort is made to calculate the forces using derivatives of the splines interpolating the potentials, etc. Therefore care should be taken to build a dense enough interpolation grid to avoid large inconsistencies between the potential and its forces and derivatives.

Interpolated potentials implemented through interpRZPotential can be used wherever Potential instances can be used anywhere in galpy. This includes functions that use C to speed up computations; to use the latter feature one needs to specify enable_c=True when setting up the interpRZPotential instance, which stores the potential and force tabulations in a way in which it can be used in C and sets the hasC attribute (see Section 3.2) to True to register that the instance has a C interface.

galpy contains a number of non-axisymmetric potentials that can be used to investigate the effect of non-axisymmetry on the dynamical structure of galaxies. Currently, all of these are two-dimensional potentials with a single exception. Primarily, these potentials represent simple parametric forms to represent perturbations to the potential of disk galaxies. galpy contains a general cos mϕ perturbation as CosmphiDiskPotential with potential Φ(R, ϕ)∝R^p cos (m (ϕ − ϕ_b)). The special cases corresponding to a lopsided disk (LopsidedDiskPotential; m = 1) and an elliptical disk (EllipticalDiskPotential; m = 2; Kuijken & Tremaine 1994) are also included; these two are also implemented in C, while the more general version currently is not. These models can all be turned on smoothly in time by specifying a time tform when the perturbation starts growing and a time period tsteady over which it is fully grown; the growth function is that of Dehnen (2000).

galpy also contains models for the time-dependent perturbations coming from spiral structure and the bar. For spiral structure this is a simple logarithmic spiral model Φ(R, ϕ, t)∝cos (α ln R − m (ϕ − Ω_st − γ)) that can either be steady (SteadyLogSpiralPotential) or transient (TransientLogSpiralPotential). The steady model can be turned on slowly in a similar way as the cos mϕ perturbations above, to allow for adiabatic growth of the spiral. The transient spiral has an amplitude ∝exp (− [t − t₀]²/2σ²). The bar potential that is contained in galpy is the simple rotating quadrupole of Dehnen (2000): DehnenBarPotential.

The one non-axisymmetric potential that is not limited to two dimensions is the potential corresponding to a moving object, MovingObjectPotential. This potential is initialized by giving an integrated Orbit instance (see below) and a mass. This potential can then be added to the potential that the moving object was integrated in to include the gravity from the moving object on test particles. This can be used, for example, to simulate the impact of molecular clouds on stellar orbits in the disk.

The galpy.potential module also contains a model MWPotential2014 for the Milky Way's gravitational potential that is designed to provide a simple, easy-to-use model in cases where a realistic model for the Milky Way is required. This model is fit to some of the dynamical data on the Milky Way as described in this section. However, it is important to keep in mind that this model is merely provided for convenience and that it is not designed to provide the best possible current model for the Milky Way's gravitational forces. Any application that is significantly affected by the current uncertainty in the potential should use the various tools in the galpy library to explore the dependence on different Milky Way potential parameters and functional forms.

We perform a fit that is very similar to that described in Section 5 of Bovy & Rix (2013), but uses some additional data to provide a more realistic description of the Milky Way on small and large scales. The potential model is a simplified version of that in Bovy & Rix (2013) and consists of a bulge modeled as a power-law density profile that is exponentially cut-off (PowerSphericalPotentialwCutoff) with a power-law exponent of −1.8 and a cut-off radius of 1.9 kpc; a MiyamotoNagaiPotential disk; and a dark-matter halo described by an NFWPotential. The relative amplitudes of these three components, the scale length and height of the disk potential, and the scale radius of the NFW halo are fit to the following data:

• 1.  
  the velocity dispersion σ_b = 117 ± 15 km s⁻¹ in Baade's window (Dehnen & Binney 1998; Binney & Tremaine 2008);
• 2.  
  the vertical force at the solar circle at 1.1 kpc from the plane |F_Z| = 67 ± 6 (2πG M_☉ pc⁻²) and the local visible surface density Σ = 55 ± 5 M_☉ pc⁻² from Zhang et al. (2013);
• 3.  
  the vertical force measurements at 1.1 kpc from the plane of Bovy & Rix (2013);
• 4.  
  the terminal-velocity measurements of Clemens (1985) and McClure-Griffiths & Dickey (2007), modeled in the same way as in Bovy & Rix (2013);
• 5.  
  the measurement of the mid-plane density at the solar circle of Holmberg & Flynn (2000): ρ(R₀, Z = 0) = 0.10 ± 0.01 M_☉ pc⁻³;
• 6.  
  the constraint on the logarithmic slope of the rotation curve from Bovy et al. (2012), represented in the same way as in Equation (41) in Bovy & Rix (2013);
• 7.  
  the measurement of the total mass within 60 kpc from Xue et al. (2008): M(r < 60 kpc) = 4.0 ± 0.7 × 10¹¹ M_☉.

In addition to these constraints, we set the solar distance to the Galactic center to R₀ = 8 kpc and the circular velocity at the Sun to V₀ = 220 km s⁻¹ (Bovy et al. 2012). We fit the potential model to these data and then adjust the parameters to a close, simple number. The parameters of the resulting MWPotential2014 are shown in Table 1.

Table 1. Parameters and Properties of MWPotential2014

Parameter	MWPotential2014	Constraint
R0 ( kpc)	8	Fixed
vc(R0) ( km s−1)	220	Fixed
fb	0.05	⋅⋅⋅
fd	0.60	⋅⋅⋅
fh	0.35	⋅⋅⋅
$\mathrm{Bulge\ power{\rm -}law\ exponent}$	−1.8	Fixed
$\mathrm{Bulge\ cut{\rm -}off\ radius}\,(\,\mathrm{kpc})$	1.9	Fixed
a ( kpc)	3.0	⋅⋅⋅
b ( pc)	280	⋅⋅⋅
Halo rs ( kpc)	16	⋅⋅⋅
σb ( km s−1)	109	117 ± 15
FZ(R0, 1.1 kpc) (2πG M☉ pc−2)	72	67 ± 6
Σvis(R0) (M☉ pc−2)	53	55 ± 5
FZ scale length ( kpc)	3.2	2.7 ± 0.1
ρ(R0, z = 0) (M☉ pc−3)	0.10	0.10 ± 0.01
$({d}\ln v_c/{d}\ln R)|_{R_0}$	−0.10	−0.2 to 0
M(r < 60 kpc) (1011 M☉)	4.1	4.0 ± 0.7
Mb (1010 M☉)	0.5	⋅⋅⋅
Md (1010 M☉)	6.8	⋅⋅⋅
Rd ( kpc)	2.6	⋅⋅⋅
ρDM(R0) (M☉ pc−3)	0.008	⋅⋅⋅
Mvir (1012 M☉)	0.8	⋅⋅⋅
rvir ( kpc)	245	⋅⋅⋅
Concentration	15.3	⋅⋅⋅
vesc(R0) ( km s−1)	513	⋅⋅⋅

Notes. The top part of this table gives the explicit parameters of the MWPotential2014 model: the relative contribution from the bulge (f_b), disk (f_d), and halo (f_h) to the radial force at R₀; the power-law exponent and the exponential cut-off radius of the bulge density; the scale length a and scale height b of the Miyamoto–Nagai disk model; and the scale radius r_s of the NFW halo model. The second section displays MWPotential2014's values for some of the constraints that the parameters of the model are fit to. The final part of this table gives some additional properties of the potential.

Download table as:  ASCIITypeset image

We compare the properties of MWPotential2014 to the data that it was fit to in the second part of Table 1 and in Figure 7. The "F_Z scale length" that is shown in Table 1 is an approximate exponential scale length of the vertical force at 1.1 kpc between 4 and 9 kpc, which is compared to the scale length obtained from an exponential fit to the data of Bovy & Rix (2013), which are also displayed in Figure 7.

Zoom In Zoom Out Reset image size

Other properties of MWPotential2014 are shown in Figures 8 and 9 and in the third part of Table 1. Figure 8 displays the rotation curve and its decomposition into bulge, disk, and dark-halo contributions. The masses of the three components are given in Table 1; the virial properties (radius, mass, and concentration) of the NFW halo are also presented there. R_d is an approximate exponential scale length from a fit to the Miyamoto–Nagai-disk density, which is compared to the determination of R_d = 2.15 ± 0.14 kpc of Bovy & Rix (2013). The local dark-mater density ρ_DM = 0.008 M_☉ pc⁻³ is in good agreement with current constraints (e.g., Bovy & Tremaine 2012; Zhang et al. 2013; Piffl et al. 2014) as is the escape velocity (e.g., Smith et al. 2007). Figure 9 shows contours of the potential and the density of MWPotential2014 in the region close to the disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

MWPotential2014 does not contain the supermassive black hole at the center of the Milky Way. If the black hole's gravity needs to be included, e.g., for tracing the orbits of stars or pulsars kicked out from the Galactic center (e.g., Dexter & O'Leary 2013), a KeplerPotential can be added with a mass of 4 × 10⁶ M_☉ (e.g., Gillessen et al. 2009) as

Finally, note that galpy contains an older version of a similar potential: MWPotential. This model was not fit to data on the Milky Way (although it is close to a good model) and is now superseded by MWPotential2014.

The integration and characterization of orbits is an essential part of galpy. Orbit integration is supported through a number of different integrators (see Section 4.2 below) and for phase-space dimensionalities from two through six. These five flavors of orbits are (1) linearOrbit for the integration of orbits in one-dimensional potentials, (2) planarOrbit and planarROrbit for orbits in non-axisymmetric and axisymmetric two-dimensional potentials, and (3) FullOrbit and RZOrbit for three-dimensional orbits. The linearOrbits are useful for investigating motions perpendicular to the mid-plane of a disk galaxy, when these motions are assumed to decouple from the motions in the plane. The axisymmetric two- and three-dimensional orbits planarROrbit and RZOrbit assume conservation of angular momentum and do not keep track of the azimuthal angle. planarOrbit and FullOrbit do keep track of the azimuthal angle and integrate the equations of motions without assuming any symmetry.

While the different flavors of orbits are all implemented as different classes deriving from a superclass OrbitTop, this structure is entirely hidden from the user and all public interfacing with instances of these classes is through a general Orbit class. All types of orbits are instantiated through this class by providing the initial conditions of the orbit and the type of orbit is determined from the dimensionality of this initial condition. Typically the initial condition is specified by an array of positions and velocities in the cylindrical Galactocentric coordinate frame (R, v_R, v_T, z, v_z, ϕ) or lower-dimensional projections of this: (x, v_x) for linearOrbit instances; (R, v_R, v_T[, ϕ]) for planarOrbit and planarROrbit objects (the latter do not specify ϕ); and (R, v_R, v_T, z, v_z[, ϕ]) for FullOrbit and RZOrbit objects. For data analysis in the Milky Way, Orbit instances can also be specified in close-to-observable coordinates: (α, δ, D, μ_{α, *}, μ_δ, v_los) or RA, Dec, distance, proper motions in (α, δ), and line-of-sight velocity (where μ_{α, *} = μ_α cos (δ)); initial conditions may also similarly be specified using Galactic coordinates and velocities can be given as (U, V, W). To use this functionality, the coordinate transformation between heliocentric and Galactocentric coordinates needs to be specified by the user by giving the Sun's distance to the Galactic center R₀ and the mid-plane z₀ and the Sun's velocity with respect to the Galactic center. The latter is split into the Sun's motion with respect to the circular velocity v_c(R₀) and v_c(R₀) itself. This way R₀ and v_c(R₀) can be used to normalize the coordinates to galpy's natural coordinates (see Section 2.2); this is done automatically. An example Orbit instantiation from observed coordinates looks like

where the solar motion is such that the Sun's rotational velocity with respect to Sgr A* is 245 km s⁻¹ (Bovy et al. 2012).

The initial conditions can be integrated in time in any of the potentials by calling the integrate method. Figure 10 demonstrates orbit integration in galpy: the two orbits displayed in the top panel of this figure have the same energy and angular momentum and are integrated in the MWPotential2014 potential; the python code that produces the top panels is given in the bottom panel as an illustration of the use of galpy. This figure is similar to Figure 3.4 in Binney & Tremaine (2008).

Zoom In Zoom Out Reset image size

After orbit integration, the orbit's characteristics and time dependence can be accessed through a variety of instance methods, which are illustrated in Figure 11. One property of Orbit instance methods is that they can produce output in physical coordinates rather than natural coordinates by specifying a distance and velocity scale through ro= and vo=, respectively. If these scales are set as keywords during the Orbit initialization, they are automatically used for any output; if not they can be specified for each individual method, which also overrides the values set at initialization. This behavior can be turned off by calling the turn_physical_off() method; this is necessary in particular for orbits initialized from observed coordinates, as these always require ro= and vo= to be set when initializing.

Zoom In Zoom Out Reset image size

Another example of how galpy's orbit routines can be used is shown in Figure 12, where the Poincaré section of the two orbits displayed in Figure 10 is computed to demonstrate that these orbits have a third integral of motion in addition to the energy and angular momentum. Calculating Poincaré sections is not currently supported by galpy, but the code in the bottom panel of Figure 12 makes clear how simple it is to extend galpy to calculate other properties of orbits.

Zoom In Zoom Out Reset image size

Orbit integration in galpy is supported using eight different integration methods, specified using the method= keyword of the integrate method. Two of these are pure python based methods. The first of these is odeint, which corresponds to scipy's ordinary-differential-equation solver of the same name. This solver uses the lsoda routine in the FORTRAN library odepack (Hindmarsh 1983). The second is leapfrog, which is a custom implementation of a leapfrog integrator, a second-order symplectic integrator (e.g., Binney & Tremaine 2008). These integrators can be used to integrate orbits in any potential implemented in the galpy Potential framework.

For faster orbit integration, higher-order solvers are implemented in C and these can be used to integrate orbits in all galpy potentials that have C implementations (see Section 3.2). Three of these integrators are Runge–Kutta solvers: the classical fourth-order Runge–Kutta method rk4_c, a fifth-order Dormand–Prince 5(4) method dopr54_c (Dormand & Prince 1980), and a sixth-order method rk6_c. galpy also contains C implementations of three symplectic integrators. The first is the same as leapfrog, but coded in C: leapfrog_c. The others are a fourth-order symplectic integrator from Forest & Ruth (1990; that corresponding to their Equation (4.9)) and the SI6A sixth-order symplectic integrator from Kinoshita et al. (1991). All of the integrators except for odeint use a constant stepsize that is set at the initial condition to produce a relative and absolute error smaller than 10⁻⁸ (combining the relative and absolute error into an overall error as abs. err. + rel. err. |(x, v)|). For this combination of position and velocity to be meaningful, use of galpy's natural coordinates (see Section 2.2) is again recommended.

The relative-energy error during the orbit integration in MWPotential2014 of the orbit displayed on the right in Figure 10 is shown in Figure 13. The energy error remains small for all of the integrators for thousands of orbits, with only odeint and rk4_c approaching relative-energy errors of 10⁻⁵ after 2000 periods. The energy error for the symplectic integrators does not grow with time as expected for phase-space-volume-conserving integrators. Because most of the phase-space of galaxies is at most a few hundred orbital times old at the present time, these energy errors are innocuous.

Zoom In Zoom Out Reset image size

galpy also has the ability to integrate phase-space volumes (Δx, Δv) rather than just positions (x, v) through the Orbit method integrate_dxdv. The phase-space volume is integrated directly, that is, by integrating the equations of motion for the phase-space volume. This requires the second derivatives of the potential to be implemented for the potential in which the volume is being integrated. This functionality is currently only supported for two-dimensional orbits, because it is used for calculating the properties of two-dimensional DFs (see Section 6.1). Extending this functionality to three-dimensional orbits would be straightforward and is planned for a future version of the code. As an example of this, Figure 14 shows how well the phase-space volume is conserved (Liouville's theorem) for the four non-symplectic integrators (symplectic integrators cannot be used for integrating (Δx, Δv) as the forces involved are not conservative). What is shown is the determinant of the Jacobian of the transformation between (Δx, Δv)(t = 0) and (Δx, Δv)(t); this determinant should be one if phase-space volume is conserved by the integrator. The fourth- and sixth-order Runge–Kutta order methods perform best for the orbit displayed in Figure 14, with much faster-growing deviations for the dopr54_c and odeint solvers. The latter in particular performs poorly.

Zoom In Zoom Out Reset image size

The various integrators contained in galpy are accessed in a uniform way (particularly the C solvers) and it is therefore easy to implement alternative methods. If these routines are cast in the same form as the currently implemented methods, adding them to the galpy framework is only a matter of adding and editing a few lines of code.

The calculation of action–angle coordinates is an integral part of galpy. They can be calculated for Orbit instances using a variety of approximations for different kinds of potentials and are used as part of some of the DFs described in Section 6. Action–angle routines are available in the galpy.actionAngle module. Each different method is implemented as a class, e.g., actionAngleSpherical (see below) that has the following methods: (1) __call__, which computes the actions only, (2) actionsFreqs for calculation of the actions and the orbital frequencies, and (3) actionsFreqsAngles, which also computes the angles. The methods are arranged in this way, because the calculation of the frequencies typically involves the computation of the actions and, similarly, computing the angles typically requires the frequencies. Grouping the outputs as is done in the three methods therefore minimizes unnecessary duplication in computations. Unless stated otherwise, all of the different types of action–angle calculations described below support all three of these methods.

The input to all three basic action–angle methods is the phase-space position (x, v). This can be specified in two different ways. The first is to pass an Orbit instance. The initial condition of this instance will be used as the phase-space position, unless a time is specified as well; in that case, the phase-space position at that time (if the orbit was integrated) will be used instead. Phase-space positions can also be specified directly in cylindrical Galactocentric coordinates. In this way, multiple phase-space positions can be passed at the same time, which is especially useful for those action–angle routines that are implemented in C, for which passing multiple points at once is more efficient.

galpy currently does not contain any methods for calculating positions and velocities for a given set of action–angle coordinates in a given potential.

Action–angle coordinates and orbital frequencies for the isochrone potential, which can be calculated analytically, and for spherical potentials, which can be calculated using a few simple numerical integrals, are implemented in galpy. The isochrone routines are contained in the actionAngleIsochrone class, which is initialized using a specific isochrone potential. This can be done either by specifying the scale parameter of the isochrone potential, in which case the potential is assumed to be in natural coordinates (i.e., the circular velocity of the isochrone potential is 1 at R = 1). Alternatively, an instance of IsochronePotential can be passed.

Action–angle coordinates for spherical potentials are contained in the actionAngleSpherical class. Instances of this class are initialized by specifying a spherical potential. Note that whether or not the potential is actually spherical is not explicitly checked by the code, so care must be taken not to use it with axisymmetric potentials. The numerical integrals can be computed either by using scipy's integrate.quad function or by using scipy's integrate.fixed_quad routine; the latter uses Gaussian quadrature and is much faster.

The performance of the isochrone and spherical action–angle methods is demonstrated in Figures 15 and 16. Figure 15 shows the radial action computed using actionAngleIsochrone and actionAngleSpherical for a similar orbit as the one displayed on the right in Figure 10. The orbit is similar in that it has the same initial condition, but it is integrated in an isochrone potential with scale parameter 0.3; this isochrone potential is similar to MWPotential2014 over the radial range covered by the orbit in Figure 10. The radial action is conserved to 1 part in 10⁻⁸ using both approximations. Figure 16 displays the error in the radial and vertical angles. What is shown is the difference between the angle at time t and the angle predicted from the linear increase of initial angle at time t = 0 with the frequency calculated as the average frequency over the orbit. Thus, this error shows how consistent angles calculated at a later time are with the initial angle and frequency. The angle error itself is ≈10⁻⁸; the frequency error of ≈2 × 10⁻⁹ leads to a linear rise in the angle error after more than a few orbital times, when frequency errors have had time to build up.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

galpy also contains multiple approximations to the actions, angles, and frequencies for axisymmetric potentials. The first of these is the adiabatic approximation (e.g., Binney 2010). In this approximation, the motion of a star near the z = 0 symmetry plane of an axisymmetric potential is approximated as being decoupled oscillators in the radial and vertical part of the effective potential. This is a good approximation for stars on nearly circular orbits that do not stray too far from the mid-plane of the potential; it is implemented as the actionAngleAdiabatic class, which is initialized using a potential instance and the specification of the γ "fudge" factor of Binney & McMillan (2011).

The second approximation for axisymmetric potentials is the Stäckel approximation of Binney (2012). This method calculates action–angle coordinates by locally approximating the potential as a Stäckel potential specified by the focal length δ of a prolate spheroidal coordinate system. As shown by Binney (2012), by specifying δ, one can calculate approximate actions, frequencies, and angles without explicitly fitting a Stäckel potential to the axisymmetric potential in question. This allows action–angle coordinates to be calculated as efficiently as for the adiabatic approximation, using just a few one-dimensional numerical integrals. This Stäckel method is implemented as the actionAngleStaeckel class, which is instantiated by giving a potential instance and a focal length δ. See Section 3.1 of Bovy & Rix (2013) for further information on the Stäckel approximation and how it is implemented in galpy.

Both the adiabatic and Stäckel approximations are implemented (partially) in python as well as in C. The actions can be calculated in python using scipy's integrate.quad function or by using scipy's integrate.fixed_quad routine; the latter is again much faster. The actions for both methods can also be computed in C, using Gaussian quadrature with 20 points. The frequencies and angles are currently not implemented for actionAngleAdiabatic, but they are available for actionAngleStaeckel. However, they are only implemented in C, and can therefore only be used with potentials that have C implementations. An interpolated potential that can be passed to C can be built using interpRZPotential (see Section 3.3) if frequencies and angles are desired for axisymmetric potentials without C implementations.

The performance of the axisymmetric approximations are demonstrated in Figures 17 and 18 for the orbit shown on the right in Figure 10 in the MWPotential2014 potential. For this orbit, actionAngleAdiabatic conserves the radial and vertical actions to a few percent, while the more precise actionAngleStaeckel display action fluctuations that are less than one percent. The consistency of the angles along the orbit with the initial angle plus the linear frequency×time increase is given in Figure 18. The angle error itself is ≈10⁻³; the inconsistency between the initial angle and the angle at time t grows with time, but only becomes of the order of one after a few hundred periods. These errors are typical for disk orbits near the Sun in the Milky Way.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As described above, actionAngleStaeckel requires the user to specify the focal length of a prolate spheroidal coordinate system. An appropriate value for this focal length can be calculated using Equation (9) in Sanders (2012), which gives the focal length of a true Stäckel potential in terms of the forces and second derivatives of the potential. If this is computed for a non-Stäckel potential, then an approximate value is found. A function estimateDeltaStaeckel included in galpy.actionAngle calculates this approximate focal length for a given position in an axisymmetric potential. If multiple positions are given, e.g., positions along an orbit, then the median of the individual focal lengths is returned. This function can be used to establish an approximate focal length for the Stäckel approximation for a given (orbit,potential) pair. In Figure 19, the approximate focal length thus computed for a wide range of orbits in MWPotential2014 is shown as a function of energy and angular momentum. This figure can be used pick δ when using actionAngleStaeckel for MWPotential2014. Very close to and very far from the center, the potential is approximately spherical, because the bulge and halo are represented with spherical models; therefore the focal length is close to zero. In regions where the disk dominates, δ ≈ 0.3 to 0.6. It is clear that δ does not vary strongly between nearby orbits, such that a single δ often suffices for a wide range of orbits near a given position.

Zoom In Zoom Out Reset image size

To speed up action calculations, both the adiabatic and the Stäckel methods can be used to tabulate the values of the radial and vertical actions on a grid in (approximate) integrals of the motion that can be interpolated for subsequent action evaluations. These grid-based methods are available as actionAngleAdiabaticGrid and actionAngleStaeckelGrid.

For actionAngleAdiabaticGrid, two separate grids for J_z and J_R are formed. The former is a grid in (R, E_z), where E_z is the vertical energy (see the caption of Figure 11 for a precise definition), because J_z only depends on R and E_z in the adiabatic approximation. The grid in R is regular over 0.01 ⩽ R ⩽ R_max. At each R_i, we then calculate E_{z, max}(R_i) ≡ E_z(z_max; R_i), evaluate J_z on a regular grid between zero and E_z(z_max; R), and divide by J_{z, max}(R_i) ≡ J_z(R_i, E_{z, max}(R_i)). We then employ third-degree bivariate spline interpolation of this normalized J_z and third-degree one-dimensional spline interpolation to interpolate the functions E_{z, max}(R) and J_{z, max}(R) as a function of R. To evaluate J_z for a given orbit (essentially [$\tilde{R},\tilde{z},\tilde{v}_z$]) using this grid, we compute $\tilde{E}_z(\tilde{R},\tilde{z},\tilde{v}_z)$, divide by the appropriate $E_{z,\mathrm{max}}(\tilde{R})$, and evaluate $J_z(\tilde{R},\tilde{E}_z/E_{z,\mathrm{max}}(\tilde{R}))/J_{z,\mathrm{max}}(\tilde{R})$ using the bivariate spline. By multiplying this by $J_{z,\mathrm{max}}(\tilde{R})$, $J_z(\tilde{R},\tilde{z},\tilde{v}_z)$ is finally obtained. The keyword parameters Rmax, zmax, nR, and nEz set the parameters R_max and z_max and the size of the grid.

J_R only depends on L_z and E_R—the radial energy, see the caption of Figure 11—in the adiabatic method. We build a regular grid in L_z over 0.01  ⩽  L_z  ⩽  v_c(R_max)R_max, using the same R_max as above. We also calculate the radius R_L of a circular orbit with angular momentum L_z at each grid point L_{z, i}. Similar to the J_z grid above, at each L_{z, i} we then compute E_{R, min, i}(L_z) ≡ E_R(v_R = 0, v_T = L_{z, i}/R_{L, i}, R = R_{L, i}) and E_{R, max, i}(L_z) ≡ E_R(v_R = 0, v_T = L_{z, i}/R_max, R = R_max) and form a regular grid in E_R at each L_{z, i} between E_{R, min, i}(L_{z, i}) and E_{R, max, i}(L_{z, i}). We subsequently evaluate J_R on this regular grid and divide by J_{R, max}(L_{z, i}) ≡ J_R(L_{z, i}, E_{R, max, i}). We again use third-degree bivariate spline interpolation to interpolate the normalized J_R and third-degree one-dimensional spline interpolation to interpolate E_{R, min}(L_z), E_{R, max}(L_z), and J_{R, max}(L_{z, i}) as functions of L_z. To evaluate J_R for a new orbit (essentially [$\tilde{R},\tilde{v}_R,\tilde{v}_T$]), we calculate $\tilde{L}_z$—as $\tilde{R}\tilde{v}_T + \gamma J_z$ using the fudge γ—and $\tilde{E}_R$ and then follow a procedure similar to that above for J_z. The size of the grid is controlled by the parameters nEr and nLz.

To interpolate the Stäckel method, we roughly follow the procedure described in Section 2.2 of Binney (2012) with some tweaks for efficiency. Because it is beyond the scope of this paper to fully describe the Stäckel method, I will employ the terminology of Binney (2012); the reader should refer to that paper for full details. We build a regular three-dimensional grid, starting with L_z over 0.01 ⩽ L_z ⩽ v_c(R_max)R_max, where R_max is again set as an input parameter Rmax. We also calculate the radius R_L of a circular orbit with angular momentum L_z at each grid point L_{z, i}. At each L_{z, i}, we then compute the energy of a circular orbit E_c(L_{z, i}) and the energy E_max(L_{z, i}) of an orbit with z = v_z = v_R = 0 and R = 25 similar to how this is done for the adiabatic grid above (this arbitrary R = 25 constant assumes that we are working in galpy's natural coordinates; see Section 2.2) and we build a logarithmic grid between these extreme energies to accurately capture orbits ranging from very-close-to-circular to radial. The Stäckel algorithm uses a reference value u₀ for the equivalent of the radial coordinate in the prolate spheroidal coordinate system. While the value of u₀ does not matter for direct evaluations of the Stäckel algorithm, it does matter for the gridded evaluation below. This is essentially because we build the grid in approximate integrals of the motion and for a good value of u₀, these integrals are better. We determine u₀ by minimizing the expression in Equation (20) of Binney (2012) on the grid in L_z and E described above and interpolate its natural logarithm using third-degree bivariate spline interpolation.

We then determine the speed w that a phase-space point at (L_z, E) and (u, v) = (u₀, π/2) has (where v is the "vertical" coordinate in the spheroidal coordinate system) and calculate the radial and vertical action on a three-dimensional grid in (L_z, E, ψ), where L_z and E are the same grid as above and ψ is a regular grid between 0 and π/2. These actions are calculated for a phase-space point (R, vR, v_T, z, v_z) = (Δsinh u₀, wcos ψ, L_z/Δ/sinh u₀, 0, wsin ψ). Similar to the adiabatic grid above, we divide each J_R and J_z two-dimensional sub-grid at L_{z, i} by its maximum and we then use three-dimensional third-degree spline interpolation to interpolate these normalized actions and one-dimensional spline interpolation to interpolate the maxima as a function of L_z. This gridding approach is slightly different from that of Binney (2012), who uses a grid in (L_z, E, E_r), where E_r is a radial energy given by Equation (18) of Binney (2012).

To evaluate the actions for a new position $(\tilde{R},\tilde{v}_R,\tilde{v}_T,\tilde{z},\tilde{v_z})$, we first calculate $\tilde{L}_z$ and $\tilde{E}$ and use the interpolation grid for u₀ to determine $\tilde{u}_0$. We then need to determine the appropriate value of $\tilde{\psi }$ to use the interpolation grid in (L_z, E, ψ). We determine two separate values for $\tilde{\psi }$ for the purpose of evaluating J_R and J_z. We compute $\tilde{\psi }_R$ by calculating $\tilde{E}_r$ and the velocity $\tilde{w}$ of the phase-space point at $(\tilde{L}_z,\tilde{E})$ and $(u,v) = (\tilde{u}_0,\pi /2)$ and determining $\tilde{\psi }_R$ from

$$\begin{equation} \cos ^2\tilde{\psi }_R = \frac{2 \tilde{E}_r}{\tilde{w}^2\cosh ^2\tilde{u}_0}. \end{equation} \tag{ 1 }$$

$\tilde{J}_R$ is then found using the interpolated J_R grid evaluated at $(\tilde{L}_z,\tilde{E},\tilde{\psi }_R)$.

For $\tilde{J}_z$ we proceed similarly, but we use a vertical energy E_z, similar to the radial energy E_r, defined as

$$\begin{equation} E_z = \frac{p^2_v}{2\Delta ^2} + \frac{L_z^2}{2\Delta ^2}\left(\frac{1}{\sin ^2 v-1}\right)-E(\sin ^2v-1)-\delta V(v), \end{equation} \tag{ 2 }$$

using the same notation as Binney (2012). We compute this $\tilde{E}_z$, find $\tilde{\psi }_z$ as

$$\begin{equation} \sin ^2\tilde{\psi }_z = \frac{2\tilde{E_z}}{\tilde{w}^2\cosh ^2\tilde{u}_0}, \end{equation} \tag{ 3 }$$

and use $(\tilde{L}_z,\tilde{E},\tilde{\psi }_z)$ to evaluate J_z.

The actions along the orbit displayed on the right in Figure 10 computed with the direct adiabatic and Stäckel methods and with their grid-based-interpolation versions are demonstrated in Figure 20.

Zoom In Zoom Out Reset image size

If the potential were an exact Stäckel potential, E_r and E_z would both be exact integrals of the motion and ψ_R would be exactly equal to ψ_z. Typical galactic potentials are not exact Stäckel potentials, so this will only hold approximately. The reason that ψ_z is a better ψ to use in the interpolation of J_z is because E_z is closely related to the constant of the motion I₃ + V(π/2) used in the direct calculation of J_z using the Stäckel method. E_r has the same relation to the constant I₃ + U(u₀), used in the direct calculation of J_R. For non-Stäckel potentials, these two constants are not interchangeable, such that using ψ_z in the interpolation gives a more accurate value of J_z. As a specific example, if we had used ψ_R instead of ψ_z to compute J_z in the bottom panel of Figure 20, the interpolated J_z would follow a straight line from the top-left to the bottom-right, because errors in J_R would directly propagate into errors in J_z, leading to a complete degeneracy.

A general method for calculating action–angle coordinates for a static potential is contained in the actionAngleIsochroneApprox class. This method works by calculating action–angle coordinates in an auxiliary isochrone potential at each point along the orbit. Actions, frequencies, and angles are then obtained by fitting—implicitly for the actions—a generating function between the (incorrect) action–angle coordinates of the isochrone potential and the (correct) action–angle coordinate of the potential of interest (see the Appendix of Bovy 2014 for full details). The actionAngleIsochroneApprox methods are built entirely on the framework provided by galpy.potential, galpy.orbit, and galpy.actionAngle.actionAngleIsochrone.

The performance of actionAngleIsochroneApprox is demonstrated in Figures 17 and 18 for the orbit in MWPotential2014 on the right in Figure 10. The actions are conserved to better than 1 part in 10⁻³. The error in the angle with respect to the initial angle plus linear increase remains small for hundreds of periods. This is essentially due to the fact that the frequency is fit simultaneously with the initial angle to the angles in the auxiliary isochrone potential over many orbital periods. The stability of the actionAngleIsochroneApprox algorithm makes it ideal for investigating the behavior of the angles over many orbital periods, such as for tidal streams (Bovy 2014).

The actionAngleIsochroneApprox method for calculating frequencies and angles is not implemented for non-axisymmetric potentials in the current version of galpy.

To illustrate the use of the action–angle routines in practice, Figure 21 demonstrates that the actions are conserved when adiabatically changing the potential. The code is given in the bottom part of the figure; it makes use of TimeInterpPotential of Figure 4 to perform the smooth change between two isochrone potentials. The energy, mean radius, and maximum height reached above the plane clearly change, but the actions are conserved.

Zoom In Zoom Out Reset image size

galpy contains both the Shu and Dehnen DFs described in Dehnen (1999) under galpy.df; both are implemented as subclasses of a general galpy.df.diskdf. These are two-dimensional, steady-state, axisymmetric DFs that use E and L_z as the sole arguments of the DF. They are both warmed up versions of a kinematically cold disk DF consisting of circular orbits only with some surface-density profile Σ(R); they differ in how this cold DF is warmed-up, i.e., generalized to include non-circular orbits.

The classical Shu DF is given by

$$\begin{equation} f_{\mathrm{Shu}}(E,L_z) = \frac{\Omega (R_L)\Sigma (R_L)}{\pi \kappa (R_L)\sigma _R^2(R_L)}\,\exp \left[\frac{E_c(L_z)-E}{\sigma _R^2(R_L)}\right], \end{equation} \tag{ 4 }$$

where Ω is the rotational frequency, κ is the epicycle frequency, Σ(· ) and σ_R(· ) are scale surface-density and radial-velocity-dispersion profiles; all of these functions are evaluated at R_L, the radius of a circular orbit with angular momentum equal to L_z. The quantity E_c(L) is the energy of this circular orbit. The Shu DF is available as galpy.df.shudf.

The Dehnen DF has the following form

$$\begin{eqnarray} f_{\mathrm{Dehnen}}(E,L_z) \ &=& \ \frac{\Omega (R_E)\Sigma (R_E)}{\pi \kappa (R_E)\sigma _R^2(R_E)}\,\nonumber\\\ms\ms && \times \exp \left[\frac{\Omega (R_E)\left(L-L_c(E)\right)}{\sigma _R^2(R_E)}\right]. \end{eqnarray} \tag{ 5 }$$

Here, R_E is the radius of a circular orbit with energy E and all of the functions appearing in this DF, which are the same as those appearing in the Shu DF, are evaluated at R_E rather than at R_L. The Dehnen DF is available as galpy.df.dehnendf.

We refer the reader to Dehnen (1999) for a detailed discussion of the advantages and disadvantages of these two DFs and of their properties. As discussed by Dehnen (1999), the surface-density Σ_DF and radial-velocity-dispersion profiles σ_{R, DF} corresponding to the Shu and Dehnen DFs differ from the scale profiles Σ(R) and σ_R(R) that appear in the definitions of these DFs. In Section 3.2 of Dehnen (1999), a procedure is given to correct the scale profiles to obtain a desired set of (Σ_out(R), σ_{R, out}(R)²)—typically these are exponential profiles characterized by a scale length and normalization at R₀—using an iterative procedure that starts with Σ(R) = Σ_out and σ_R = σ_{R, out} and applies multiplicative corrections to these based on the difference Σ_DF/Σ_out and σ_{R, DF}/σ_{R, out}. This procedure is implemented in galpy in a DFcorrection class. The user can specify the number of iterations to use and correction factors are calculated and stored automatically for re-use. Corrections are stored, because the iterative procedure is quite slow.³ A large number of correction factors stored in the galpy format are available for download online for a few profiles (Σ_out(R), σ_{R, out}(R)²). Instructions on how to download and install these are given in the online documentation.

Figure 22 demonstrates how the procedure of correcting the Σ(R) and σ_R(R) profiles in the DF leads to the desired set of (Σ_out(R), σ_{R, out}(R)²) for Σ_out(R)∝exp [ − R/h_R] and σ_{R, out}(R) = σ_R(R₀)exp [ − (R − R₀)/h_σ], with h_R = R₀/3, h_σ = R₀, and σ_R(R₀) = 0.2 (in natural coordinates). This figure shows that the correction factors are highly effective in making the differences Σ_DF/Σ_out and σ_{R, DF}/σ_{R, out} extremely small.

Zoom In Zoom Out Reset image size

The initialization of a shudf or a dehnendf instance requires one to specify the potential (through the beta= keyword that sets the logarithmic slope of the rotation curve, see below), the surface-density and radial-velocity-dispersion profiles Σ(R) and σ_R(R), and whether or not to apply corrections. When applying corrections, some keywords related to the iterative calculation of the correction factors can also be specified.

Many properties of the Shu and Dehnen DFs can be calculated by galpy. This includes the surface-density, mean velocities, velocity dispersions, as well as higher-order moments of the velocity DF. We can also calculate the Oort functions A(R), B(R), C(R), and K(R) (see, e.g., Kuijken & Tremaine 1994 for a definition of these in terms of mean velocities and their derivatives; this definition can be applied to kinematically warm populations). These are calculated by integration over the DF and radial derivatives of the DF. As an axisymmetric DF, C(R) and K(R) should be exactly zero everywhere, but they are calculated by explicit integration as a check on the DF implementation. We can further sample from the DF in a variety of ways: sampling (1) (E, L_z) or full (R, v_R, v_T, ϕ) phase-space coordinates with or without a restriction in radius (following the procedure given in Section 5 of Dehnen 1999), (2) distances along a given line of sight, (3) (v_R, v_T) velocities at a given position, or (4) full (R, v_R, v_T, ϕ) phase-space coordinates along a given line of sight. Most of the methods available for dehnendf and shudf instances are demonstrated in Figure 23.

Zoom In Zoom Out Reset image size

A limitation of the current diskdf implementation is that it only supports power-law or logarithmic potentials, i.e., power-law rotation curves (including a flat rotation curve). While it would be relatively straightforward to generalize the code to use any Potential instance, this has not been done so far.⁴ As the main purpose of the Shu and Dehnen DFs is to provide simple models for the kinematics of a disk galaxy, this is a not a serious limitation.

As an example of the functionality in the galpy.df.diskdf module, I compute the Oort functions A(R) and B(R) at the solar radius⁵ R₀ for populations with different radial-velocity dispersions and investigate how they depend on the asymmetric drift, i.e., the offset between the mean rotational velocity of a kinematically warm population and the circular velocity (Binney & Tremaine 2008). This can be achieved in galpy by running

for a fiducial model of a dehnendf with an exponential surface-density profile with scale length h_R = R₀/3, an exponential radial-velocity-dispersion profile with scale h_σ = R₀, and a flat rotation curve (beta=0 for β in v_c(R) = v_c(R₀) (R/R₀)^β). As usual, we work in natural units in which R₀ and v₀ = v_c(R₀) are both equal to one. This specific example calculates A and B for a population with a radial velocity dispersion at the Sun of σ_R(R₀) = 0.1 v₀ and by varying this parameter we can investigate A[σ_R(R₀)].

I compare A[σ_R(R₀)] and B[σ_R(R₀)] with the Oort constants, that is, A[σ_R(R₀) = 0] and B[σ_R(R₀) = 0], as a function of the asymmetric drift of the population. For the power-law rotation curve used here A[σ_R(R₀) = 0] = (1 − β)/2 and B[σ_R(R₀) = 0] = −(1 + β)/2. This comparison is presented in Figure 24 for DFs with different surface-density and velocity-dispersion profiles and with different forms (the shudf).

Zoom In Zoom Out Reset image size

It is clear from Figure 24 that both A and B depend linearly on the asymmetric drift and that A is much more affected by the kinematic temperature of the population than B. The effect of varying the DF form or parameters is much less pronounced on A than it is on B, with the only parameter making a large difference for A being h_σ. The behavior of ΔB for the Shu DF is more erratic than for the Dehnen DF, but ΔB remains small.

We can understand all of this behavior by calculating the dependence of A and B on the asymmetric drift using the radial, cylindrical Jeans equation. This equation for exponential surface-density and velocity-dispersion profiles can be written as

$$\begin{equation} v_c^2 - \bar{v}_T^2 = \sigma _R^2\,\left(\gamma ^2 + R\left[\frac{1}{h_R}+\frac{2}{h_\sigma }\right]-1\right), \end{equation} \tag{ 6 }$$

where $\gamma ^2 = \sigma _T^2/\sigma _R^2$. Taking the derivative of this, we find using (1) that γ² − 1 is typically much smaller than the other term in the parentheses and (2) that $v_c + \bar{v}_T \approx 2 v_c$

$$\begin{equation} \frac{{d}\bar{v}_T}{{d}R} = -v_a\,\left(\frac{1}{R}-\frac{2}{h_\sigma }\right)+\frac{{d}v_c}{{d}R}\,. \end{equation} \tag{ 7 }$$

Therefore, we find that

$$\begin{eqnarray} \Delta A & \approx -\frac{v_a}{h_\sigma }, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \Delta B & \approx \frac{v_a}{R}\left(1-\frac{R}{h_\sigma }\right). \end{eqnarray} \tag{ 9 }$$

These simple relations are to the author's knowledge not present in the literature. Olling & Dehnen (2003) use the Jeans equation to estimate ΔA and ΔB, but they stop short of deriving the relation above.⁶ Their value of ΔA = −0.14 v_a, where v_a is in units of  km s⁻¹ and A in units of  km s⁻¹ kpc⁻¹, computed for h_σ = 0.9R₀ agrees with Equation (8), as in our units their value is ≈0.14 × R₀ = 1.12 and Equation (8) gives 1/h_σ = 1.11. Kuijken & Tremaine (1991) derived expressions for ΔA and ΔB to leading order in $\sigma _R^2/v_c^2$; it is straightforward to demonstrate that their expressions are almost equivalent to Equations (8) and (9) (up to small terms).

Equation (8) explains why A is robust to all changes to the DF except for h_σ: in the simple approximation h_σ is the only parameter that affects ΔA. From Equation (8) we expect a linear relation between ΔA and v_a, that in the particular case of our fiducial model is one-to-one. We also expect ΔA to be negative in all cases. All of these analytical predictions are borne out in Figure 24.

Equation (9) shows why ΔB is much smaller than ΔA: for our fiducial model, ΔB vanishes in this approximation. When h_σ ≠ R₀, ΔB is much larger, as is seen in Figure 24. As for the fiducial model ΔB vanishes in the simple approximation, the actual (small) value of ΔB is determined by other factors such as the slope of the rotation curve (which affects γ²) or the scale length. The sign of ΔB can therefore be both positive and negative. In the simple approximation above, |ΔB| < |ΔA| as long as h_σ < 2R₀. While the velocity-dispersion profile in the Milky Way is likely quite flat (e.g., Bovy et al. 2012), h_σ probably does not actually exceed 2 R₀, such that we should expect ΔB to be smaller than ΔA. This implies that measurements of B in the solar neighborhood are much less affected by the effects due to non-circular motions than determinations of A (see also Olling & Dehnen 2003 who reach a similar conclusion).

A three-dimensional disk DF that is included in galpy is the quasi-isothermal DF (qDF) of Binney (2010), including the improvements made by Binney & McMillan (2011). This qDF is a steady-state, axisymmetric DF that is a function of all three of the orbital actions (J_r, L_z, J_z), where we denote the azimuthal action J_ϕ by L_z, as this action is equal to the z-component of the angular momentum for an axisymmetric potential (with rotational symmetry around the $\hat{\mathbf {z}}$ unit vector).

The quasi-isothermal DF is given by

$$\begin{equation} \mathrm{qDF}(J_r,L_z,J_z)= f_{\sigma _R}(J_r,L_z) \times \frac{\nu }{2\pi \sigma _z^2}\,\exp \left(-\frac{\nu J_z}{\sigma _z^2(R_L)}\right), \end{equation} \tag{ 10 }$$

where $f_{\sigma _R}$ is given by

$$\begin{eqnarray} f_{\sigma _R}(J_r,L_z) & =& \left.\frac{\Omega n(R_L)}{\pi \sigma _R(R_L)^2 \kappa }\right|_{R_L}\times \left[1+\tanh (L_z/L_0)\right] \nonumber\\ && \times \,\exp \left(-\frac{\kappa J_r}{\sigma _R^2(R_L)}\right). \end{eqnarray} \tag{ 11 }$$

The functions κ, Ω, and ν are the epicycle, circular, and vertical frequencies (Binney & Tremaine 2008) evaluated at R_L, where R_L is again the radius of the circular orbit with angular momentum L_z. We include the factor in Equation (11) containing the tanh  to eliminate stars on counter-rotating orbits following Binney & McMillan (2011), but there is also an initialization option to explicitly set the qDF to zero for counter-rotating orbits. The functions n, σ_R, and σ_z are free functions of R_L, which indirectly determine the radial profiles of the tracer density, the radial velocity dispersion, and the vertical dispersion. However, these are merely scale profiles, as opposed to the actual, physical profiles that can be calculated by taking the appropriate moments of the DF (similar to the difference between the scale profiles and the actual profiles for the two-dimensional DFs above). In principle these three functions can take any form, but currently galpy only supports setting each of these to an exponential an exponential

$$\begin{eqnarray} n(R_L) & \propto \exp \left(-R_L/h_R\right), \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} \sigma _R(R_L) & = \sigma _{R,0}\,\exp \left(-\left[R_L-R_0\right]/h_{\sigma _R}\right), \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} \sigma _z(R_L) & = \sigma _{z,0}\,\exp \left(-\left[R_L-R_0\right]/h_{\sigma _z}\right). \end{eqnarray} \tag{ 14 }$$

The qDF is available as galpy.df.quasiisothermaldf.

The quasiisothermaldf can be used with any Potential instance or list of such instances in galpy. The necessary action–angle calculations are performed using an instance of one of the subclasses of actionAngle that is provided at the initialization of the quasiisothermaldf instance. Any actionAngle instance can be used, but for the disk orbits that are represented by a qDF, primarily the actionAngleAdiabatic and actionAngleStaeckel as well as their gridded versions are useful.

The quasiisothermaldf implementation was used extensively in the analysis of Bovy & Rix (2013). That paper contains some details of the implementation of some of the instance methods. In Figure 25, I provide examples of some of the instance methods to illustrate the use of quasiisothermaldf instances. These include methods to evaluate the qDF itself, the DF marginalized over one or two of the velocity components, and moments of the DF such as the density, mean velocities, and velocity dispersions. A method for sampling velocities from the DF at a given position is also included to aid in generating mock data from the qDF.

Zoom In Zoom Out Reset image size

To demonstrate the use of the galpy modules with a larger example, I describe in this section the implementation in galpy of a general method for calculating the response of the stellar disk to non-axisymmetric perturbations. This method was proposed by Dehnen (2000) and consists of assuming that the DF was axisymmetric and time-independent with a known DF f₀ some time in the past and that this initially axisymmetric was then acted on by the non-axisymmetric perturbations. With these assumptions, the DF at a phase-space point (R, ϕ, v_R, v_T)—all of the discussion here is limited to two-dimensional DFs in the mid-plane of a galaxy—can be evaluated by integrating an orbit launched at this phase-space point backward in time until the time at which the DF was axisymmetric was reached; at this time, the orbit is at the point $(R^{\prime },\phi ^{\prime },v_R^{\prime },v_T^{\prime })$. By Liouville's theorem, we can then evaluate the DF f_na(R, ϕ, v_R, v_T) today as

$$\begin{equation} f_{\mathrm{na}}(R,\phi,v_R,v_T) = f_0(R^{\prime },\phi ^{\prime },v_R^{\prime },v_T^{\prime }). \end{equation} \tag{ 15 }$$

This simple procedure is implemented in the evolveddiskdf class in galpy.df. In its current implementation, it uses one of the galpy.df.diskdfs, i.e., shudf or dehnendf as f₀ and is therefore limited to power-law models of the axisymmetric part of the potential (see discussion in Section 6.1). The expression in Equation (15) allows for straightforward evaluation of the DF and its moments at a particular location, such as the mean velocities, velocity dispersions, and the vertex deviation. evolveddiskdf also contains methods to evaluate the Oort functions at a given position. In general, the Oort functions can be defined in terms of radial and azimuthal derivatives of the mean velocity (e.g., Kuijken & Tremaine 1994) and we therefore need to evaluate derivatives of f_na. These are calculated directly using the chain rule. For example,

$$\begin{eqnarray} && \frac{ \partial f_{\mathrm{na}}(R,\phi,v_R,v_T)}{\partial R} = \frac{\partial f_0(R^{\prime },\phi ^{\prime },v_R^{\prime },v_T^{\prime })}{\partial R^{\prime }}\,\frac{\partial R^{\prime }}{\partial R} \nonumber\\ && \quad +\frac{\partial f_0(R^{\prime },\phi ^{\prime },v_R^{\prime },v_T^{\prime })}{\partial v_R^{\prime }}\,\frac{\partial v_R^{\prime }}{\partial R} +\frac{\partial f_0(R^{\prime },\phi ^{\prime },v_R^{\prime },v_T^{\prime })}{\partial v_T^{\prime }}\,\frac{\partial v_T^{\prime }}{\partial R}\,, \quad\quad \end{eqnarray} \tag{ 16 }$$

where there is no ∂f₀/∂ϕ' term, because the initial DF is axisymmetric. The expression for ∂f_na/∂ϕ is similar, replacing ∂R with ∂ϕ. The partial derivatives of f₀ with respect to R, v_R, and v_T can be computed directly from the form of the initial DF. The derivatives ∂R'/∂R, etc., are computed by integrating a small phase-space volume (dx, dv) along the orbit, using the integrate_dxdv Orbit method, described in Section 4.2.

The methods of evolveddiskdf objects are similar to those of diskdf objects. One major difference is that f_na and its derivatives can be easily evaluated on a grid in (v_R, v_T) at a given (R, ϕ), allowing for quick evaluation of the moments and Oort functions.

As an example of the functionality of evolveddiskdf objects, we investigate the response of a stellar disk to an non-rotating m = 2 mode (an elliptical perturbation; Kuijken & Tremaine 1994). This is a perturbation of the form

$$\begin{equation} \phi (R,\phi) = \phi _0(R) + \frac{\epsilon (R)v^2_c(R)}{2}\,\cos 2(\phi -\phi _b), \end{equation} \tag{ 17 }$$

where ϕ₀(R) is the axisymmetric part of the potential, here modeled as a power law or a logarithm (for a flat rotation curve). With this definition, the equipotential contours are approximately elliptical with axis ratio 1 − . When assuming a radial profile for the perturbation, we use (R)∝R^{p − 2β}, where β is the power-law exponent of the rotation curve. The fiducial model below has p = β = 0.

The response of a cold disk, i.e., one with σ_R = 0, is straightforward to calculate using the method in Section 6.2.2 of Binney & Tremaine (2008). This was first done by Kuijken & Tremaine (1994) and we repeat the "cold response" here, as all of our σ_R ≠ 0 results will be with respect to the cold response. The cold response is

$$\begin{eqnarray} \bar{v}_R & = -\frac{1}{2}\left(\frac{p+2}{1-\beta }\right)\epsilon v_c \sin 2(\phi -\phi _b), \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \bar{v}_T-v_c & = -\left(\frac{1+p(1+\beta)/4}{1-\beta }\right)\epsilon v_c \cos 2(\phi -\phi _b), \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} l_v & = \frac{(p+2)[4+(1+\beta)(3p-2\beta)]}{4(1-\beta)^2(1+2\beta)}\epsilon \sin 2(\phi -\phi _b), \end{eqnarray} \tag{ 20 }$$

where l_v is the vertex deviation.

Figure 26 displays the result from a numerical calculation of the response at R₀ using the methods in evolveddiskdf. The fiducial model used in this figure is that of a flat perturbation (p = 0) with an amplitude of = 0.05 in a background potential with a flat rotation curve (β = 0) that is slowly grown in an initially steady-state disk with a Dehnen DF with h_R = R₀/3 and h_σ = R₀. From Equations (18)–(20), we can calculate that the cold response has an amplitude in $\bar{v}_R$, $\bar{v}_T-v_c$, and l_v of , , and 2 = 573, respectively. The bluest line, which has σ_R = 0.0125v_c, is very close to this cold prediction, with the only exception that the vertex deviation does not behave exactly sinusoidally. The response of warmer populations is smaller than the cold response in $\bar{v}_R$ and $\bar{v}$ and larger in l_v. The response in $\bar{v}_T$ for populations with σ_R > 0 is measured against the mean rotational velocity $\bar{v}_T^0$ of the initial DF, i.e., taking into account the asymmetric drift. Figure 26 clearly demonstrates that the response in the mean velocities remains sinusoidal, while the l_v response also keeps approximately the same shape (see below).

Zoom In Zoom Out Reset image size

To further characterize the response of a warm disk to an elliptical perturbation, we have repeated the calculation displayed in Figure 26 for variations on the fiducial model. The results from this are shown in Figure 27. This figure shows the main Fourier modes of the response normalized by the cold response as a function of σ_R. As discussed above, for $\bar{v}_R$ and $\bar{v}_T$ the only significant modes are the sin 2(ϕ − ϕ_b) and cos 2(ϕ − ϕ_b) modes, respectively. For l_v we find that we also need to include the sin 4(ϕ − ϕ_b) component, as this typically contributes about 10%.

Zoom In Zoom Out Reset image size

From linear perturbation theory, we expect the dependence of the response on σ_R to be of the form $\mathcal {F}= 1-f(\sigma _R/v_c)^2$ to lowest order in σ_R and we fit this to the fiducial model. This gives

$$\begin{eqnarray} \mathcal {F}_{v_R} & = 1- 8(\sigma _R/v_c)^2, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \mathcal {F}_{v_T} & = 1- 2.5(\sigma _R/v_c)^2, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \mathcal {F}_{l_v} & = 1+20(\sigma _R/v_c)^2, \end{eqnarray} \tag{ 23 }$$

where the index on $\mathcal {F}$ indicates the quantity that this factor applies to. For the mean velocities, reasonable variations in the DF and rotation curve only leads to corrections that are typically ≲ 10%; for $\bar{v}_R$ variations in the radial profile of the ellipticity parameterized by p are also small. The warm response in l_v is typically larger than the cold response and the magnitude of this increase depends more strongly on all of the parameters of the model. Its higher-order Fourier mode (that ∝sin 4[ϕ − ϕ_b]) also strongly depends on the model and is not linear in , as shown in the bottom two panels of the fourth column. Interpretations of the vertex deviation in terms of non-axisymmetric perturbations should therefore be treated with caution.

Finally, we also calculate the dependence of the Oort functions at R₀ on ϕ − ϕ_b and σ_R. The cold response can be computed from Equations (18)–(20) and is given explicitly in Kuijken & Tremaine (1994); we do not repeat it here. For the fiducial model, the response in all Oort functions has an amplitude of /2 and the azimuthal dependence is as cos 2(ϕ − ϕ_b) for A and B and as sin 2(ϕ − ϕ_b) for C and K. To calculate the warm response for A and B, we again subtract the axisymmetric expectation, which is not exactly equal to the axisymmetric Oort constants that would be measured for a cold population (see Section 6.2).

The result from this calculation is displayed in Figure 28. This figure demonstrates that the response for warmer populations stays sinusoidal, as expected from the sinusoidal velocity response of warm populations. The warm response in A and B is suppressed compared to the cold response and somewhat more strongly so for B. The warm response in C is stronger than the cold response. Remarkably, the response in K is independent of σ_R. Using the response factors in Equations (21)–(23) and assuming that the only radial dependence of the warm-response/cold-response is in the $\sigma _R^2$ term, we can understand these trends. The warm-response factors for the Oort functions that one computes in this way are $\mathcal {F}_A = 1\hbox{--}8.5(\sigma _R/v_c)^2$, $\mathcal {F}_B = 1\hbox{--}18.5(\sigma _R/v_c)^2$, $\mathcal {F}_C = 1+19(\sigma _R/v_c)^2$, and $\mathcal {F}_K = 1\hbox{--}13(\sigma _R/v_c)^2$. For A, B, and C these agree reasonably well with the numerical calculations displayed in Figure 28, but this estimate entirely fails to explain the lack of σ_R dependence for K. This is not due to unmodeled radial derivatives of $\mathcal {F}_R$, as from Figure 28 it can be seen that these would predict an even stronger suppression. It is likely due to higher-order terms in , because these do not need a large prefactor to become important.

Zoom In Zoom Out Reset image size

As discussed by Minchev et al. (2007) and Minchev & Quillen (2007), the dependence of the Oort constants on σ_R—foremost C and K, as these are zero for an axisymmetric disk—can be used to constrain the type of non-axisymmetry perturbing the velocity field. Minchev & Quillen (2007) demonstrate that |C| for spiral structure is smaller for warmer populations, while Minchev et al. (2007) find that |C| is larger in the case of a bar (a rotating m = 2 perturbation). Here I have explicitly demonstrated that this is also the case for a non-rotating m = 2 mode. Tentative evidence that |C| is larger for warmer populations was presented by Olling & Dehnen (2003). Gaia will allow the Oort functions to be measured with exquisite accuracy for various stellar populations and the tools in evolveddiskdf will be useful for determining the expected response and its dependence on σ_R of different non-axisymmetric agents.

As a final example of the use of evolveddiskdf, I investigate the response of the stellar disk to a weak bar. The bar model that we consider here is that of Dehnen (2000), who demonstrated that this model can explain the Hercules stream seen in the local velocity distribution. The response of the disk was further discussed in Mühlbauer & Dehnen (2003). I do not attempt as detailed a characterization as in the previous section on non-rotating m = 2 modes, but just show the response for a fiducial model.

The bar potential is modeled as a simple rotating quadrupole in DehnenBarPotential:

$$\begin{eqnarray} \Phi _b(R,\phi) &=& A_b(t) \cos [2(\phi - \Omega _b t)] \nonumber\\ && \times \left\lbrace \begin{array}{@{}ll}-(R_b/R)^3\,, & \mathrm{for}\ R \ge R_b,\\ (R/R_b)^3-2\,, & \mathrm{for}\ R \le R_b.\end{array} \right. \end{eqnarray} \tag{ 24 }$$

Here Ω_b is the pattern speed of the bar and R_b is the bar radius, chosen to be 80 percent of the bar's corotation radius. The bar amplitude is parameterized by the dimensionless parameter α, defined as

$$\begin{equation} \alpha \equiv \frac{3\,A_b}{v_c^2}\,\left(\frac{R_b}{R}\right)^3, \end{equation} \tag{ 25 }$$

and this amplitude can be increased smoothly from zero to a given bar strength. We again use a Dehnen DF with parameters h_R = R₀/3 and h_σ = R₀ to model the initial state and a logarithmic background potential with v_c(R) = v₀. We only consider the case of a disk with σ_R = 0.2 v_c. We consider the bar model with Ω_b = 1.9 Ω₀, α = 0.01, and a bar angle of 25°.

Figure 29 shows the response in $\bar{v}_R$, $\bar{v}_T$, and l_v over a large part of the disk when the bar is grown adiabatically. The radial velocity responds more strongly than the rotational velocity, with a maximum response of 0.04 v₀ for $\bar{v}_R$ versus a maximum of 0.03 v₀ for $\bar{v}_T$, both of which are obtained near the outer Lindblad resonance (R_OLR = 0.9 R₀). The vertex deviation is typically ≲ 15°, with the largest values near the outer Lindblad resonance and near R = 1.3 R₀.

Zoom In Zoom Out Reset image size

Figure 30 shows the same as Figure 29, but for a bar that is grown over only two bar periods to mimic a more realistic bar growth scenario. The response is similar to the adiabatically grown case, but the rapid bar growth induces spirality in the response, especially in $\bar{v}_R$. These figures can be compared to the results from Mühlbauer & Dehnen (2003), who used a forward-integration Monte Carlo technique to calculate the moments of the velocity distribution. The agreement between these results and the results from Mühlbauer & Dehnen (2003) is good (note that the captions of their Figures 6 and 7 appear to be switched).

Zoom In Zoom Out Reset image size

The models shown here are useful for determining the influence of the bar on observations of the mean velocity field of stars in the Milky Way. Both the elliptical–disk and bar perturbations were used by Bovy et al. (2012) to estimate the effect of non-axisymmetry on their measurement of the Milky Way's rotation curve.

All of galpy's source code is public and has been so since it was first started, in 2010 July. The source code is currently hosted in a git repository on GitHub at http://github.com/jobovy/galpy. There is no private development version. GitHub provides issue tracking that is used extensively to keep track of open issues and feature requests. Users are encouraged to report bugs through this issue tracker. As a git repository, it is easy for other users to copy the repository and create their own version of the code (and store it on GitHub); changes to the main repository are then merged through "pull requests," which is the preferred method for contributing to the codebase. galpy is released under the three-clause BSD license.

Overall, galpy consists currently of about 23,000 lines of code, with about 16,750 lines in python and 6,000 lines in C. The sizes of the major submodules are quite similar, with about 5,000 lines for each of the five main modules (actionAngle, df, orbit, potential, and util). While not incredibly large, a codebase of this size inevitably has hard-to-follow dependencies. Below I describe the automated testing framework that makes sure that changes to one part of the code do not break other parts.

galpy has an extensive set of documentation. This documentation is produced using Sphinx, which generates documentation in various format, and hosted by Read the Docs, which automatically builds the documentation upon each push to GitHub. Therefore, the documentation is always up to date. About 14,000 lines of documentation are contained in the source code, most of which are documentation of function inputs and outputs that are automatically placed in an API section of the online documentation by Sphinx. Apart from these, galpy contains another ≈5, 500 lines of reStructuredText in various tutorials. A PDF version of the online documentation currently weighs in at 283 pages.

Much effort has been put into creating a useful and meaningful test suite to ensure the reliability of the galpy code and to maintain the integrity of the code as it is extended. To do this, galpy uses the nose test environment for python. The test suite, which we discuss in more detail below, is automatically run using the Travis CI continuous integration service. This online service detects pushes to the GitHub repository and downloads the code upon each push. galpy's dependencies and galpy itself are then compiled and the test suite is run. This automatic procedure ensures that the code always compiles properly (at least on the Travis CI setup) and that all tests of the code are satisfied. The status of the code (whether it currently compiles and passes the test suite) can be checked on the galpy GitHub page.

The test suite currently consists of 500 test functions consisting of about 11,000 lines of code. These tests overall check the truth of about 20,000 assertions about the behavior of the code. All but 1,600 of these assertions are in tests of the potential and orbit-integration capabilities, because some of these tests are performed on grids of a large number of points.

Tests of galpy.potential functions are, for example, that the forces are correctly implemented as the negative derivatives of the potential and that the second derivatives of the potential are consistent with the forces (this is done using numerical differentiation). I also test that the Poisson equation is satisfied in cases where the density is implemented independently from the potential derivatives. All of these tests are automatically performed for all potentials contained in galpy. User-contributed new potentials that are registered in galpy.potential are also automatically subjected to these tests (without any need for the user to do anything). The potential test module further contains a variety of tests related to the mass, that the circular velocity and epicycle frequencies behave as expected, that the MWPotential2014 of Section 3.5 is implemented as specified in that section, etc.

Tests of the orbit-integration routines make up a significant chunk of the galpy test suite. I test energy conservation for all time-independent potentials and for all integrators and test that the Jacobi integral is conserved for the potentials with a fixed pattern speed. I further test that the energy error does not increase secularly for the symplectic integrators, by integrating orbits in a Kepler potential for 1600 periods and asserting that the slope of a linear fit to the energy errors is close to zero. I also test that Liouville's theorem is satisfied to test integrate_dxdv. Further there are various tests of the peri/apocenters, eccentricity, and maximum height of orbits in many of the potentials and of their analytical calculation, of the approximate conservation of the radial and vertical energy for orbits close to the galactic plane, of the interpolation of the orbit, and of the projection of the orbit in various coordinate and unit systems.

The methods of actionAngle instances are tested by checking that actions are conserved along orbits in static potentials, that the radial and vertical actions are very small for close-to-circular orbits, that the frequencies of close-to-circular orbits are approximately the epicycle frequencies (radial, rotational, and vertical), that the angles increase linearly in time, and that the slope of this increase is given by the frequencies. We further test that the action–angle coordinates computed by all methods agree for an isochrone potential, by checking that the actions, frequencies, and angles are the same when calculated by the method in question and actionAngleIsochrone. Thus, we know that all different methods are implemented consistently.

The methods in galpy.df.diskdf are tested by calculating moments of the DF and the Oort functions for very cold populations and comparing them to analytical predictions. This is done at a few different positions and for rotation curves that are flat, falling, and rising. We also check that the surface-density routines work as expected, both for cold and warm populations, and that the sampling methods return Monte Carlo samples that are consistent with the moments of the DF. The methods in galpy.df.quasiisothermaldf are also tested against analytical predictions for an axisymmetric DF and, for orbits close to the plane, against the calculations in galpy.df.diskdf (which should be similar). We also test that the sampling routine is consistent with the moments of the DF and that the methods that compute the DF marginalized over one or two velocity components are consistent with a simple Riemann sum of the DF. The methods in galpy.df.evolveddiskdf are primarily tested for potentials that are very close to axisymmetric: for such potentials the moments of the DF should agree with those of the initial diskdf. We also check the computations for an elliptical-disk perturbation applied to a very cold population against the analytical predictions (see Section 7 above).

The utility functions in galpy.util are also well tested. All of the coordinate transformations contained in bovy_coords are tested for various points (such as the North Galactic Pole for radec_to_lb). The conversion factors in bovy_conversion are tested by checking that they scale correctly with the distance and velocity scales. The plotting routines in bovy_plot are not tested, because it is hard to meaningfully test plotting routines. Plotting routines in the rest of the galpy module are tested by making sure that they run without raising exceptions.

Finally, tests are run that all of the code examples given in this paper produce the output given here, such that these examples will remain valid in future updates to the code (or at least that any changes will be known).

The amount of the galpy code that is covered by all of these tests is tracked using coverage tools (coverage.py for the python code and gcov for the C code) and automatically reported online after Travis CI runs the test suite. The coverage of the code is currently 99.6%, that is, 14,166 of 14,220 relevant lines are touched by the test code (relevant lines here indicates individual statements and combines multi-line statements into one relevant line; this is one of the main reasons for the discrepancy between this number of lines and the number of 23,000 lines of code quoted above for the whole galpy module). All python files are covered at more than 99%. The only lines that are not covered in the whole codebase are related to special cases that are very unlikely to come up in practice and paths through the code that are currently impossible to access through the user interface (but are left for completeness).