1. Installation

1.1. Requirements

Installation of celmech will require a few additional Python packages.

1.2. GitHub

In order to install celmech you will first need to clone the GitHub repository located here. This can be accomplished from the terminal by first navigating to the desired target directory and then issuing the command:

git clone https://github.com/shadden/celmech.git

After you have cloned the git repository, navigate into the top celmech directory and issue the terminal command:

python setup.py install

in order to install the Python package.

1.3. Pip

celmech is available through PyPI and can be installed with the terminal command:

pip install celmech

2. A First Example

Now that celmech is installed, we’ll run through a short example of how to use it. In our example we’ll use it in conjunction with the rebound N-body integrator to build and integrate a simple Hamiltonian model for the dynamics of a pair of planets in a mean motion resonance. The example presented below is also available as a Jupyter notebook on GitHub.

2.1. Setup

We’ll start by importing the requisite packages.

import numpy as np
import rebound as rb
from celmech import Poincare, PoincareHamiltonian
from sympy import init_printing
init_printing() # This will typeset symbolic expressions in LaTeX

Now we’ll initialize a rebound simulation containing a pair of Earth-mass planets orbiting a central solar-mass star with a 3:2 period ratio commensurability.

sim = rb.Simulation()
sim.add(m=1,hash='star')
sim.add(m=3e-6,P = 1, e = 0.03,l=0)
sim.add(m=3e-6,P = 3/2, e = 0.03,l=np.pi/5,pomega = np.pi)
sim.move_to_com()
rb.OrbitPlot(sim,periastron=True)

After running the above code, we should see a represntation of our planetary system:

_images/quickstart_orbit_plot.png

We’d like a construct a Hamiltonian that can capture the dynamical evolution of the system. To start, we’ll initialize a Poincare object in order to represent our system and a PoincareHamiltonian object to model its dynamical evolution. These objects can be initialized directly from the rebound simulation that we created above:

pvars = Poincare.from_Simulation(sim)
pham = PoincareHamiltonian(pvars)

Now that we’ve now generated a PoincareHamiltonian object, we’ll examine the symbolic expression for the Hamiltonian governing our system:

pham.H

which should display:

\[- \frac{G^{2} M_{2}^{2} m_{2}^{3}}{2 \Lambda_{2}^{2}} - \frac{G^{2} M_{1}^{2} m_{1}^{3}}{2 \Lambda_{1}^{2}}\]

This expression is the just Hamiltonian of two non-interacting Keplerian orbits expressed in canonical variables used by celmech. The canonical momenta for the \(i\)-th planet are defined 1 in terms of the planet’s standard orbital elements \((a_i,e_i,I_i,\lambda_i,\varpi_i,\Omega_i)\) and mass parameters \(\mu_i\sim m_i\) and \(M_i \sim M_*\):

\[\begin{split}\begin{align*} \Lambda_i &= \mu_i \sqrt{G M_i a_i}\\ \kappa_i &= \sqrt{2\Lambda_i(1-\sqrt{1-e_i^2})}\cos\varpi_i\\ \sigma_i &= \sqrt{2\Lambda_i\sqrt{1-e_i^2}(1-\cos I_i)}\cos\Omega_i \end{align*}\end{split}\]

and their conjugate coordinates are:

\[\begin{split}\begin{align*} \lambda_i & \\ \eta_i &= -\sqrt{2\Lambda_i(1-\sqrt{1-e_i^2})}\sin\varpi_i\\ \rho_i &= -\sqrt{2\Lambda_i\sqrt{1-e_i^2}(1-\cos I_i)}\sin\Omega_i \end{align*}\end{split}\]

When a PoincareHamiltonian is first initialized, it will only contain the ‘Keplerian’ terms of the Hamiltonian and will not contain any terms representing gravitaional interactions between the planets. This will result in quite boring dynamical evolution: the planets’ mean longitudes, \(\lambda_i\), will simply increase linearly with time at a rate of \(n_i = \frac{G^{2} M_{2}^{2} m_{i}^{3}}{\Lambda_{i}^{3}}\), while all other orbital elements remain constant.

In order explore more interesting dynamics, we need to add terms to Hamiltonian that capture pieces of the gravitational interactions between planets. Since our planet pair is near a 3:2 MMR, terms associated with this resonance are a natural choice to explore. For a pair of co-planar planets, these terms will all involve linear combinations of the two resonant angles

\[\begin{split}\theta_1 = 3\lambda_2-2\lambda_1 - \varpi_1 \\ \theta_2 = 3\lambda_2-2\lambda_1 - \varpi_2\end{split}\]

In fact, at lowest order in the planets’ eccentricities, there are just two such terms, \(\propto e_1\cos\theta_1\) and \(\propto e_2\cos\theta_2\). The method add_MMR_terms() provides a convenient method for adding these terms to our Hamiltonian:

pham.add_MMR_terms(3,1,max_order=1)
pham.H

which should now display

\[- \frac{C^{0,0,0,0;(1,2)}_{0,0,0,0,0,0} G^{2} M_{2}^{2} m_{1}}{\Lambda_{2}^{2} M_{1}} m_{2}^{3} - \frac{C^{0,0,0,0;(1,2)}_{3,-2,-1,0,0,0} G^{2} M_{2}^{2} m_{1}}{\Lambda_{2}^{2} M_{1}} m_{2}^{3} \left(\frac{\eta_{1}}{\sqrt{\Lambda_{1}}} \sin{\left (2 \lambda_{1} - 3 \lambda_{2} \right )} + \frac{\kappa_{1}}{\sqrt{\Lambda_{1}}} \cos{\left (2 \lambda_{1} - 3 \lambda_{2} \right )}\right) - \frac{C^{0,0,0,0;(1,2)}_{3,-2,0,-1,0,0} G^{2} M_{2}^{2} m_{1}}{\Lambda_{2}^{2} M_{1}} m_{2}^{3} \left(\frac{\eta_{2}}{\sqrt{\Lambda_{2}}} \sin{\left (2 \lambda_{1} - 3 \lambda_{2} \right )} + \frac{\kappa_{2}}{\sqrt{\Lambda_{2}}} \cos{\left (2 \lambda_{1} - 3 \lambda_{2} \right )}\right) - \frac{G^{2} M_{2}^{2} m_{2}^{3}}{2 \Lambda_{2}^{2}} - \frac{G^{2} M_{1}^{2} m_{1}^{3}}{2 \Lambda_{1}^{2}}\]

This somewhat cumbersome expression is just equivalent to

\[- \frac{GM_*m_1}{2 a_1} - \frac{GM_*m_2}{2 a_2} - \frac{Gm_1m_2}{a_2}\left(C^{0,0,0,0;(1,2)}_{3,-2,-1,0,0,0}e_1\cos(3\lambda_2-2\lambda_1-\varpi_1) + C^{0,0,0,0;(1,2)}_{3,-2,0,-1,0,0} e_2\cos(3\lambda_2-2\lambda_1-\varpi_2)\right)\]

but expressed in the canonical variables used by celmech. 2

2.2. Integration

Now that we have a Hamiltonain model, we’ll integrate it and compare the results to direct \(N\)-body. First, we’ll set up some preliminary python dictionaries and arrays to hold the results of both integrations.

# Here we define the times at which we'll get simulation outputs
Nout = 150
times = np.linspace(0 , 3e3, Nout) * sim.particles[1].P

# These are the quantites we'll track in our rebound and celmech integrations
keys = ['l1','l2','pomega1','pomega2','e1','e2','a1','a2']

# These dictionaries will hold our results
rebound_results= {key:np.zeros(Nout) for key in keys}
celmech_results= {key:np.zeros(Nout) for key in keys}

# These are the lists of particles in both simulations
# for which we'll save quantities.
rb_particles = sim.particles
cm_particles = pvars.particles

The PoincareHamiltonian class inherits the method integrate() that can be used to evolve the system forward in much the same way as rebound’s rebound.Simulation.integrate() method. Below is the main integration loop where we’ll integrate our system and store the results:

for i,t in enumerate(times):
    sim.integrate(t) # advance N-body
    pham.integrate(t) # advance celmech
    for j,p_rb,p_cm in zip([1,2],rb_particles[1:],cm_particles[1:]):
        # store N-body results
        rebound_results["l{}".format(j)][i] = p_rb.l
        rebound_results["pomega{}".format(j)][i] = p_rb.pomega
        rebound_results["e{}".format(j)][i] = p_rb.e
        rebound_results["a{}".format(j)][i] = p_rb.a

        # store celmech results
        celmech_results["l{}".format(j)][i] = p_cm.l
        celmech_results["pomega{}".format(j)][i] = p_cm.pomega
        celmech_results["e{}".format(j)][i] = p_cm.e
        celmech_results["a{}".format(j)][i] = p_cm.a

Finally, we’ll plot the simulation results in order to compare them:

# First, we compute resonant angles for both sets of results
for d in [celmech_results,rebound_results]:
    d['theta1'] = np.mod(3 * d['l2'] - 2 * d['l1'] - d['pomega1'],2*np.pi)
    d['theta2'] = np.mod(3 * d['l2'] - 2 * d['l1'] - d['pomega2'],2*np.pi)

# Now we'll create a figure...
import matplotlib.pyplot as plt
fig,ax = plt.subplots(3,2,sharex = True,figsize = (12,8))
for i,q in enumerate(['theta','e','a']):
    for j in range(2):
        key = "{:s}{:d}".format(q,j+1)
        ax[i,j].plot(times,rebound_results[key],'k.',label='$N$-body')
        ax[i,j].plot(times,celmech_results[key],'r.',label='celmech')
        ax[i,j].set_ylabel(key,fontsize=15)
        ax[i,j].legend(loc='upper left')

#... and make it pretty
ax[0,0].set_ylim(0,2*np.pi);
ax[0,1].set_ylim(0,2*np.pi);
ax[2,0].set_xlabel(r"$t/P_1$",fontsize=15);
ax[2,1].set_xlabel(r"$t/P_1$",fontsize=15);

This should produce a figure that looks something like this:

_images/quickstart_example_plot.png

Not too bad! Our celmech model reproduces the libration amplitudes and frequencies observed in the \(N\)-body results quite successfully.

2.3. Next steps

For a more detailed example (discussed in Sec. 2 of the paper) that incorporates additional corrections to better approximate the N-body evolution, see AddingDisturbingFunctionTerms.ipynb.

You can also check out the rest of the documentation for many other features. Another good starting point for exploration are the numerous Jupyter notebook examples available on GitHub.

1

The precise definitions of the orbital elements and mass parameters \(\mu_i,M_i\) depend on the adopted coordinate system. By default celmech uses canonical heliocentric coordinates.

2

The \(C\) coefficients used by celmech are defined in Disturbing Function. For those familiar with the notation of Murray & Dermott (1999), \(C^{0,0,0,0;(1,2)}_{3,-2,-1,0,0,0} = f_{27}(\alpha)\) and \(C^{0,0,0,0;(1,2)}_{3,-2,0,-1,0,0} = f_{31}(\alpha)\) evaluated at \(\alpha\approx (2/3)^{2/3}\).