Why?

>>> import astropy.units as u
>>> import numpy as np
>>> import gala.integrate as gi
>>> import gala.dynamics as gd
>>> import gala.potential as gp
>>> from gala.units import galactic

Commonly used Galactic gravitational potentials

Create potential objects that know how to compute the gradient, energy, etc. at a given position:

>>> pot = gp.IsochronePotential(m=1E10*u.Msun, b=15.*u.kpc, units=galactic)
>>> pot.value([8.,6.,7.]*u.kpc) 
<Quantity [-0.00131039] kpc2 / Myr2>
>>> pot.gradient([8.,6.,7.]*u.kpc) 
<Quantity [[  1.57857635e-05],
           [  1.18393226e-05],
           [  1.38125430e-05]] kpc / Myr2>

Extensible and easy to define new potentials

New potentials can be easily defined by subclassing the base potential class, PotentialBase. For faster orbit integration and computation, you can also define potentials with functions that evaluate its derived quantities in C by subclassing CPotentialBase. For fast creation of potentials for quick testing, you can also create a potential class directly from an equation that expresses the potential:

>>> SHOPotential = gp.from_equation("1/2*k*x**2", vars="x", pars="k",
...                                 name='HarmonicOscillator')

(note: this requires sympy).

Classes created this way can then be instantiated and used like any other PotentialBase subclass:

>>> pot = SHOPotential(k=1.)
>>> pot.value([1.])
<Quantity [ 0.5]>

Extremely fast orbit integration

Most of the guts of the built-in potential classes are implemented in C, enabling extremely fast orbit integration for single or composite potentials:

>>> pot = gp.IsochronePotential(m=1E10*u.Msun, b=15.*u.kpc, units=galactic)
>>> w0 = gd.CartesianPhaseSpacePosition(pos=[7.,0,0]*u.kpc,
...                                     vel=[0.,50.,0]*u.km/u.s)
>>> import timeit
>>> timeit.timeit(lambda: pot.integrate_orbit(w0, dt=0.5, n_steps=10000), number=100) / 100. 
0.0028513244865462184

For a composite potential:

>>> bulge = gp.IsochronePotential(m=2E10*u.Msun, b=0.5*u.kpc, units=galactic)
>>> disk = gp.MiyamotoNagaiPotential(m=6E10*u.Msun, a=3*u.kpc, b=0.26*u.kpc, units=galactic)
>>> pot = gp.CCompositePotential(bulge=bulge, disk=disk)
>>> timeit.timeit(lambda: pot.integrate_orbit(w0, dt=0.5, n_steps=10000), number=100) / 100. 
0.0031369362445548177

Precise integrators

The default orbit integration routine uses LeapfrogIntegrator, but the high-order Dormand-Prince 853 integration scheme is also implemented as DOPRI853Integrator:

>>> orbit = pot.integrate_orbit(w0, dt=0.5, n_steps=10000,
...                             Integrator=gi.DOPRI853Integrator)

Easy visualization

Numerically integrated orbits can be easily visualized using the plot() method:

>>> orbit.plot()

Astropy units support

All functions and classes have Astropy unit support built in: they accept and return Quantity objects wherever possible. In addition, this package uses an experimental new UnitSystem class for storing systems of units and default representations.

Astropy coordinates support

Gala also contains functionality for transforming velocities between certain Astropy coordinate frames. See Coordinate Systems (gala.coordinates) for more information.