API Reference

CubicSplineCache{T<:AbstractFloat}

Precomputed LU factorization for natural cubic spline interpolation on a fixed time grid. Stores the interval widths, their reciprocals, and the Thomas factors for the tridiagonal second-derivative system.

After constructing the cache once from the knot vector t, repeated interpolation of new data vectors dₖ requires only an O(N) forward/backward sweep to compute the second derivatives d̈ₖ, plus O(1) per evaluation point.

Fields

• h: interval widths h[j] = t[j+1]−t[j], length N−1
• h⁻¹: reciprocals 1/h[j], length N−1
• l: Thomas sub-diagonal factors lᵢ = h[i+1]/uᵢ for i=1..N−3
• u⁻¹: reciprocals of the modified diagonal 1/uᵢ for i=1..N−2

CubicSplineCache(t)

Precompute the Thomas LU factorization of the natural-cubic-spline tridiagonal system for knot vector t. t must be strictly increasing with at least 4 elements.

Natural spline system

For N knots, the N−2 interior second derivatives d̈[2..N-1] satisfy the symmetric tridiagonal system

h[j-1]*d̈[j-1] + 2(h[j-1]+h[j])*d̈[j] + h[j]*d̈[j+1] = rₖ[j]

with d̈[1]=d̈[N]=0 (natural boundary conditions) and

rₖ[j] = 6*( (d[j+1]−d[j])/h[j] − (d[j]−d[j-1])/h[j-1] ).

Denoting the system Ad̈ₖ = rₖ, the Thomas LU factorization gives

A = LU,  L unit lower bidiagonal with factors lᵢ,
          U upper bidiagonal with diagonal uᵢ and super-diagonal sᵢ = h[i+1].

Thomas forward sweep

With the cache built, the forward sweep for a new data vector d is

z[2] = r[2] * u⁻¹[1]
z[j] = ( r[j] − h[j-1] * z[j-1] ) * u⁻¹[j-1]   for j = 3..N-1

and the backward sweep recovers d̈ from z:

d̈[N-1] = z[N-1]
d̈[j]   = z[j] − l[j-1] * d̈[j+1]                for j = N-2:-1:2

Evaluation

On interval [t[j], t[j+1]] at offset τ = t_query − t[j]:

c = (d̈[j+1]−d̈[j]) * h⁻¹[j] / 6
b = d̈[j] / 2
a = h⁻¹[j]*(d[j+1]−d[j]) − h[j]/6*(2d̈[j]+d̈[j+1])
S = d[j] + τ*(a + τ*(b + τ*c))

DataComponents{C, εᴵ}

Encodes a fixed set of waveform data components at the type level. C is an NTuple{N,Symbol} whose elements are drawn from ValidDataComponents:

(:ψ₀, :ψ₁, :ψ₂, :ψ₃, :ψ₄, :σ, :h, :News, :φ₀, :φ₁, :φ₂)

The parameter εᴵ is the sign of the time direction for which these components are defined: +1 for $ℐ⁺$ (outgoing) and -1 for $ℐ⁻$ (incoming). The default value is +1, since most applications will be for $ℐ⁺$.

The inputs may alternatively be strings; any reasonable spelling will be parsed into the canonical symbol form. For example, DataComponents("Psi_3", "psi4", "sigma") will be parsed into DataComponents(:ψ₃, :ψ₄, :σ). The constructor will throw an error if the input components are invalid.

Examples:

DataComponents(:ψ₄)                          # gravitational waves only
DataComponents(:ψ₄, :ψ₃, :ψ₂)                # top three Weyl components
DataComponents(:ψ₀, :ψ₁, :ψ₂, :ψ₃, :ψ₄, :σ)  # full Weyl set with strain
DataComponents(:φ₀, :φ₁, :φ₂)                # Faraday components
DataComponents(:φ₀, :φ₁, :φ₂; εᴵ=-1)         # Faraday components on ℐ⁻

aberration(RRₚᵢ, v⃗; emitted=true)

Transform from the rotor RRₚᵢ in the boosted frame to the corresponding rotor in the rest frame, given the boost velocity vector v⃗. This is derived in Appendix C of [5]. It is also described in Penrose-Rindler Vol. 1, around Eq. (1.3.5).

A rotor on the sphere encodes not only a point (direction) but also an orientation of the tangent frame at that point — the full bundle on which spin-weighted functions are defined. This function implements the complete transformation of that bundle under a boost: both the direction aberration and the rotation of the tangent frame (the spin-weight phase factor).

The input RRₚᵢ = R * Rₚᵢ is the product of the overall frame rotation R and the pixel rotor Rₚᵢ. The rotation is applied before the boost correction, matching the standard BMS transformation order (translate → rotate → boost). The boost correction B′ is then computed from the rotated direction n̂′ = RRₚᵢ(𝐤).

All rotors — input and output — are expressed in the same standard spatial coordinate basis; the "boosted frame" label refers to which physical direction the input rotor points to, not to any change of coordinates.

The emitted keyword

The emitted keyword selects which null infinity the field lives on:

• emitted = true (default): future null infinity ℐ⁺ — outgoing radiation emitted by the source. A wave propagating at angle Θ in the rest frame appears at the larger angle Θ̑ > Θ in the boosted frame.
• emitted = false: past null infinity ℐ⁻ — incoming radiation received by the observer. A wave arriving from direction n̂ in the rest frame appears shifted toward the boost axis in the boosted frame, so Θ̑ < Θ.

Penrose-Rindler work on the past celestial sphere (emitted = false), which reverses the sign compared to the emitted = true convention used throughout this package; hence the sign difference noted at Eq. (1.3.5) there.

Formula

The function returns B′ * RRₚᵢ, where B′ is the aberration rotor. B′ is computed as exp(n̂′ × v⃗ ⋅ δ/2) where δ = (Θ̑−Θ)/|n̂′ × v⃗| absorbs the normalization. Here Θ̑ is the polar angle of n̂′ relative to v⃗ in the boosted frame and Θ is the corresponding angle in the rest frame, related by

\[\tan(Θ/2) = e^{-εφ} \tan(Θ̑/2), \quad φ = \operatorname{atanh}(β), \quad ε = \begin{cases}+1 & \text{emitted} \\ -1 & \text{received}\end{cases}\]

Θ̑ is computed via atan(|n̂′ × v⃗|, v⃗ ⋅ n̂′) (a scale-invariant atan2 that requires no division by |v⃗|), keeping the computation well-conditioned.

Numerical stability

When β sinΘ̑ = |n̂′ × v⃗| is small — either because β → 0 or because n̂′ is nearly (anti-)parallel to v⃗ — the exact formula is ill-conditioned. In that regime we use

\[\exp\!\left(\frac{n̂′ \times v⃗}{|n̂′ × v⃗|}\,\frac{Θ̑ - Θ}{2}\right) = \cos\!\left(\frac{Θ̑ - Θ}{2}\right) + \frac{n̂′ \times v⃗}{|n̂′ × v⃗|}\,\sin\!\left(\frac{Θ̑ - Θ}{2}\right),\]

and expand the cos term and the ratio $\sin((Θ̑-Θ)/2)\,/\,(β \sin Θ̑)$ as Taylor series in εβ to fifth order. Because the error in a degree-5 truncation is $O((β \sin Θ̑)^6)$, the Taylor branch is used when $β \sin Θ̑ < \epsilon^{1/3}$; the resulting error is $O(\epsilon^2)$, well below floating-point precision.

The implementation is compatible with arbitrary floating-point types, including dual numbers for automatic differentiation. Branch selection is based on value(β sinΘ̑) so that the derivative of the chosen branch is always evaluated, with no discontinuity in derivatives at the branch boundary.

component_index(dc::DataComponents{C}, ::Val{S})

Return the index of component S within dc, or nothing if absent. The return type is inferred as a compile-time constant when dc has a concrete type.

compute_t′(t, αₚ, Rₚ, v⃗)

Compute the new time samples t′ corresponding to the input time samples t after a BMS transformation with supertranslation αₚ and boost velocity v⃗. The Rₚ describe the locations of the pixels.

The objective is to create a new time grid that has the same number of samples as t and has roughly the same spacing, while accounting for the fact that some parts of the cylinder (the subset of ℐ⁺) on which we have data will need to be dropped because the new slices at the ends will be "tilted", and thus won't have a complete sphere of data on which to compute modes.

The result is a simple rescaling:

scale = (t′ₘₐₓ - t′ₘᵢₙ) / (tₘₐₓ - tₘᵢₙ)
t′ = @. t′ₘᵢₙ + scale * (t - tₘᵢₙ)

For convenience, this function also returns tᵪ, the crossover time at which t′=t (the fixed point of the t ↦ t′ map), which we can derive from the above formula as

tᵪ = (t′ₘᵢₙ - scale * tₘᵢₙ) / (1 - scale)

conformal_weight(::Val{S})

Return the conformal weight (spin weight + boost weight) of field component S, which is the power of the conformal factor $k$ in the BMS transformation law.

Note that we are assuming that these fields represent the asymptotic values of the physical fields at null infinity, so they have already been rescaled by the appropriate power of the conformal factor to be finite and nonzero at null infinity. The transformation law for the physical fields (finite-radius Weyl and Faraday spinors) would have a different conformal weight, but the asymptotic fields are the ones we are transforming.

More specifically, the asymptotic Weyl spinor $ψ$ and Faraday spinor $φ$ are related to the finite-radius Weyl spinor $Ψ$ by $ψ = ωΨ$ and the finite-radius Faraday spinor $Φ$ by $φ = ωΦ$. The factor $ω$ is the conformal factor that goes to zero (but has nonzero derivative) at null infinity, which transforms as $ω′=kω$, so we pick up a factor of $k⁻¹$ in the transformation laws for $ψ$ and $φ$ compared to $Ψ$ and $Φ$. Since $Ψ$ and $Φ$ are the physical quantities, they do not change under coordinate transformations.

Meanwhile the basis spinors each transform with a factor of $1/√k$. The Weyl components $ψₙ$ are defined by contracting the Weyl spinor with four basis spinors, so they pick up a factor of $k⁻²$ from the basis spinors and an additional factor of $k⁻¹$ from the conformal factor, for a total of $k⁻³$. Similarly, the Faraday components $φₙ$ are defined by contracting the Faraday spinor with two basis spinors, so they pick up a factor of $k⁻¹$ from the conformal factor and an additional factor of $k⁻¹$ from the basis spinors, for a total of $k⁻²$.

diagnostics(data, data_components)

Compute power monitors for the input data.

has_component(dc::DataComponents{C}, ::Val{S})

Return true if component S is present in dc.

impose_reality(αᵢₙ, ℓₘₐₓ, εᵅ)

Given a set of mode weights αᵢₙ of a spin-weight 0 field, return a new set of mode weights that satisfy the reality condition $α_{ℓ,-m} = (-1)^m ᾱ_{ℓ,m}$. Simultaneously, pad the output array with zeros up to ℓₘₐₓ. The input αᵢₙ is expected to be ordered by increasing ℓ, starting from 0, then by increasing m within each ℓ. The output array is ordered in the same way, and has length (ℓₘₐₓ + 1)^2.

mix_components!(dataᵢⱼ, k⁻¹, ðt′╱k, ð²α, dc)

Apply the BMS component-mixing transformation to dataᵢⱼ. k⁻¹ is the inverse conformal factor for this pixel, ðt′╱k is the eth-derivative of the retarded time in the new frame divided by $k$, and ð²α is the sign-adjusted second anti-eth-derivative of the supertranslation (used for the strain/shear component).

Note that Julia specializes on the concrete type of dc. This means that the indexes into dataᵢⱼ for the various components are known at compile time, and the branches for which components are present or absent will be resolved at compile time. The result is that — even though the function body looks unwieldy and slow — it compiles down to minimal code with no branches and only the necessary components, making it very fast in practice. The @inline annotation encourages this specialization and inlining, especially when just a few components are being processed.

ncomponents(dc::DataComponents)

Return the number of components in dc.

spin_weight(::Val{S})

Return the spin weight of field component S.

spline_eval(dⱼ, dⱼ₊₁, d̈ⱼ, d̈ⱼ₊₁, hⱼ, h⁻¹ⱼ, τ)

Evaluate the natural cubic spline on interval j at offset τ = t_query − t[j], given the knot values dⱼ, dⱼ₊₁, the second derivatives d̈ⱼ, d̈ⱼ₊₁, the interval width hⱼ, and its reciprocal h⁻¹ⱼ. Horner form.

transform!(data, t, v⃗, R, αᵢₙ, dc, εᵅ=+1)

Transform the mode weights in data — sampled at times t — from the rest frame to the BMS-transformed frame. The BMS transformation is specified by the boost velocity v⃗, the overall rotation R, the supertranslation αᵢₙ, and the DataComponents descriptor dc. The transformation is performed in-place, modifying the input data array.

The complex data array is expected to have dimensions (Nᵗ, Nᵐ, Nᵈ), where Nᵗ is the number of time samples, Nᵐ is the number of modes, and Nᵈ is the number of data components (e.g., strain and/or Newman-Penrose Weyl components). The modes are expected to be ordered by increasing ℓ, then by increasing m within each ℓ. For data with spin weight $s ≠ 0$, the modes with $ℓ < |s|$ are expected to be present, but will be ignored. The maximum ℓ value is determined by the size of the second dimension of data as ℓₘₐₓ = √Nᵐ - 1. The array must have complex type, with the underlying real type being at least as wide as the types of the other inputs.

The t array is expected to have length Nᵗ, matching the first dimension of data. The v⃗ and R inputs are expected to be of types QuatVec and Rotor.

The αᵢₙ array must be a complex vector, and must be ordered as described above for the final dimension of data, though it may have a smaller ℓₘₐₓ. (A larger ℓₘₐₓ cannot be allowed, because the result would have higher angular dependence than data, which is not possible since we are transforming data in-place.) This array is expected to represent a supertranslation in the rest frame, which is a real-valued function (with spin weight 0) on the sphere. The reality condition is that the mode weights satisfy $α_{ℓ,-m} = (-1)^m ᾱ_{ℓ,m}$. This function will automatically impose this condition by averaging each mode with its complex-conjugate partner. This is done in a copy of the array for simplicity, rather than being done in place.

The dc argument is a DataComponents value specifying which field components are stored in data, in the order they appear along its third dimension. Because of the hierarchical nature of the BMS transformation, any Weyl component $ψᵢ$ must be accompanied by all higher-index components $ψⱼ$ for $j > i$. Note that DataComponents includes a sign indicating whether data represents data on $ℐ⁺$ if εᴵ = +1 or $ℐ⁻$ if εᴵ = -1.

The optional keyword εᵅ represents the sign in the time-transformation law $t′ = t - εᵅ α$.

transform!(data, t, v⃗, R, αᵢₙ; data_components=nothing, εᵅ=+1, εᴵ=+1)

Backward-compatible keyword-argument form. See the main docstring for details.

The data_components argument may be a DataComponents value, a tuple of symbols such as (:ψ₄, :ψ₃), or a sequence of strings that indicate those symbols. The strings are parsed in a flexible way, so that, for example, "psi4", "Psi_4", and "PSI₄" all indicate the same component :ψ₄. Alternatively, if the argument is nothing (the default), the first Nᵈ of (:σ, :ψ₄, :ψ₃, :ψ₂, :ψ₁, :ψ₀) will be chosen — though a warning will be issued.