11. Miscellaneous

A number of useful odds and ends are gathered into the celmech.miscellaneous module. These include:

1. Frequency analysis: The frequency_modified_fourier_transform() function provides a wrapper for the FMFT algorithm of Šidlichovský & Nesvorný (1996). A jupyter notebook example is available here.

2. AMD: Functions related to angular momentum deficit (AMD) and the analytic stability criteria (see Laskar & Petit (2017)):

• compute_AMD()

• AMD_stable_Q()

• AMD_stability_coefficients()

• critical_relative_AMD()

• critical_relative_AMD_resonance_overlap()

3. The function holman_weigert_stability_boundary() which computes a stability criterion for circum-binary planets based on Holman & Wiegert (1999)

4. The functions sk() and Dsk() which give approximate disturbing function coefficient amplitudes as and their derivatives with respect to eccentricity. See Hadden & Lithwick (2018)

5. Convenience functions get_symbol() and get_symbol0() for generating sympy symbols from LaTeX strings.

6. The function truncated_expansion() for expanding and truncating sympy expressions after assigning variables (possibly different) orders in a book-keeping parameter.

11.1. API

celmech.miscellaneous.AMD_stability_coefficients(sim, overlap=False)

Compute AMD stability coefficients of the successive adjacent planet pairs of a planetary system.

A planet pair’s AMD stability coefficicent is defined as the total planetary system’s AMD divided by the critical AMD required for the pair’s orbits to cross. (Equation 58 of Laskar & Petit (2017))

Parameters:

• sim (rebound.Simulation) – Simulation object to copmute AMD coefficients for.

• overlap (bool, optional) – If True, planet pairs’ critical AMD values are computed as the critical AMD value for resonance overlap. By default the critical values are computed as the value required for orbit crossing.

Returns:

Values of \(\beta = \frac{C}{\Lambda'C_c}\) for planet pairs.

Return type:

ndarray

celmech.miscellaneous.AMD_stable_Q(sim)

Test the AMD-stability of a planetary system. Returns True if a planetary system is AMD-stable and False if not.

Parameters:

sim (rebound.Simulation) – Simulation object to copmute stability criterion for.

Returns:

True if the sytem is AMD-stable, otherwise False.

Return type:

bool

celmech.miscellaneous.Dsk(k, y, tol=1.49e-08, rtol=1.49e-08, maxiter=50, miniter=1)

Derivative of disturibing function coefficient s_k with respect to argument y. Coefficients are described in Hadden & Lithwick (2018)

Quadrature routine based on scipy.quadrature.

Parameters:

• k (int) – Order of resonance

• y (float) – e / ecross; must be y<1

• tol (float, optional) – Control absolute and relative tolerance of integration. Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

• rtol (float, optional) – Control absolute and relative tolerance of integration. Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

• maxiter (int, optional) – Maximum order of Gaussian quadrature.

• miniter (int, optional) – Minimum order of Gaussian quadrature

Returns:

val – Gaussian quadrature approximation of s_k(y)

Return type:

float

celmech.miscellaneous.EulerMatrix(Omega, inc, omega)

The Euler 3D rotation matrix for Euler angles (Omega,inc,omega). The (3,1,3) convention is followed.

Parameters:

• Omega (float) – First Euler angle (ascending node)

• inc (float) – Second Euler angle (inclination)

• omega (float) – Third Euler angle (argument of periapsis)

Returns:

3 x 3 rotation matrix

Return type:

numpy.array

celmech.miscellaneous.compute_AMD(sim)

Compute total AMD of a planetary system.

The angular momentum deficit (AMD) of a planetary system is the difference between the angular momentum of a hypothetical system with the same masses and semi-major axes but with circular, coplanar orbits and the actual angular momentum of a planetary system. It is a conserved quantity of the purely secular dynamics of a system.

Parameters:

sim (rebound.Simulation) – A REBOUND simulation of a planetary system.

Returns:

AMD – The value of the systems angular momentum deficit.

Return type:

float

celmech.miscellaneous.critical_relative_AMD(alpha, gamma)

The critical value of ‘relative AMD’, \({\cal C} = C/\Lambda_\mathrm{out}\), of a planet pair above which intersecting orbits are allowed.

See Equation 29 of Laskar & Petit (2017)

Parameters:

• alpha (float) – The semi-major axis ratio, \(\alpha=a_\mathrm{in}/a_\mathrm{out}\) of the planet pair.

• gamma (float) – The mass ratio of the planet pair, \(\gamma = m_\mathrm{in}/m_\mathrm{out}\).

Returns:

Ccrit – The value of the the critical AMD (\(C_c(\alpha,\gamma)\) in the notation of Laskar & Petit (2017)

Return type:

float

celmech.miscellaneous.critical_relative_AMD_resonance_overlap(alpha, gamma, mutot)

The critical value of ‘relative AMD’, \({\cal C} = C/\Lambda_\mathrm{out}\), of a planet pair above which resonance overlap can occur based on the resonance overlap criterion of Hadden & Lithwick (2018)

Parameters:

• alpha (float) – The semi-major axis ratio, \(\alpha=a_\mathrm{in}/a_\mathrm{out}\) of the planet pair.

• gamma (float) – The mass ratio of the planet pair, \(\gamma = m_\mathrm{in}/m_\mathrm{out}\).

• mutot (float) – The total mass of the planet pair relative to the star, i.e., \((\mu_\mathrm{in} + \mu_\mathrm{out}) / M_*\)

Returns:

Ccrit – The value of the the critical AMD (\(C_c(\alpha,\gamma)\) in the notation of Laskar & Petit (2017)

Return type:

float

celmech.miscellaneous.frequency_modified_fourier_transform(time, z, Nfreq, method_flag=3, min_freq=None, max_freq=None)

Apply the frequency-modified Fourier transfrorm algorithm of Šidlichovský & Nesvorný (1996) to a time series to determine the series’ principle Fourier modes. This function simply proivdes a wrapper to to C implementation written by D. Nesvorný available at www-n.oca.eu/nesvorny/programs.html.

Parameters:

• time (ndarray, shape (N,)) – Times of input data values.

• z (complex ndarray, shape (N,)) – Input data time series in the form.

• Nfreq (int) – Number of Fourier modes to determine.

• method_flag (int) –

  The FMFT algorithm

  Basic Fourier Transform algorithm if flag = 0; not implemented Modified Fourier Transform if flag = 1; Frequency Modified Fourier Transform if flag = 2; FMFT with additional non-linear correction if flag = 3

celmech.miscellaneous.getOmegaMatrix(n)

Get the 2n x 2n skew-symmetric block matrix

\[\begin{split}\begin{pmatrix} 0 & I_n \\ -I_n & 0 \end{pmatrix},\end{split}\]

where \(I_n\) is the \(n \times n\) identity matrix, that appears in Hamilton’s equations.

Parameters:

n (int) – Determines matrix dimension

Return type:

numpy.array

celmech.miscellaneous.get_symbol(latex, subscript=None, **kwargs)

Get a sympy sympy based on an input LaTeX string. Valid keyword arguments for the function sympy.symbols can also be passed.

Parameters:

• latex (string) – LaTeX expression to render as a sympy symbol

• subscript (string or int, optional) – A subscript for the sympy symbol

Return type:

sympy symbol

celmech.miscellaneous.get_symbol0(latex, subscript=None, **kwargs)

Same as get_symbol(), but appends a “0” to the subscript.

celmech.miscellaneous.holman_weigert_stability_boundary(mu, e, Ptype=True)

Compute the critical semi-major axis represnting an approximate stability boundary for circumbinary planets in P- or S-type orbits. Formulas for critical semi-major axes are taken from Holman & Wiegert (1999)

Parameters:

• mu (float) –

  The mass-ratio of the binary,

  \[\mu = \frac{m_B}{m_A+m_B}\]

  where math:m_A and math:m_B are the component masses of the binary.

• e (float) – The eccentricity of the binary.

• Ptype (bool, optional) – If True (default) orbit is assumed to be a P-type circumbinary orbit. If False, a S-type circum-primary/secondary orbit is considered.

Returns:

aC – The critical semi-major axis marking the stability boundary

Return type:

float

celmech.miscellaneous.levin_method_integrate(fvec_fn, wvec_fn, Amtrx_fns, a, b, N=32)

Evlauate integrals of the form

\[I(a,b) = \int_{a}^{b} \vec{f}(x)\cdot\vec{w}(x) dx\]

where the functions \(\vec{w}(x)\) satisfy a linear differential equation of the form \(\frac{d}{dx}\vec{w}(x) = A(x) \cdot \vec{w}(x)\). Evaluation is done using Levin’s method.

Parameters:

• fvec_fn (function) – A function that returns the vector \(\vec{f}(x)\).

• wvec_fn (function) – A function that returns the vector \(\vec{w}(x)\).

• Amtrx_fns (list of functions) – A list of functions giving the entries of matrix \(A(x)\) appearing in the differential equation obeyed by \(\vec{w}(x)\). Input should use nested lists to match the structure of the matrix.

• a (float) – Lower integration limit

• b (float) – Upper integration limit

• N (int) – Number of quadrature points to use

Return type:

float

celmech.miscellaneous.levin_method_integrate_adaptive(fvec_fn, wvec_fn, Amtrx_fns, a, b, N0=32, Nmax=128, rtol=1.49e-08, atol=1.49e-08)

Evlauate integrals of the form

\[I(a,b) = \int_{a}^{b} \vec{f}(x)\cdot\vec{w}(x) dx\]

where the functions \(\vec{w}(x)\) satisfy a linear differential equation of the form \(\frac{d}{dx}\vec{w}(x) = A(x) \cdot \vec{w}(x)\). Evaluation is done using Levin’s method.

Method is applied adaptively with increasing number of quadrature points until the estimated error, \(\delta\) satisfies \(\delta < \epsilon_\mathrm{rel}|I(a,b)| + \epsilon_\mathrm{abs}\).

Parameters:

• fvec_fn (function) – A function that returns the vector \(\vec{f}(x)\).

• wvec_fn (function) – A function that returns the vector \(\vec{w}(x)\).

• Amtrx_fns (list of functions) – A list of functions giving the entries of matrix \(A(x)\) appearing in the differential equation obeyed by \(\vec{w}(x)\). Input should use nested lists to match the structure of the matrix.

• a (float) – Lower integration limit

• b (float) – Upper integration limit

• N0 (int) – Initial number of Gauss-Lobatto quadrature points to use.

• Nmax (int) – Maximum number of quadrature points to use.

• rtol (float) – Relative tolerance

• atol (float) – Absolute tolerance

Return type:

float

celmech.miscellaneous.linking_l(orbit1, orbit2)

Computes the linking coefficient defined by Kholshevnikov and Vassiliev (1999) for a pair of Keplerian orbits.

Parameters:

• orbit1 (rebound.Orbit) – First orbit

• orbit2 (rebound.Orbit) – Second orbit

Returns:

The linking coefficient, \(l_1\), defined by Equation (1) of Kholshevnikov and Vassiliev (1999)

Return type:

float

celmech.miscellaneous.linsolve(A, y)

Solve linear system of equations

\[A \cdot x = y\]

for y.

celmech.miscellaneous.poisson_bracket(f, g, re_varslist, complex_varslist)

Calculate the Poisson bracket

\[[f,g] = \sum_{i=1}^N \frac{\partial f}{\partial q_i} \frac{\partial g}{\partial p_i} - \frac{\partial f}{\partial p_i} \frac{\partial g}{\partial q_i} - i \sum_{j=1}^{M} \frac{\partial f}{\partial z_j} \frac{\partial g}{\partial \bar{z}_j} - \frac{\partial f}{\partial \bar{z}_j} \frac{\partial g}{\partial {z}_i}\]

where re_varslist is \(=(q_1,...,q_N,p_1,...,p_N)\) and complex_varslist is \(=(x_1,...,x_M,\bar{x}_1,...,\bar{x}_M)\).

Parameters:

• f (sympy expression) – Function appearing in Poisson bracket.

• g (sympy expression) – Other function appearing in Poisson bracket.

• re_varslist (list of sympy symbols) – List of real canonical variables in the form \((q_1,...,q_N,p_1,...,p_N)\)

• complex_varslist (list of sympy symbols) – List of complex canonical variables in the form \((x_1,...,x_M,\bar{x}_1,...,\bar{x}_M)\)

Return type:

sympy expression

celmech.miscellaneous.sk(k, y, tol=1.49e-08, rtol=1.49e-08, maxiter=50, miniter=1)

Approximate disturibing function coefficient described in Hadden & Lithwick (2018)

Quadrature routine based on scipy.quadrature.

Parameters:

• k (int) – Order of resonance

• y (float) – e / ecross; must be y<1

• tol (float, optional) – Control absolute and relative tolerance of integration. Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

• rtol (float, optional) – Control absolute and relative tolerance of integration. Iteration stops when error between last two iterates is less than tol OR the relative change is less than rtol.

• maxiter (int, optional) – Maximum order of Gaussian quadrature.

• miniter (int, optional) – Minimum order of Gaussian quadrature

Returns:

val – Gaussian quadrature approximation of s_k(y)

Return type:

float

celmech.miscellaneous.sundman_convergence(sim)

Determine whether the expansion of the disturbing function in eccentricity converges for all pairs of orbits in a planetary system based on the Sundman criterion. See Ferraz-Mello (1994) for details.

Parameters:

sim (rebound.Simulation) – A REBOUND simulation of a planetary system.

Returns:

True if the expansion of the disturbing function in eccentricity converges for all pairs of orbits in the system, otherwise False.

Return type:

bool

celmech.miscellaneous.sundman_factor(p1, p2)

Determine whether the expansion of the disturbing function in eccentricity converges for a pair of orbits based on the Sundman criterion. See Ferraz-Mello (1994) for details.

Parameters:

• p1 (rebound.Particle or celmech.PoincareParticle) – One of the particles in the pair for which to compute the Sundman factor.

• p2 (rebound.Particle or celmech.PoincareParticle) – The other particle in the pair for which to compute the Sundman factor.

Returns:

Value of the Sundman factor for the pair of orbits. If this value is less than 1, the expansion of the disturbing function in eccentricity converges. If it is greater than 1, the expansion does not converge.

Return type:

float

celmech.miscellaneous.truncated_expansion(exprn, order_rules, max_order)

Expand a sympy expression up to a maximum order in a small book-keeping parameter after assigning variables appearing in the expression a given order using the order_rules argument.

Parameters:

• exprn (sympy expression) – The original expression from which to calculate expansion.

• order_rules (dict) –

  A dictionary specifying what order various variables should be assumed to have in the book-keeping parameter. Each key-value pair {n:[x_1,x_2,..,x_m]} in order_rules specifies that a set of variables

  \[(x_1,...,x_m) \sim \mathcal{O}(\epsilon^n)\]

  where \(\epsilon\) is the book-keeping parameter.

• max_order (int) – The order at which the resulting series expansion in the book-keeping parameter \(\epsilon\) should be truncated.

Return type:

sympy expression