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()
sim.add("Sun")
sim.add("Jupiter")
sim.add("Saturn")

```
```

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():
print(orbit)

```
```

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]:
```sim.add("Churyumov-Gerasimenko")

```
```

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

```
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):
sim.integrate(time)
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)
ax.set_xlim([-6,6])
ax.set_ylim([-6,6])
plt.plot(x[0], y[0]);
plt.plot(x[1], y[1]);

```
```

```
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))

```
```

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

```
Out[10]:
```