Utilities

PNTerm{T,PNOrder,c⁻¹Exponent}

This object represents a single term in a PNExpansion. It has a single field: coeff, which is the coefficient of the term. The type parameter T is the type of the coefficient. The type parameter PNOrder is a half-integer (just as in PNSystems) representing the PN order of the expansion. And the type parameter c⁻¹Exponent is an integer representing the exponent of the PN expansion parameter $1/c$.

PNTerms can be multiplied and divided by scalars and exponentiated by integers, to produce another PNTerm. They can also be added to other PNTerms to produce a PNExpansion.

A simple way to define a PNTerm or a PNExpansion is to define the PN expansion parameter

c = PNExpansionParameter(pnsystem)

and use that naturally in formulas, as in

e = 1 + (v/c)^2 * (-ν/12 - 3//4) + (v/c)^4 * (-ν^2/24 + 19ν/8 - 27//8)

Any exponent higher than the desired PNOrder will be automatically set to zero.

Useful facts:

• v has order 1/c
• x has order 1/c^2
• γ has order 1/c^2
• 1/r has order 1/c^2

PNExpansionParameter(pnsystem)

Create a PNTerm object representing the post-Newtonian expansion parameter $c$. This can be used to automatically create more complicated PNTerms, which combine to form a PNExpansion. This is a simple but effective way to write PN formulas while automatically tracking the PN order of each term.

PNExpansion{N,T,NMax}

This object can be multiplied by a scalar or another PNExpansion object, and contains a tuple of coefficients. The coefficients are stored in the order of the expansion, with the zeroth-order coefficient first. Multiplication by a scalar just multiplies each of the elements. Multiplication by another PNExpansion object is done by convolving the two tuples of coefficients.

Blanchet (2014) defines the post-Newtonian expansion parameter as follows:

This parameter represents essentially a slow motion estimate $ϵ ∼ 𝑣/𝑐$, where $𝑣$ denotes a typical internal velocity. By a slight abuse of notation, following Chandrasekhar et al. [...], we shall henceforth write formally $ϵ ≡ 1/𝑐$, even though $ϵ$ is dimensionless whereas $𝑐$ has the dimension of a velocity. Thus, $1/𝑐 ≪ 1$ in the case of post-Newtonian sources. The small post-Newtonian remainders will be denoted $𝒪(1/𝑐^𝑛)$. Furthermore, [...] we shall refer to a small post-Newtonian term with formal order $𝒪(1/𝑐^𝑛)$ relative to the Newtonian acceleration in the equations of motion, as $\frac{𝑛}{2}\text{PN}$.

Therefore, we consider the coefficients of the PNExpansion to be a polynomial in $1/𝑐$. Here, the type parameter N corresponds to the number of elements actually present in the tuple of coefficients, and T is the type of the coefficients. The NMax parameter is the maximum number of elements, related to the usual PN order by

\[\text{pn_order} = \frac{\texttt{NMax}-1} {2}.\]

The N parameter is not related to the PN order; it is just used by Julia to know how many elements are currently in the coefficients, but is required to be 1 ≤ N ≤ NMax.

@pn_expansion [pnsystem] expansion

Gather terms in expansion by the powers of $1/c$ involved, zeroing out any terms with powers of $1/c$ higher than the pnsystem's PNOrder parameter, and combine the terms using the PNExpansionReducer specified in argument of the function that includes this macro call.

This expansion is achieved by setting — inside a let block created by this macro —

Note that the pnsystem argument can be inserted automatically by @pn_expression.

@pn_expression [arg_index=1] func

This macro takes the function func, looks for various symbols inside that function, and if present defines them appropriately inside that function.

The first argument to this macro is arg_index, which just tells us which argument to the function func is a PNSystem. For example, the variables defined in PostNewtonian.FundamentalVariables all take a single argument of pnsystem, which is used to compute the values for those variables; this macro just needs to know where to find pnsystem.

Once it has this information, there are five types of transformations it will make:

1. Adds a keyword argument pn_expansion_reducer::Val{PNExpansionReducer}=Val(sum) to the function signature. This is used to determine how to reduce the PN expansion terms. The default is Val(sum), which will just return a single number, but Val(identity) can be used to return the expansion. This should be used inside the function as PNExpansionReducer, and will be automatically used inside any @pn_expansion.
2. For every fundamental or derived variable, the name of that variable used in the body of func will be replaced by its value when called with pnsystem. For example, you can simply use the symbols M₁ or μ in your code, rather than calling them as M₁(pnsystem) or μ(pnsystem) every time they appear.
3. Every Irrational defined in Base.MathConstants or PostNewtonian.MathConstants will be transformed to the eltype of pnsystem. This lets you naturally use such constants in expressions like 2π/3 without automatically converting to Float64.
4. Each of a short list of functions given by unary_funcs in utilities/macros.jl will first convert their arguments to the eltype of pnsystem. In particular, you can use expressions like √10 or ln(2) without the result being converted to a Float64.
5. Insert the pnsystem argument as the first argument to each occurrence of @pn_expansion that needs it.

To be more explicit, the first three are achieved by defining the relevant quantities in a let block placed around the body of func, so that the values may be used efficiently without recomputation.

If you need to use one of the fundamental- or derived-variable functions as arguments of values other than those encapsulated in pnsystem, you'll need to address them explicitly with the module name — as in PostNewtonian.v(;Ω, M).

type_converter(pnsystem, x)

Convert x to a type appropriate for the float type of pnsystem.

This is sort of an expansion of the convert function, but with nicer syntax for types from this package, including the ability to do really weird things for SymbolicPNSystem. It's needed to ensure that the types of variables and constants are correct when we use them in expressions, rather than just assuming everything is a Float64.

decreasing_v_terminator([quiet])

Construct termination criterion to stop integration when v is decreasing.

Note that some systems may truly have decreasing v as physical solutions — including eccentric systems and possibly precessing systems. You may prefer to implement another solution, like detecting when v decreases below some threshold, or detecting when v is decreasing too quickly. See this function's source code for a simple

If this terminator is triggered while v is less than 0.35, a warning will always be issued; otherwise an info message will be issued only if the quiet flag is set to false.

dtmin_terminator(T, [quiet])

Construct termination criterion to terminate when dt drops below √eps(T).

Pass force_dtmin=true to solve when using this callback. Otherwise, the time-step size may decrease too much within a single time step, so that the integrator itself will quit before reaching this callback, leading to a less graceful exit.

If this terminator is triggered while v is less than 0.35, a warning will always be issued; otherwise an info message will be issued only if the quiet flag is set to false.

nonfinite_terminator()

Construct termination criterion to stop integration when any NaN or Inf is found in the data after an integration step.

If this terminator is triggered, a warning will always be issued.

termination_backwards(v₁, [quiet])

Construct termination criteria of solving PN evolution backwards in time

These criteria include checking that the masses are positive and the dimensionless spins are less than 1, as well as ensuring that the evolution will terminate at v₁.

The optional quiet argument will silence informational messages about reaching the target value of v₁ if set to true, but warnings will still be issued when terminating for other reasons.

termination_forwards(vₑ, [quiet])

Construct termination criteria of solving PN evolution forwards in time

These criteria include checking that the masses are positive and the dimensionless spins are less than 1, as well as ensuring that the evolution will terminate at vₑ.

The optional quiet argument will silence informational messages about reaching the target value of vₑ if set to true, but warnings will still be issued when terminating for other reasons.

combine_solutions(sol₋, sol₊)

Combine ODESolutions

This function is internal to this package. It is not entirely general, but allows us to combine the backwards- and forwards-in-time solutions of the PN orbital-evolution ODE equations into a single ODESolution object that should behave just as if it were the result of solve. In particular, indexing, interpolation, and iterations should behave exactly as described in the DifferentialEquations docs.

γₑ

Euler's constant (also known as the Euler–Mascheroni constant) is defined as the limit as $n \to \infty$ of the difference between the $n$th partial sum of the harmonic series and $\log(n)$. This is OEIS sequence A001620.

julia> PostNewtonian.γₑ
γₑ = 0.5772156649015...

julia> n=10_000_000; sum(1 ./ (1:n))-log(n)
0.5772157149015307

ζ3
apery

Apéry's constant is defined as $ζ(3)$, where $ζ$ is the Riemann zeta function. This is OEIS sequence A002117.

julia> PostNewtonian.apery
ζ3 = 1.2020569031595...

julia> PostNewtonian.ζ3
ζ3 = 1.2020569031595...

julia> sum((1:10_000_000).^-3)
1.2020569031595896

ln2
log2

The natural logarithm of 2. This is OEIS sequence A002162.

julia> PostNewtonian.ln2
ln2 = 0.6931471805599...

julia> exp(PostNewtonian.ln2)
2.0

julia> exp(big(PostNewtonian.ln2))
2.0

ln3
log3

The natural logarithm of 3. This is OEIS sequence A002391.

julia> PostNewtonian.ln3
ln3 = 1.0986122886681...

julia> exp(PostNewtonian.ln3)
3.0000000000000004

julia> exp(big(PostNewtonian.ln3))
3.0

ln5
log5

The natural logarithm of 5. This is OEIS sequence A016628.

julia> PostNewtonian.ln5
ln5 = 1.6094379124341...

julia> exp(PostNewtonian.ln5)
4.999999999999999

julia> exp(big(PostNewtonian.ln5))
5.0

ln³╱₂
log³╱₂
log3halves

The natural logarithm of 3//2. This is OEIS sequence A016578.

julia> PostNewtonian.ln³╱₂
ln³╱₂ = 0.4054651081081...

julia> exp(PostNewtonian.ln³╱₂)
1.5

julia> exp(big(PostNewtonian.ln³╱₂))
1.5

ln⁵╱₂
log⁵╱₂
log5halves

The natural logarithm of 5//2. This is OEIS sequence A016579.

julia> PostNewtonian.ln⁵╱₂
ln⁵╱₂ = 0.9162907318741...

julia> exp(PostNewtonian.ln⁵╱₂)
2.5

julia> exp(big(PostNewtonian.ln⁵╱₂))
2.5

truncated_series_inverse(a)
truncated_series_inverse!(b, a)

Given the coefficients a of a series, find the coefficients b of the multiplicative inverse of that series, up to the order of the original series. That is, if

\[A \colonequals \sum_{i=0}^n a_{i-1} v^i,\]

then we return the coefficients b of the series

\[B \colonequals \sum_{i=0}^n b_{i-1} v^i\]

such that

\[A\, B = 1 + \mathcal{O}(v^{n+1}).\]

See lagrange_inversion for the compositional inverse (a.k.a. reversion), which returns the coefficients of $f^{-1}$ such that $f^{-1}(f(v)) = v + \mathcal{O}(v^{n+1})$.

Note that this function returns the coefficients of the inverse, rather than its value. This is relevant for use in truncated_series_product and truncated_series_ratio — the latter of which just combines the former with this function.

For example, suppose the input coefficients represent the series

\[A \colonequals \sum_{i=0}^n a_{i-1} v^i.\]

(Remember that Julia's indexing is 1-based, so we subtract 1 to get the index of the coefficient of $v^i$.) Then we return the coefficients b of the series

\[B \colonequals \sum_{i=0}^n b_{i-1} v^i\]

such that

\[A\, B = 1 + \mathcal{O}(v^{n+1}).\]

Explanation

The inverse coefficients can be computed fairly easily by induction. Start by defining

\[b_0 = 1/a_0.\]

Now, assuming that we've computed all coefficients up to and including $b_{i}$, we can compute $b_{i+1}$ from the condition that the term proportional to $v^{i+1}$ in the product of the series and its inverse must be zero. This gives

\[b_{i+1} = -b_0\sum_{j=1}^{i} a_j b_{i-j}.\]

truncated_series_product(a, b, v)

Evaluate the truncated product of the series a and b, which are expanded in powers of v.

Note that this function returns the value of the summation, rather than its coefficients.

Here we define the series in terms of the coefficients a and b as

\[A \colonequals \sum_{i=0}^n a_{i-1} v^i \qquad B \colonequals \sum_{i=0}^n b_{i-1} v^i,\]

and return the value of the product $A\, B$ truncated at $v^n$.

Internally, the sums are performed using evalpoly.

See also truncated_series_ratio.

truncated_series_ratio(a, b, v)

Evaluate the truncated ratio of the series a and b, which are expanded in powers of v.

Note that this function returns the value of the summation, rather than its coefficients.

Here we define the series in terms of the coefficients a and b as

\[A \colonequals \sum_{i=0}^n a_{i-1} v^i \qquad B \colonequals \sum_{i=0}^n b_{i-1} v^i,\]

and return the value of the ratio $A / B$ truncated at $v^n$.

This function simply combines truncated_series_product and truncated_series_inverse.

truncated_series_ratio(a, b)

Evaluate the truncated ratio of the series a and b, evaluated at expansion value 1. This is relevant when the expansion is not in the dynamic variable v, for example, but in powers of $1/c$ as in post-Newtonian expansions. (That is, when the v dependence is already include in the input coefficients.)

lagrange_inversion(a)

Compute the compositional inverse (a.k.a. reversion) of the power series expansion

\[f(x) = a[1]*x + a[2]*x^2 + \ldots + a[n-1]*x^{n-1}\]

about 0, where a is an NTuple. Note, in particular, that there is no constant term. The result is a similar NTuple b allowing us to write

\[f^{-1}(y) = b[1]*y + b[2]*y^2 + \ldots + b[n-1]*y^{n-1},\]

such that $f^{-1}(f(x)) = f(f^{-1}(y)) = x$ mod $x^n$.

See truncated_series_inverse for the multiplicative inverse.

When the constant coefficient $a_0$ is nonzero, the result must be expanded about a different point, which is done by evaluating the output as $f^{-1}(y-a_0)$. Similarly, if the original expansion is about a point $x_0 ≠ 0$, the result must be shifted by adding $x_0$ to the output.

Johansson (2015) summarizes this basic form of the algorithm nicely:

Our setting is the ring of truncated power series $R[[x]]/\langle x^n \rangle$ over a commutative coefficient ring $R$ in which the integers $1,...,n−1$ are cancellable (i.e., nonzero and not zero divisors). [...] If $f(x) = x/h(x)$ where $h(0)$ is a unit, then the compositional inverse or reversion $f^{-1}(x)$ satisfying $f(f^{-1}(x)) = f^{-1}(f(x)) = x$ exists and its coefficients are given by

\[[x^k] f^{-1}(x) = \frac{1}{k} [x^{k-1}] h(x)^k.\]

Note that Johansson also presents a pair of asymptotically faster algorithms for computing the compositional inverse. Because of the low orders of the power series expansions we typically work with, it is not clear if the improved scaling of those algorithms would actually be beneficial in practice, so we stick with the basic algorithm here — though they would not be too difficult to implement if needed.

x╱f_mod_xⁿ⁻¹(a)

Compute the truncated power series expansion of $x/f$ mod $x^{n-1}$, where $f$ is the power series expansion of a function $f(x)$

\[f(x) = a[1] x + a[2] x^2 + a[3] x^3 + \ldots + a[n-1] x^{n-1}.\]

Note, in particular that there is no constant term. The result is a power series expansion

\[h(x) = h[1] + h[2] x + h[3] x^2 + \ldots + h[n-1] x^{n-2},\]

which notably does have a constant term.

This function is essentially a helper function for the lagrange_inversion function.

hⁱ✖h_mod_xⁿ⁻¹(h, i)

Compute the truncated power series expansion of $h^i$ mod $x^{n-1}$, where $h$ is the power series expansion of a function $h(x)$

\[h(x) = h[1] + h[2] x + h[3] x^2 + \ldots + h[n-1] x^{n-2}.\]

The result is a power series expansion

\[h^i(x) = h^i[1] + h^i[2] x + h^i[3] x^2 + \ldots + h^i[n-1] x^{n-2}.\]

This function is essentially a helper function for the lagrange_inversion function.

apply_to_first_add!(expr, func)

Apply func to the first sub-expression found in a "prewalk"-traversal of expr that satisfies isadd. If func acts in place, so does this function. In either case, the expression should be returned.

find_symbols_of_type(mod, T)

Given a module mod (not just its name, but the actual imported module), find all objects inside that module that are instances of the given type T. The returned quantity is a vector of Symbols.

flatten_binary!(expr, symbols)

Flatten nested binary operations — that is, apply associativity repeatedly.

isadd(x)

Return true if the Expr x is a call to (+) or :+.

iscall(x, symbols)

Return true if the Expr x is a call to any element of symbols.

ismul(x)

Return true if the Expr x is a call to (*) or :*.

value(x)

Return x or the value wrapped by the Dual number x