11. N-body simulation utilities

In addition to the main features designed for creating and working with symbolic expressions, celmech provides a number of utility functions for working with REBOUND simulations.

The get_simarchive_integration_results function provides a convenient means of reading the data stored in a rebound SimulationArchive. The method takes a simulation archive file as input and returns a dictionary containing the times of each snapshot along with the orbital elements of each planet in the simulation at those snapshots. We illustrate this with a short example below.

First, we run a short simulation of a pair of planets, saving results to a simulation archive:

import rebound as rb
sim = rb.Simulation()
sim.add(m=1)
sim.add(m=1e-4,P=1)
sim.add(m=1e-4,P=2,e=0.1)
sim.automateSimulationArchive("save.sa",interval = 10.)
sim.integrate(1e4)

Now, we’ll use get_simarchive_integration_results to read the results of our simulation.

from celmech.nbody_simulation_utilities import get_simarchive_integration_results
results = get_simarchive_integration_results("save.sa")

To see the contents of results, we can print the dictionary keys:

print(results.keys())

which should produce

dict_keys(['time', 'P', 'e', 'l', 'inc', 'pomega', 'omega', 'Omega', 'a', 'Energy'])

Finally, we’ll use the contents of our results dictionary to plot the planets’ eccentricities versus time:

import matplotlib.pyplot as plt
plt.plot(results['time'],results['e'][0],label="$e_1$")
plt.plot(results['time'],results['e'][1],label="$e_2$")
plt.xlabel("Time")
plt.ylabel("Eccentricity")
plt.legend()

and we get the following plot:

Other convenient simulations utilities include:

1. set_time_step sets the time step of a simulation to a user-specified fraction of the shortest perihelion passage timescale.

2. calculate_mutual_inclinations computes the mutual inclinations between pairs of planets in a simulation.

3. align_simulation aligns a rebound simulation so that the total angular momentum is oriented along the \(z\) axis.

11.1. API

celmech.nbody_simulation_utilities.add_canonical_heliocentric_elements_particle(m, elements, sim)

Add a new particle to a rebound simulation by specifying its mass and canonical heliocentric orbital elements.

Parameters

• m (float) – Physical mass of particle to add.

• elements (dict) – Dictionary of orbital elements for particle. Dictionary must contain valid set of orbital elements accepted by REBOUND.

• sim (rebound.Simulation) – Simulation to add particle to.

celmech.nbody_simulation_utilities.align_simulation(sim)

Change particle positions and velocities of a rebound simulations so that the z-axis corresponds with the direction of the angular momentum.

Parameters

sim (rebound.Simulation) –

celmech.nbody_simulation_utilities.calculate_mutual_inclinations(s)

Calculate the mutual inclination between pairs of bodies in a simulation archive

Parameters

s (rebound.Simulation or rebound.SimulationArchive) –

Returns

imutuals – Dictionary of mutual inclinations. Keys are pairs of particle numbers, (i,j), where i and j are integers and values are arrays of mutual inclination values for that key pair.

Return type

dictionary

celmech.nbody_simulation_utilities.get_canonical_heliocentric_orbits(sim)

Compute orbital elements in canonical heliocentric coordinates (e.g., Laskar & Robutel 1995), in the center-of-mass frame.

11. Arguments:

simrb.Simulation

simulation for which to compute orbits

returns

Orbits of particles in canonical heliocentric coordinates.

rtype

list of rebound.Orbits

celmech.nbody_simulation_utilities.get_simarchive_integration_results(sa, coordinates='jacobi')

Read a simulation archive and store orbital elements as arrays in a dictionary.

Parameters

• sa (rebound.SimulationArchive or str) – The simulation archive to read or the file name of the simulation archive file. Can also be a reboundx simulation archive.

• coordinates (str) – The coordinate system to use for calculating orbital elements. This can be: | - ‘jacobi’ : Use Jacobi coordinates (including Jacobi masses) | - ‘barycentric’ : Use barycentric coordinates. | - ‘heliocentric’ : Use canonical heliocentric elements. | The canonical cooridantes are heliocentric distance vectors. | The conjugate momenta are barycentric momenta.

Returns

sim_results – Dictionary containing time and orbital elements at each snapshot of the simulation archive.

Return type

dict

celmech.nbody_simulation_utilities.set_time_step(sim, dtfactor)

Set the time step of a simulation to a fraction,``dtfactor``, of the minimum perihelion passage timescale, defined as

\[\tau_{p} = P\sqrt{\frac{(1-e)^3}{1+e}}~.\]

Generally, the time step of WHFAST should be set to \(\lesssim\tau_p/20\) in order to obtain reliable results (see Wisdom 2015)