Getting started with ASSIST: Integrating Asteroid (84100) Farnocchia¶
Here we integrate Asteroid (84100) Farnocchia, a main belt asteroid. This is a good notebook to get started with ASSIST. ASSIST provides a set of python-wrapped C functions that implement an ephemeris quality-integrator. ASSIST can integrate massless test particles in the field of the Sun, planets, Moon, and 16 massive asteroids. It also includes the J2, J3, and J4 gravitational harmonics of the Earth, the J2 gravitational harmonic of the Sun, and the solar general relativistic correction terms. ASSIST makes use of the REBOUND integrator package. The underlying numerical integrator is IAS15 (Rein & Liu 2015), a 15th order predictor-corrector integrator with an adaptive step-size.
The positions of the massive bodies come from two binary files, both provided by NASA JPL. The first is for the Sun, planets, and Moon, with the latest DE441 ephemeris. The other is for the asteroids, corresponding to DE441. To use assist, you need to download these files and store them on your computer. They can be found at these URLs:
- https://ssd.jpl.nasa.gov/ftp/eph/planets/Linux/de441/linux_m13000p17000.441
- https://ssd.jpl.nasa.gov/ftp/eph/small_bodies/asteroids_de441/sb441-n16.bsp
Note that these are large files, more than 3 GB in total. Instead of the DE441 file, you can also download the DE440 file which is only 100 MB and covers a shorter timespan:
Once you have installed ASSIST and downloaded these files, you're ready to go ahead with this notebook.
First, let's import REBOUND, ASSIST as well as numpy and matplotlib.
import numpy as np
import matplotlib.pyplot as plt
import rebound
import assist
#au_km = 149597870.700
Next, we create an ephemeris structure. ASSIST will use this structure to access the JPL Ephemeris files during the integrations. Depending on where you have stored the files, you might need to adjust the paths.
ephem = assist.Ephem("../data/linux_p1550p2650.440", "../data/sb441-n16.bsp")
ASSIST measures time relative to some reference time to assure high accuracy during the integrations. The default reference time is (in Julian Days):
print(ephem.jd_ref)
2451545.0
We can now create a REBOUND simulation and "attach" ASSIST. This will point the forces calculation routine of the REBOUND simulation to ASSIST and setup a few data structures needed for the integration.
sim = rebound.Simulation()
ex = assist.Extras(sim, ephem)
We are now ready to add particles to our simulation. In contrast to a normal REBOUND simulation, the coordinate frame and units are not flexible when using ASSIST. The coordinate frame is the equatorial International Celestial Reference Frame (ICRF). This is also the native coordinate system for the JPL binary files. Note that this is equatorial rather than ecliptic which means that orbits with zero inclination have a finite z value. In addition, the native coordinates are barycentric, rather than heliocentric.
For units we use solar masses, astronomical units, and days. The time coordinate is Barycentric Dynamical Time (TDB) in Julian days.
sim.add(x = 9.572786624350362E-01, y=-2.101947659454845E+00, z=-7.418444809938682E-01,
vx = 9.715385994801106E-03, vy=6.153629531316274E-03, vz=1.549521077070705E-03)
The above coordinates correspond to the positions and velocities of asteroid (84100) Farnocchia at time 2458849.5 (2020-Jan-01). We need to set the time in our simulation accordingly (remember, time is always relative to jd_ref
):
sim.t = 2458849.5 - ephem.jd_ref
We can now integrate the trajectory 10 000 days forward in time with reb_integate
:
sim.integrate(sim.t + 10_000)
We can now print out the final position and velocity of the asteroid.
print(sim.particles[0])
<rebound.particle.Particle object at 0x11302aa40, m=0.0 x=1.810998201450716 y=-1.2540881214261608 z=-0.5062323305041747 vx=0.006256088191987634 vy=0.010224229837940246 vz=0.0030657114679789417>