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.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