$s$-SHT Transformations

Supertype of storage for spin-spherical-harmonic transforms

SSHT(s, ℓₘₐₓ; [method="Direct"], [T=Float64], [kwargs...])

Construct an SSHT object to transform between spin-weighted spherical-harmonic mode weights and function values — performing an $s$-SHT.

This object behaves similarly to an AbstractFFTs.Plan object — specifically in the ability to use the semantics of algebra to perform transforms. For example, if the function values are stored as a vector f, the mode weights as f̃, and the SSHT as 𝒯, then we can compute the function values from the mode weights as

f = 𝒯 * f̃

or solve for the mode weights from the function values as

f̃ = 𝒯 \ f

The first dimensions of f̃ must index the mode weights (as usual, for ℓ∈abs(s):ℓₘₐₓ and m∈-ℓ:ℓ) and the first index of f must index the locations at which the function is evaluated. Any following dimensions will be broadcast over. Note that certain types will broadcast using Julia threads, while others will broadcast using BLAS threads. The relevant number of threads must be set appropriately.

Certain SSHT types (currently, only Minimal and Direct) also have an option to always act in place — meaning that they simply re-use the input storage, even when used in an expression like 𝒯 \ f. The option must be passed as the inplace argument to the constructors, and part of the type of the resulting object. Regardless of the value of that option, for those types where the option exists, it is also possible to use mul! and ldiv! from the LinearAlgebra package to force operation in place.

Note that different algorithms require different "pixelizations", or sets of Rotors on which to evaluate the function. These can be obtained from the SSHT object using the pixels and rotors functions.

pixels(𝒯)

Return the spherical coordinates (θ, ϕ) on which the spin-weighted spherical harmonics are evaluated. See also rotors, which provides the actual Rotors on which they are evaluated.

rotors(𝒯)

Return the Rotors on which the spin-weighted spherical harmonics are evaluated. See also pixels, which provides the corresponding spherical coordinates.

SSHTDirect(s, ℓₘₐₓ; decomposition=LinearAlgebra.qr, T=Float64, Rθϕ=golden_ratio_spiral_rotors(s, ℓₘₐₓ, T), inplace=true)

Construct an $s$-SHT object that uses the "Direct" method; see ₛ𝐘 for details about the method and optional arguments. Also see SSHT for general information about how to use these objects.

By default, this uses precisely optimal sampling — meaning that the number of points on which the function is evaluated, represented by Rθϕ, is equal to the number of modes. However, it is equally possible to evaluate on more points than there are modes. This can be useful, for example, when processing multiple fields with different spin weights; the function could be evaluated on points appropriate for the lowest value of $|s|$, and therefore could also be used to solve for fields of all other spin weights.

Note that in-place operation is possible for this type when the length of the input Rθϕ is equal to the number of modes given s and ℓₘₐₓ — and is the default behavior when possible. See SSHT for description of in-place operation.

This method is typically better than other current implementations for $ℓₘₐₓ ≲ 24$, both in terms of speed and accuracy. However, this advantage quickly falls away. A warning will be issued if ℓₘₐₓ is greater than about 64, because this method is not likely to be the most efficient or most accurate choice.

Storage for Minimal spin-spherical-harmonic transform

The Minimal algorithm was described by Elahi et al., and allows for the minimal number of function samples.

SSHTMinimal(s, ℓₘₐₓ; kwargs...)

Construct a SSHTMinimal object directly. This may also be achieved by calling the main SSHT function with the same keywords, along with method="Minimal".

This object uses the algorithm described by Elahi et al.

The basic floating-point number type may be adjusted with the keyword argument T, which defaults to Float64.

The SSHs are evaluated on a series of "rings" at constant colatitude. Their locations are specified by the θ keyword argument, which defaults to sorted_rings(s, ℓₘₐₓ, T). The first element of θ is the colatitude of the smallest ring (containing $2s+1$ elements), and so on to the last element of θ, which is the colatitude of the largest ring (containing $2ℓₘₐₓ+1$ elements).

Whenever T is either Float64 or Float32, the keyword arguments plan_fft_flags and plan_fft_timelimit may also be useful for obtaining more efficient FFTs. They default to FFTW.ESTIMATE and Inf, respectively. They are passed to AbstractFFTs.plan_fft.

Note that, because this algorithm achieves optimal dimensionality, the transformation will be performed in place by default. If this is not desired, pass the keyword argument inplace=false. This will cause the algorithm to copy the input and perform in-place transformation on that copy.

Storage for spin-spherical-harmonic transform

The algorithm was described in by Reinecke and Seljebotn.

SSHTRS(s, ℓₘₐₓ; kwargs...)

Construct a SSHTRS object directly. This may also be achieved by calling the main SSHT function with the same keywords, along with method="RS".

This object uses the algorithm described by Reinecke and Seljebotn.

The basic floating-point number type may be adjusted with the keyword argument T, which defaults to Float64.

The SSHs are evaluated on a series of "rings" at constant colatitude. Their locations are specified by the θ keyword argument, which defaults to fejer1_rings(2ℓₘₐₓ+1, T). If this is changed, the user should also provide the corresponding quadrature_weights argument — the default being fejer1(length(θ), T).

On each of these rings, an FFT is performed. To reach the band limit of $m = ± ℓₘₐₓ$, the number of points along each ring must therefore be at least $2ℓₘₐₓ+1$, but may be greater. For example, if $2ℓₘₐₓ+1$ does not factorize neatly into a product of small primes, it may be preferable to use $2ℓₘₐₓ+2$ points along each ring. (In that case, whenever ℓₘₐₓ is 1 less than a power of 2, the number of points will be exactly a power of 2, which is usually particularly efficient.) The number of points on each ring can be modified independently, if given as a vector with the same length as θ, or as a single number which is assumed to be the same for all rings.

Whenever T is either Float64 or Float32, the keyword arguments plan_fft_flags and plan_fft_timelimit may also be useful for obtaining more efficient FFTs. They default to FFTW.ESTIMATE and Inf, respectively. They are passed to AbstractFFTs.plan_fft.

clenshaw_curtis_rings(N, [T=Float64])

Values of the colatitude coordinate ($θ$) appropriate for quadrature by the Clenshaw-Curtis rule, using weights provided by clenshaw_curtis.

Note that the first argument to this function is N, rather than the ℓₘₐₓ used in some other functions. For spin-weighted spherical harmonics, you may want to use N=2ℓₘₐₓ+1.

fejer1_rings(N, [T=Float64])

Values of the colatitude coordinate ($θ$) appropriate for quadrature by Fejér's first rule, using weights provided by fejer1.

Note that the first argument to this function is N, rather than the ℓₘₐₓ used in some other functions. For spin-weighted spherical harmonics, you may want to use N=2ℓₘₐₓ+1.

fejer2_rings(N, [T=Float64])

Values of the colatitude coordinate ($θ$) appropriate for quadrature by Fejér's second rule, using weights provided by fejer2.

Note that the first argument to this function is N, rather than the ℓₘₐₓ used in some other functions. For spin-weighted spherical harmonics, you may want to use N=2ℓₘₐₓ+1.

golden_ratio_spiral_pixels(s, ℓₘₐₓ, [T=Float64])

Cover the sphere 𝕊² with pixels generated by the golden-ratio spiral. Successive pixels are separated by the azimuthal angle $Δϕ = 2π(2-φ)$, and are uniformly distributed in $\cos θ$.

This is also known as the "Fibonacci sphere" or "Fibonacci lattice".

Visually, this is a very reasonable-looking pixelization, with fairly uniform distance between neighbors, and approximate isotropy. No two pixels will share the same values of either $θ$ or $ϕ$. Also note that no point is present on either the North or South poles.

The returned quantity is a vector of 2-SVectors providing the spherical coordinates of each pixel. See also golden_ratio_spiral_rotors for the corresponding Rotors.

golden_ratio_spiral_rotors(s, ℓₘₐₓ, [T=Float64])

Cover the sphere 𝕊² with pixels generated by the golden-ratio spiral.

See golden_ratio_spiral_pixels for more detailed explanation. The quantity returned by this function is a vector of Rotors providing each pixel.

sorted_ring_pixels(s, ℓₘₐₓ, [T=Float64])

Cover the sphere 𝕊² with $(ℓₘₐₓ+1)²-s²$ pixels distributed in rings provided by sorted_rings; see that function's documentation for more description.

The returned quantity is a vector of 2-SVectors containing the spherical coordinates of each pixel. See also sorted_ring_rotors for the corresponding Rotors.

sorted_ring_rotors(s, ℓₘₐₓ, [T=Float64])

Cover the sphere 𝕊² with $(ℓₘₐₓ+1)²-s²$ pixels distributed in rings provided by sorted_rings; see that function's documentation for more description.

The returned quantity is a vector of Rotors. See also sorted_ring_rotors for the corresponding spherical coordinates.

sorted_rings(s, ℓₘₐₓ, [T=Float64])

Compute locations of a series of rings labelled by $j ∈ |s|:ℓₘₐₓ$ (analogous to $ℓ$), where each ring will contain $k = 2j+1$ (analogous to $m$) pixels distributed evenly around the ring. These rings are then sorted, so that the ring with the most pixels ($j = ℓₘₐₓ$) is closest to the equator, and the next-largest ring is placed just above or below the equator (depending on the sign of $s$), the next just below or above, and so on. This is generally a fairly good first guess when minimizing the condition number of matrices used to solve for mode weights from function values. In particular, I use this to initialize the Minimal algorithm, which is then fed into an optimizer to fine-tune the positions of the rings.

This function does not provide the individual pixels; it just provides the colatitude values of the rings on which the pixels will be placed. The pixels themselves are provided by sorted_ring_pixels.

clenshaw_curtis(n, [T])

Compute n weights for the Clenshaw-Curtis rule, corresponding to n evenly spaced nodes from 0 to π inclusive. That is, the nodes are located at

\[\theta_k = k \frac{\pi}{n-1} \quad k=0, \ldots, n-1.\]

This function uses Waldvogel's method.

The type T may be any AbstractFloat, but defaults to Float64.

fejer1(n, [T])

Compute n weights for Fejér's first rule, corresponding to n evenly spaced nodes from 0 to π inclusive. That is, the nodes are located at

\[\theta_k = k \frac{\pi}{n-1} \quad k=0, \ldots, n-1.\]

This function uses Waldvogel's method.

The type T may be any AbstractFloat, but defaults to Float64.

fejer2(n, [T])

Compute n weights for Fejér's second rule, corresponding to n evenly spaced nodes between 0 and π exclusive. That is, the nodes are located at

\[\theta_k = k \frac{\pi}{n+1} \quad k=1, \ldots, n.\]

This function uses Waldvogel's method. However, contrary to Waldvogel's notation, this routine does not include the weight corresponding to the ϑ=0 or π nodes, which both have weight 0.

The type T may be any AbstractFloat, but defaults to Float64.