3. Poincare Module

3.1. Introduction

The Poincare module is used to represent planetary systems in terms of Poincare canonical variables. These canonical variables are related to the Delauney elliptic variables,

\[\begin{split}\begin{eqnarray} \Lambda_i &=& \mu_i\sqrt{G{ M}_ia_i}\nonumber\\ \Gamma_i &=& \mu_i\sqrt{G{ M}_ia_i}(1-\sqrt{1-e_i^2})\nonumber\\ Q_i &=& \mu_i\sqrt{{M}_ia_i}\sqrt{1-e_i^2}(1-\cos(I_i)) \end{eqnarray}\end{split}\]

where \(\mu_i\sim m_i\) and \(M_i\sim M_*\). (The exact definitions of \(\mu_i\sim m_i\) and \(M_i\sim M_*\) are described in greater detail below.) The associated canonical conjugate coordinates are the mean longitudes \(\lambda_i, \gamma_i=-\varpi_i,~\mathrm{and}~q_i=-\Omega_i\) where \(\varpi_i\) is the longitude of periapse and \(\Omega_i\) is the longitude of ascending node. Instead of the variables \((q_i,Q_i)~\mathrm{and}~(\gamma_i,\Gamma_i)\), celmech formulates equations of motion in terms of ‘Cartesian-style’ canonical coordinate-momentum pairs

\[\begin{split}\begin{eqnarray} (\eta_i,\kappa_i)&=&\sqrt{2\Gamma_i} \times (\sin\gamma_i,\cos\gamma_i)\\ (\sigma_i,\rho_i)&=&\sqrt{2Q_i} \times (\sin q_i,\cos q_i) \end{eqnarray}\end{split}\]

which have the advantage of being well-defined for \(e_i\rightarrow 0\) and \(I_i\rightarrow 0\).

celmech uses the Poincare class to represent a set of canonical Poincare variables for a planetary system. The Poincare.from_Simulation and Poincare.to_Simulation methods provide easy integration with REBOUND.

3.2. Constructing A Hamiltonian

The PoincareHamiltonian class is used to construct Hamiltonians or normal forms governing the dynamics of a system as well as integrate the corresponding equations of motion. A PoincareHamiltonian instance is initialized by passing the Poincare instance whose dynamics will be governed by the resulting Hamiltonian.

Upon intialization, a PoincareHamiltonian consists only of the Keplerian components of the Hamiltonian, i.e.,

\[H = -\sum_{i=1}^{N}\frac{G^2M_i^2\mu^3}{2\Lambda_i^2}\]

Individual disturbing funtion terms from planet pairs’ interactions can then be added with the add_cosine_term() method. For example, to add the leading order term \(\propto \cos(3\lambda_2 - \lambda_1 - \varpi_1-\varpi_2)\) to the Hamiltonian represented by a PoincareHamiltonian instance, pham, you would call

pham.add_cosine_term([3,-1,-1,-1,0,0],indexIn=1,indexOut=2)

More details about the disturbing function expansion are described in the Disturbing Function section. The PoincareHamiltonian provides additional methods that can be used to conveniently add multiple disturbing function terms at once, including add_MMR_terms() and add_secular_terms().

3.3. Canonical Heliocentric Coordinates

celmech uses canonical heliocentric coordinates in order to formulate the equations of motion.

3.3.1. Cartesian Heliocentric Coordinates

If \(\pmb{u}_i\) are the Cartesian position vectors of an \(N\)-body system in an inertial frame and \(\tilde{\pmb{u}}_i=m_i\frac{d}{dt}\pmb{u}_i\) are their corresponding conjugate momenta, the canonical heliocentric coordinates \(\pmb{r}_i\) are given by

\[\begin{split}\begin{pmatrix} \pmb{r}_0\\ \pmb{r}_1\\ \pmb{r}_2\\ \vdots \\ \pmb{r}_{N-1} \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ -1 & 1 & 0 & \cdots & 0 \\ -1 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0\\ -1 & 0 & 0 & \cdots & 1 \end{pmatrix} \cdot \begin{pmatrix} \pmb{u}_0\\ \pmb{u}_1\\ \pmb{u}_2\\ \vdots \\ \pmb{u}_{N-1} \end{pmatrix}\end{split}\]

and their conjugate momenta are

\[\begin{split}\begin{pmatrix} \tilde{\pmb{r}}_0\\ \tilde{\pmb{r}}_1\\ \tilde{\pmb{r}}_2\\ \vdots \\ \tilde{\pmb{r}}_{N-1} \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 & \cdots & 1 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0\\ 0 & 0 & 0 & \cdots & 1 \end{pmatrix} \cdot \begin{pmatrix} \tilde{\pmb{u}}_0\\ \tilde{\pmb{u}}_1\\ \tilde{\pmb{u}}_2\\ \vdots \\ \tilde{\pmb{u}}_{N-1} \end{pmatrix}\end{split}\]

The \(N\)-body Hamiltonian is independent of \(\pmb{r}_0\) and, consequently \(\tilde{\pmb{r}}_0\), i.e., the total momentum of the system, is conserved. In the new coordinates, the Hamiltonian becomes

\[\begin{split}\begin{eqnarray} H &=& \sum_{i=1}^{N-1}H_{\mathrm{kep},i} + \sum_{i=1}^{N-1} \sum_{j=1}^{i} H_{\mathrm{int.}}^{(i,j)} \\ H_{\mathrm{kep},i} &=& \frac{1}{2}\frac{\tilde{r}_i^2}{\mu_i} - \frac{GM_i\mu_i}{r_i} \\ H_{\mathrm{int.}}^{(i,j)} &=& -\frac{Gm_im_j}{|\pmb{r}_i - \pmb{r}_j|} + \frac{\tilde{\pmb{r}}_i \cdot \tilde{\pmb{r}}_j }{M_*} \end{eqnarray}\end{split}\]

where the parameters

\[\begin{split}{\mu}_i = \frac{m_iM_*}{M_* + m_i}\\ {M}_i = {M_* + m_i}~.\end{split}\]

are the same as those that appear in the definitions of celmech’s canonical variables.

Above, each \(H_{\mathrm{kep},i}\) is the Hamiltonian of a two-body probelm with reduced mass \(\mu_i\). We can therefore rewrite it in terms of the canonical Delauney variables defined above as

\[H_{\mathrm{kep},i} = -\frac{G^2M_i^2\mu^3}{2\Lambda_i^2}\]

3.3.2. Heliocentric Orbital Elements

The Delauney variables, derived via a canonical transformation from the canonical heliocentric coordinates, are defined in terms of a set of orbital elements. A set of six Keplerian orbital elements , \((a_i,e_i,I_i,M_i,\varpi_i,\Omega_i)\), define by a mapping to Cartesian positions and velocities. In the case of the elements appearing in the Delauney variables, the elements specify the heliocentric position, \(\pmb{r}_i\) and the “velocity” \(\tilde{\pmb{r}}_i/\mu_i = \frac{M_* + m_i}{M_*}\frac{d \pmb{u}_i}{dt}\). While this does not correspond to any physical velocity in the system, it ensures that the transformation from the coordinate-momentum pairs \((\pmb{r},\tilde{\pmb{r}}_i)\) to Delauney variables is canonical.

celmech provides the functions nbody_simulation_utilities.get_canonical_heliocentric_orbits to compute these ‘canonical heliocentric’ elements from a REBOUND simulation along with nbody_simulation_utilities.add_canonical_heliocentric_elements_particle to add particles to a REOBOUND simulation by specifying the particles orbit in terms of these elements.

For an example demonstrating the use of these functions, see CoordinatesConvertingToFromNbody.ipynb.

3.4. API

class celmech.poincare.Poincare(G, poincareparticles, coordinates='canonical heliocentric', t=0)[source]

A class representing a collection of Poincare particles constituting a planetary system.

property N

Return total number of bodies, including central star (analogous to REBOUND sim.N)

__init__(G, poincareparticles, coordinates='canonical heliocentric', t=0)[source]
Parameters

  • qp_vars (list of symbols) – List of symbols to use for the canonical coordiantes and momenta. The list should be of length 2*N for an integer N with entries 0,…,N-1 representing the canonical coordiante variables and entries N,…,2*N-1 representing the corresponding conjugate momenta.

  • values (array-like) – The numerical values assigned to the canonical coordiantes and momenta.

  • t (float, optional) – The current time of system represented by the phase-space state. Default value is 0.

classmethod from_Simulation(sim, coordinates='canonical heliocentric')[source]

Convert REBOUND Simulation to Poincare object, using specified canonical coordinates. Assumes the dominant mass is sim.particles[0].

Parameters

  • sim (rebound.Simulation) – Simulation to convert.

  • coordinates (str) – Specifices the canonical coordinate system. This determines the appropriate definitions of mu and M. Options: ‘canonical heliocentric’ (default): canonical heliocentric coordinates in the COM frame e.g. Laskar & Robutel 1995 ‘democratic heliocentric’: e.g. Duncan et al. 1998

Return type

Poincare object

to_Simulation()[source]

Convert Poincare object to a REBOUND simulation in COM frame.

Returns

sim

Return type

rebound.Simulation

class celmech.poincare.PoincareHamiltonian(pvars)[source]

A class representing the Hamiltonian governing the dynamical evolution of a system of particles, stored as a celmech.poincare.Poincare instance.

H

Symbolic expression for the Hamiltonian.

Type

sympy expression

NH

Symbolic expression for the Hamiltonian with numerical values of parameters substituted where applicable.

Type

sympy expression

N

Number of particles

Type

int

particles

List of :class:`celmech.poincare.PoincareParticle`s making up the system.

Type

list

state

A set of Poincare variables to which transformations are applied.

Type

celmech.poincare.Poincare

__init__(pvars)[source]
Parameters

  • H (sympy expression) – Hamiltonian made up only of sympy symbols in state.qp_pairs and keys in H_params

  • H_params (dict) – dictionary from sympy symbols for the constant parameters in H to their value

  • state (PhaseSpaceState) – Object for holding the dynamical state.

  • above (In addition to the) –

  • objects (one needs to write 2 methods to map between the two) –

  • state_to_list(self (def) – returns a list of values from state in the same order as pqpairs e.g. [P1,Q1,P2,Q2]

  • state) – returns a list of values from state in the same order as pqpairs e.g. [P1,Q1,P2,Q2]

  • update_state_from_list(self (def) – updates state object from a list of values y for the variables in the same order as pqpairs and integrator time ‘t’

  • state – updates state object from a list of values y for the variables in the same order as pqpairs and integrator time ‘t’

  • y – updates state object from a list of values y for the variables in the same order as pqpairs and integrator time ‘t’

  • t) – updates state object from a list of values y for the variables in the same order as pqpairs and integrator time ‘t’

add_Hkep_term(H, index)[source]

Add the Keplerian component of the Hamiltonian for planet ‘’.

add_MMR_terms(p, q, max_order=None, l_max=0, indexIn=1, indexOut=2, eccentricities=True, inclinations=True)[source]

Add all eccentricity and inclination disturbing function terms associated with a p:p-q mean motion resonance up to a given order.

Parameters

  • p (int) –

    Coefficient of lambdaOut in resonant argument

    j*lambdaOut - (j-k)*lambdaIn

  • q (int) – Order of the mean motion resonance.

  • max_order (int) – Maximum order of terms to add.

  • l_max (int, optional) – Maximum degree of expansion in \(\delta = (\Lambda-\Lambda_0)/\Lambda_0\) to include in cosine coefficients. Default is 0.

  • indexIn (int) – Index of inner planet.

  • indexOut (int) – Index of outer planet.

  • eccentricities (bool, optional) – By default, includes all eccentricity terms. Can set to False to exclude any eccentricity terms (e.g., fixed circular orbits).

  • inclinations (bool, optional) – By default, includes all inclination terms. Can set to False to exclude any inclination terms (e.g., co-planar systems).

add_cosine_term(k_vec, max_order=None, l_max=0, nu_vecs=None, l_vecs=None, indexIn=1, indexOut=2, eccentricities=True, inclinations=True)[source]

Add a specific cosine term to the disturbing function up to an optional max_order in the eccentricities and inclinations.

Parameters

  • k_vec (Tuple of 6 ints) – Coefficients (k1, k2, k3, k4, k5, k6) for a cosine term of the form cos(k1*lambdaOut + k2*lambdaIn + k3*pomegaIn + k4*pomegaOut + k5*OmegaIn + k6*OmegaOut) Note the ks must sum to zero, and by convention we write lambdaOut first, e.g., (3, -2, -1, 0, 0, 0) is cos(3*lambdaOut - 2*lambdaIn - pomega1)

  • max_order (int, optional) – Maximum order to go to in the (combined) eccentricities and inclinations. Can specify either max_order OR nu_vecs (see below), but not both (most users will use max_order). If neither are passed, max_order defaults to the leading order of the cosine term specified by k_vec.

  • l_max (int, optional) – Maximum degree of expansion in \(\delta = (\Lambda-\Lambda_0)/\Lambda_0\) to include in cosine coefficients. Default is 0. Can specify either l_max OR l_vecs (see below), but not both (most users will use l_max).

  • nu_vecs (List of 4-tuples, optional) – A list of the specific combinations of nu indices to include for the cosine term coefficient, e.g., [(0, 0, 0, 0), (1, 0, 0, 0), …] See paper and examples for definition and use. Can specify either max_order OR nu_vecs, but not both (max_order makes more sense for most use cases).

  • l_vecs (List of 2-tuples, optional) – A list of the specific combinations of l indices to include for the cosine term coefficient, e.g., [(0, 0), (1, 0), (2, 0), …] See paper and examples for definition and use. Can specify either l_max OR l_vecs, but not both (l_max makes more sense for most use cases). One use case for passing particular l_vecs is if one body is massless, so the other Lambda doesn’t vary (and doesn’t need expansion)

  • indexIn (int, optional) – Index of inner planet. Default 1.

  • indexOut (int, optional) – Index of outer planet. Default 2.

  • eccentricities (bool, optional) – By default, includes all eccentricity terms. Can set to False to exclude any eccentricity terms (e.g., fixed circular orbits).

  • inclinations (bool, optional) – By default, includes all inclination terms. Can set to False to exclude any inclination terms (e.g., co-planar systems).

add_gr_potential_terms(c, max_e_order=None, particles='all')[source]

Add Hamiltonian terms that capture the orbital precession caused by general relativity by adding a potential term of the form

\[\phi_\mathrm{GR} = \frac{3G^2M_*^2}{c^2a^2\sqrt{1-e^2}}\]
Parameters

  • c (float) – The speed of light in the appropriate simulation units. harmonic.

  • max_e_order (int, optional) – Maximum order of expansion in eccentricity. By default, the value is set to ‘None’ and no expansion in eccentricity is done.

  • particles (list, optional) – Which particle numbers to add \(J_2\) terms for. Default is set to all particles.

add_orbit_average_J2_terms(J2, Rin, max_ei_order=None, max_delta_order=None, particles='all', **kwargs)[source]

Add Hamiltonian terms that capture the orbit-averaged effect of a central body’s oblateness parameterized by the \(J_2\) gravitational harmonic.

Parameters

  • J2 (float) – The value of the central body’s J2 gravitational harmonic.

  • Rin (float) – The central body’s radius

  • max_ei_order (int, optional) – Maximum order of expansion in eccentricity and inclination. By default, the value is set to ‘None’ and no expansion in eccentricity and inclinaion is done.

  • max_delta_order (int, optional) – Maxmimum order in \(\delta =(\Lambda-\Lambda_0)/\Lambda_0). Default is 'None' and dependence on :math:\) is exact.

  • particles (list, optional) – Which particle numbers to add \(J_2\) terms for. Default is set to all particles.

add_secular_terms(min_order=2, max_order=2, l_max=0, indexIn=1, indexOut=2, eccentricities=True, inclinations=True)[source]

Add all eccentricity and inclination secular terms in the disturbing function from min_order to max_order in the eccentricities and inclinations.

Parameters

  • min_order (int, optional) – Minimum order terms in the eccentricities and inclinations to add. Defaults to 2.

  • max_order (int, optional) – Maximum order terms in the eccentricities and inclinations to add. Defaults to 2.

  • l_max (int, optional) – Maximum degree of expansion in \(\delta = (\Lambda-\Lambda_0)/\Lambda_0\) to include in cosine coefficients. Default is 0. to include in cosine coefficients. Default is 0.

  • indexIn (int) – Index of inner planet.

  • indexOut (int) – Index of outer planet.

  • eccentricities (bool, optional) – By default, includes all eccentricity terms. Can set to False to exclude any eccentricity terms (e.g., fixed circular orbits).

  • inclinations (bool, optional) – By default, includes all inclination terms. Can set to False to exclude any inclination terms (e.g., co-planar systems).

property df

Pretty print the terms currently included in the disturbing function in a jupyter notebook.

property df_latex

Print the latex for the terms currently included in the disturbing function (e.g., to copy-paste into a paper)

property df_raw_latex

Get the raw latex string (r”…”) for the terms currently included in the disturbing function. This property is likely not what you’re looking for. Most users will either want df to pretty print df in a jupyter notebook, or df_latex for the actual latex (with single slashes) to copy-paste into a paper.

class celmech.poincare.PoincareParticle(coordinates='canonical heliocentric', G=1.0, m=None, Mstar=None, mu=None, M=None, sLambda=None, l=None, sGamma=None, gamma=None, sQ=None, q=None, Lambda=None, Gamma=None, Q=None, a=None, e=None, inc=None, pomega=None, Omega=None, skappa=None, seta=None, ssigma=None, srho=None, kappa=None, eta=None, sigma=None, rho=None)[source]

A class representing an individual member (star, planet, or test particle) of a planetary system. The appropriate value for mu and M depends on the adopted coordinate system and Kepler splitting (see e.g., Hernandez and Dehnen 2017 for a review and comparison). celmech supports canonical heliocentric coordinates (default) and democratic heliocentric coordinates.

Parameters

  • coordinates (str) – Specifices the canonical coordinate system. This determines the appropriate definitions of mu and M. Options: ‘canonical heliocentric’ (default): canonical heliocentric coordinates in the COM frame e.g. Laskar & Robutel 1995 ‘democratic heliocentric’: e.g. Duncan et al. 1998

  • G (float) – Gravitational constant (Default: 0)

  • m (float) – Physical mass of particle.

  • Mstar (float) – Physical mass of central body.

  • mu (float) – ‘Canonical’ mass of body. mu=reduced mass for canonical heliocentric coordinates (default) mu=m for democratic heliocentric coordinates.

  • M (float) – ‘Canonical’ central mass. M=Mstar+m for canonical heliocentric coordinates (default) M=Mstar for democratic heliocentric coordinates.

  • Lambda (float) – These variables specify the semimajor axis of the orbit. Can pass any of the three, but at least one must be specified.

  • sLambda (float) – These variables specify the semimajor axis of the orbit. Can pass any of the three, but at least one must be specified.

  • a (float) – These variables specify the semimajor axis of the orbit. Can pass any of the three, but at least one must be specified.

  • l (float) – Mean longitude of the orbit. If not passed, defaults to 0.

  • Gamma (float) – These variables specify the orbital eccentricity. Any one can be passed. If none passed, defaults to 0

  • sGamma (float) – These variables specify the orbital eccentricity. Any one can be passed. If none passed, defaults to 0

  • e (float) – These variables specify the orbital eccentricity. Any one can be passed. If none passed, defaults to 0

  • gamma (float) – These variables specify the pericenter orientation. Any one can be passed. If none passed, defaults to 0

  • pomega (float) – These variables specify the pericenter orientation. Any one can be passed. If none passed, defaults to 0

  • Q (float) – These variables specify the orbital inclination. Any one can be passed. If none passed, defaults to 0

  • sQ (float) – These variables specify the orbital inclination. Any one can be passed. If none passed, defaults to 0

  • inc (float) – These variables specify the orbital inclination. Any one can be passed. If none passed, defaults to 0

  • q (float) – These variables specify the node longitude. Any one can be passed. If none passed, defaults to 0

  • Omega (float) – These variables specify the node longitude. Any one can be passed. If none passed, defaults to 0

__init__(coordinates='canonical heliocentric', G=1.0, m=None, Mstar=None, mu=None, M=None, sLambda=None, l=None, sGamma=None, gamma=None, sQ=None, q=None, Lambda=None, Gamma=None, Q=None, a=None, e=None, inc=None, pomega=None, Omega=None, skappa=None, seta=None, ssigma=None, srho=None, kappa=None, eta=None, sigma=None, rho=None)[source]

We store Cartesian components of specific actions to support test particles

class celmech.poincare.PoincareParticles(poincare)[source]
__init__(poincare)[source]
celmech.poincare.get_re_im_components(x, y, k)[source]
Get the real and imaginary components of

(x + sgn(k) * i y)^|k|