Integrate an orbit and compute uncertainties in Milky Way potential model

gala provides a simple mass model for the Milky Way based on recent measurements of the enclosed mass compiled from the literature. See the Defining a Milky Way potential model documentation for more information about how this model was defined.

In this example, we’ll use the position and velocity and uncertainties of the Milky Way satellite galaxy “Draco” to integrate orbits in a Milky Way mass model starting from samples from the error distribution over initial conditions defined by its observed kinematics. We’ll then compute distributions of orbital properties like orbital period, pericenter, and eccentricity.

In [1]:
# Some imports we'll need later:

# Third-party
import astropy.units as u
import astropy.coordinates as coord
from import ascii
import matplotlib.pyplot as plt
import numpy as np

# Gala
from gala.mpl_style import mpl_style
import gala.dynamics as gd
import gala.integrate as gi
import gala.potential as gp
from gala.units import galactic
%matplotlib inline

For the Milky Way model, we’ll use the built-in potential class in gala (see above for definition):

In [2]:
potential = gp.MilkyWayPotential()

For the sky position and distance of Draco, we’ll use measurements from Bonanos et al. 2004. For proper motion components, we’ll use the recent HSTPROMO measurements (Sohn et al. 2017) and the line-of-sight velocity from Walker et al. 2007.

In [3]:
icrs = coord.ICRS(ra=coord.Angle('17h 20m 12.4s'),
                  dec=coord.Angle('+57° 54′ 55″'),

icrs_err = coord.ICRS(ra=0*u.deg, dec=0*u.deg, distance=6*u.kpc,

Let’s start by transforming the measured values to a Galactocentric reference frame so we can integrate an orbit in our Milky Way model. We’ll do this using the velocity transformation support in `astropy.coordinates <>`__ (new to Astropy v2.0). We first have to define the position and motion of the sun relative to the Galactocentric frame, and create an `astropy.coordinates.Galactocentric <>`__ object with these parameters:

In [4]:
v_sun = coord.CartesianDifferential([11.1, 250, 7.25]*
gc_frame = coord.Galactocentric(galcen_distance=8.3*u.kpc,

To transform the mean observed kinematics to this frame, we simply do:

In [5]:
gc = icrs.transform_to(gc_frame)

That’s it! Now we have to turn the resulting Galactocentric object into orbital initial conditions, and integrate the orbit in our Milky Way model. We’ll use a timestep of 0.5 Myr and integrate the orbit backwards for 10000 steps (5 Gyr):

In [6]:
w0 = gd.PhaseSpacePosition(
orbit = potential.integrate_orbit(w0, dt=-0.5*u.Myr, n_steps=10000)

Let’s visualize the orbit:

In [7]:
fig = orbit.plot()

With the orbit object, we can easily compute quantities like the pericenter, apocenter, or eccentricity of the orbit:

In [8]:
orbit.pericenter(), orbit.apocenter(), orbit.eccentricity()
(<Quantity 46.817235533248734 kpc>,
 <Quantity 104.79117839058257 kpc>,
 <Quantity 0.38239264798627987>)

We can also use these functions to get the time of each pericenter or apocenter - let’s plot the time of pericenter, and time of apocenter over the time series of the Galactocentric radius of the orbit:

In [9]:
plt.plot(orbit.t, orbit.spherical.distance, marker='None')

per, per_times = orbit.pericenter(return_times=True)
apo, apo_times = orbit.apocenter(return_times=True)

for t in per_times:
    plt.axvline(t.value, color='#67a9cf')

for t in apo_times:
    plt.axvline(t.value, color='#ef8a62')

plt.xlabel('$t$ [{0}]'.format(orbit.t.unit.to_string('latex')))
plt.ylabel('$r$ [{0}]'.format(orbit.x.unit.to_string('latex')))
<matplotlib.text.Text at 0x7f005ff557f0>

Now we’ll sample from the error distribution over the distance, proper motions, and radial velocity, compute orbits, and plot distributions of mean pericenter and apocenter:

In [10]:
n_samples = 256
dist = np.random.normal(icrs.distance.value, icrs_err.distance.value,
                        n_samples) * icrs.distance.unit

pm_ra_cosdec = np.random.normal(icrs.pm_ra_cosdec.value,
                                n_samples) * icrs.pm_ra_cosdec.unit

pm_dec = np.random.normal(icrs.pm_dec.value,
                          n_samples) * icrs.pm_dec.unit

rv = np.random.normal(icrs.radial_velocity.value, icrs_err.radial_velocity.value,
                      n_samples) * icrs.radial_velocity.unit

ra = np.full(n_samples, *
dec = np.full(n_samples, *
In [11]:
icrs_samples = coord.ICRS(ra=ra, dec=dec, distance=dist,
                          pm_dec=pm_dec, radial_velocity=rv)
In [12]:
In [13]:
gc_samples = icrs_samples.transform_to(gc_frame)
In [14]:
w0_samples = gd.PhaseSpacePosition(
orbit_samples = potential.integrate_orbit(w0_samples, dt=-0.5*u.Myr, n_steps=10000)
In [15]:
(10001, 256)
In [16]:
pers = []
apos = []
eccs = []
for n in range(n_samples):
    orbit = orbit_samples[:,n]

pers = u.Quantity(pers)
apos = u.Quantity(apos)
eccs = u.Quantity(eccs)
In [17]:
fig, axes = plt.subplots(1, 3, figsize=(12, 4))

axes[0].hist(pers, bins='auto')
axes[0].set_xlabel('pericenter [kpc]')

axes[1].hist(apos, bins='auto')
axes[1].set_xlabel('apocenter [kpc]')

axes[2].hist(eccs, bins='auto')