12. 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)):
The function
holman_weigert_stability_boundary()
which computes a stability criterion for circum-binary planets based on Holman & Wiegert (1999)The functions
sk()
andDsk()
which give approximate disturbing function coefficient amplitudes as and their derivatives with respect to eccentricity. See Hadden & Lithwick (2018)Convenience functions
get_symbol()
andget_symbol0()
for generating sympy symbols from LaTeX strings.The function
truncated_expansion()
for expanding and truncating sympy expressions after assigning variables (possibly different) orders in a book-keeping parameter.
12.1. API¶
- celmech.miscellaneous.AMD_stability_coefficients(sim, overlap=False)[source]¶
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)[source]¶
Test the AMD-stability of a planetary system. Returns
True
if a planetary system is AMD-stable andFalse
if not.- Parameters
sim (
rebound.Simulation
) – Simulation object to copmute stability criterion for.- Returns
True
if the sytem is AMD-stable, otherwiseFalse
.- Return type
- celmech.miscellaneous.Dsk(k, y, tol=1.49e-08, rtol=1.49e-08, maxiter=50, miniter=1)[source]¶
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
- celmech.miscellaneous.compute_AMD(sim)[source]¶
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
- celmech.miscellaneous.critical_relative_AMD(alpha, gamma)[source]¶
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
- Returns
Ccrit – The value of the the critical AMD (\(C_c(\alpha,\gamma)\) in the notation of Laskar & Petit (2017)
- Return type
- celmech.miscellaneous.critical_relative_AMD_resonance_overlap(alpha, gamma, mutot)[source]¶
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
- celmech.miscellaneous.frequency_modified_fourier_transform(time, z, Nfreq, method_flag=3, min_freq=None, max_freq=None)[source]¶
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)[source]¶
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)[source]¶
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)[source]¶
Same as
get_symbol()
, but appends a “0” to the subscript.
- celmech.miscellaneous.holman_weigert_stability_boundary(mu, e, Ptype=True)[source]¶
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. IfFalse
, a S-type circum-primary/secondary orbit is considered.
- Returns
aC – The critical semi-major axis marking the stability boundary
- Return type
- celmech.miscellaneous.levin_method_integrate(fvec_fn, wvec_fn, Amtrx_fns, a, b, N=32)[source]¶
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
- 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)[source]¶
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
- celmech.miscellaneous.linsolve(A, y)[source]¶
Solve linear system of equations
\[A \cdot x = y\]for y.
- celmech.miscellaneous.poisson_bracket(f, g, re_varslist, complex_varslist)[source]¶
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)\) andcomplex_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)[source]¶
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
- celmech.miscellaneous.truncated_expansion(exprn, order_rules, max_order)[source]¶
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]}
inorder_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