This tutorial teaches you how to use the IAS15 integator (Rein and Spiegel, 2015) to simulate the orbit of 67P/Churyumov–Gerasimenko. We will download the data from NASA Horizons and visualize the orbit using matplotlib.

This tutorial assumes that you have already installed REBOUND.

NASA Horizons

If you're interested in Solar System dynamics, you have probably heard of NASA Horizons. It's a large database of Solar System objects, their orbits and physical properties. It includes planets, moons, satellites, asteroids, comets and spacecrafts. With REBOUND, you can easily import data from NASA Horizons. As an example, let's pull in the present day positions of Jupiter, Saturn and the Sun:

In [1]:
import rebound
sim = rebound.Simulation()

Searching NASA Horizons for 'Sun'... Found: Sun (10).
Searching NASA Horizons for 'Jupiter'... Found: Jupiter Barycenter (5).
Searching NASA Horizons for 'Saturn'... Found: Saturn Barycenter (6).

Now all the data is in REBOUND! Let's have a look at the orbits of the two planets.

In [2]:
for orbit in sim.calculate_orbits():

<rebound.Orbit instance, a=5.20203714989 e=0.0489258266926 inc=0.022754492909 Omega=1.75429866032 omega=-1.50402911875 f=2.42176555887>
<rebound.Orbit instance, a=9.54679703442 e=0.0539199147952 inc=0.0434070818366 Omega=1.98278944147 omega=-0.364701706059 f=-3.63607179799>

Although there are three bodies, the calculate_orbits() function only returns two objects as the orbit for the Sun would be a little boring. The function returns the orbits in Jacobi coordinates. Since we didn't specify a falue for $G$, REBOUND assumes that $G=1$. The unit of length is one astronomical unit, the unit of time is one year/$2\pi$.

Let's add something more interesting to our simulation: the comet 67P/Churyumov-Gerasimenko.

In [3]:

Searching NASA Horizons for 'Churyumov-Gerasimenko'... Found: 67P/Churyumov-Gerasimenko.
Warning: Mass cannot be retrieved from NASA HORIZONS. Set to 0.

When searching for a body by name, REBOUND takes the first dataset that Horizons offers. In this case, it's a set of parameters from 1962. You probably want to go to the Horizons website and check that the values you are using are up-to-date and appropriate for what you want to do. You can also use more complicted Horizons queries, for example, to get the most recent apparition solution for the comet, use

sim.add("NAME=Churyumov-Gerasimenko; CAP")

You can also use the IAU asteroid number for numbered asteroids, or the database record numbers from Horizons for objects not yet numbered by the IAU (but note that database record numbers can change as the database gets rearranged with new discoveries, see http://ssd.jpl.nasa.gov/?horizons_doc#sb for details). In our case the current database record number is 900647, so you could use sim.add("900647") to get the newest set of orbital parameters for Churyumov-Gerasimenko.

NASA Horizons doesn't have masses for all bodies. If REBOUND doesn't find a mass, you get a warning message (see above). In our case, we don't need the mass of the comet (it's really smal). However, it you want, you can add it manually using the syntax sim.add("Churyumov-Gerasimenko", m=5.03e-18).

Before we integrate the orbits, let's plot the instantaneous orbits using the built-in REBOUND function OrbitPlot.

In [4]:
%matplotlib inline
fig = rebound.OrbitPlot(sim, trails=True, unitlabel="[AU]")

Integration with IAS15

We will integrate backwards in time for 70 years. Because we don't know what will happen yet (hint: a close encounter) we will use the IAS15 integrator. It is fast, accurate and has adaptive timesteps to capture any potential close encounters.

To integrate backwards, we could set a negative timestep or multiply all velocities with $-1$. We'll choose the first option:

In [5]:
sim.dt = -0.01

While we're integrating, let's store the positions of Jupiter and the comet at 1000 times during the interval. We'll need to prepare a few variables to do that:

In [6]:
import numpy as np
Noutputs = 1000
year = 2.*np.pi # One year in units where G=1
times = np.linspace(0.,-70.*year, Noutputs)
x = np.zeros((2,Noutputs))
y = np.zeros((2,Noutputs))

Now we're ready to start the integration:

In [7]:
sim.integrator = "ias15" # IAS15 is the default integrator, so we actually don't need this line
sim.move_to_com()        # We always move to the center of momentum frame before an integration
ps = sim.particles       # ps is now an array of pointers and will change as the simulation runs

for i,time in enumerate(times):
    x[0][i] = ps[1].x   # This stores the data which allows us to plot it later
    y[0][i] = ps[1].y
    x[1][i] = ps[3].x
    y[1][i] = ps[3].y

Visualization with matplotlib

Let's plot the orbits of Jupiter (blue) and the comet (green) to get an idea of what was going on during our integration.

In [8]:
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,5))
ax = plt.subplot(111)
plt.plot(x[0], y[0]);
plt.plot(x[1], y[1]);

As you can see in the above image, the comet 67P had a rather strong encounter with Jupiter a few years ago. Of course, if you wanted to do a realistic simulation of that encounter, you'd need to include all the other planets and maybe even some non-gravitational effects for the comet. However, let's stick with our simplistic model and try to find out when exactly the two bodies had a close encouter. We already stored the data, so we can just plot their distance as a function of time.

In [9]:
fig = plt.figure(figsize=(12,5))
ax = plt.subplot(111)
ax.set_xlabel("time [yrs]")
ax.set_ylabel("distance [AU]")
distance = np.sqrt(np.square(x[0]-x[1])+np.square(y[0]-y[1]))
plt.plot(times/year, distance);
closeencountertime = times[np.argmin(distance)]/year
print("Minimum distance (%f AU) occured at time: %f years." % (np.min(distance),closeencountertime))

Minimum distance (0.023324 AU) occured at time: -56.546547 years.

We can see that the minimum distance occured approximately 56 years ago (as of writing this tutorial). Let's see what date that was using some python magic and the datetime module:

In [10]:
import datetime
encounterdate = datetime.datetime.today() + datetime.timedelta(days=365.25*closeencountertime)
encounterdate.strftime("%Y-%m-%d %H:%M")

'1959-02-07 19:05'

According our calculation, the close encounter happened in February 1959. If you check wikipedia it turns out we are (at least approximately) right!