H!(H, expiβ, ℓₘₐₓ, m′ₘₐₓ, H_rec_coeffs)
H!(H, expiβ, ℓₘₐₓ, m′ₘₐₓ, H_rec_coeffs, Hindex)

Compute the $H$ matrix defined by Gumerov and Duraiswami [8].

This computation forms the basis for computing Wigner's $d$ and $𝔇$ matrices via d_matrices! and D_matrices!, the spin-weighted spherical harmonics via sYlm_values!, and for transforming from values of spin-weighted spherical functions evaluated on a grid to the corresponding mode weights via map2salm.

Due to symmetries, we only need to compute ~1/4 of the elements of this matrix, so only those elements with $m ≥ |m′|$ are computed. The relevant indices of the H vector are computed based on the Hindex function — which defaults to WignerHindex, but could reasonably be WignerDindex if the input H vector contains all valid indices. However, it is assumed that the storage scheme used for H is such that the successive $m$ values are located in successive elements.

If $m′ₘₐₓ < ℓₘₐₓ$, we don't even need 1/4 of the elements, and only values with $|m′| ≤ m′ₘₐₓ$ will be computed. This is particularly useful for computing spin-weighted spherical harmonics.

Note that the recursion coefficients H_rec_coeffs should be the quantity returned by H_recursion_coefficients.

H_recursion_coefficients(ℓₘₐₓ, T)

Pre-compute constants used in Wigner H recursion.

ALFcompute(expiβ, nmax)
ALFcompute!(P̄, expiβ, nmax)
ALFcompute(expiβ, nmax, recursion_coefficients)
ALFcompute!(P̄, expiβ, nmax, recursion_coefficients)

Compute the "fully normalized" Associated Legendre Functions up to some maximum n value nmax.

These functions can take a vector P̄, to store the data, stored in order of increasing m most rapidly varying and then increasing n. If not supplied, P̄ will be constructed for you and returned.

The optional recursion_coefficients argument must be an ALFRecursionCoefficients, which stores various constant coefficients used in the recursion. This object requires more than 3x the memory and more than 20x the time to compute a single P̄ vector without this argument, but passing it will typically speed up the calculation of each P̄ by a factor of 8x or so. Thus, if you expect to compute P̄ more than a few times, it will take less time to pre-compute those factors, and pass them to this function.

Note that the base real type will be inferred from the (complex) type of expiβ. If present, the base types of P̄ and recursion_coefficients must agree.

λ_recursion_initialize(cosθ, sin½θ, cos½θ, s, ℓ, m)

This provides initial values for the recursion to find ${}_{s}\lambda_{\ell,m}$ along indices of increasing $\ell$, due to Kostelec & Rockmore Specifically, this function computes values with $\ell=m$.

\[{}_{s}\lambda_{\ell,m}(\theta) := {}_{s}Y_{\ell,m}(\theta, 0) = (-1)^m\, \sqrt{\frac{2\ell+1}{4\pi}} d^\ell_{-m,s}(\theta)\]

λ_iterator(θ, s, m)

Construct an object to iterate over ₛλₗₘ values.

The $ₛλₗₘ(θ)$ function is defined as the spin-weighted spherical harmonic evaluated at spherical coordinates $(θ, ϕ)$, with $ϕ=0$. In particular, note that it is real-valued. The return type is determined by the type of θ (or more precisely, cos½θ).

This algorithm by Kostelec & Rockmore allows us to iterate over increasing $ℓ$ values, for given fixed $s$ and $m$ values.

Note that this iteration has no inherent bound, so if you try to iterate over all values, you will end up in an infinite loop. Instead, you can zip this iterator with another:

θ = 0.1
s = -2
m = 1
λ = λ_iterator(θ, s, m)
Δ = max(abs(s), abs(m))
for (ℓ, ₛλₗₘ) ∈ zip(Δ:Δ+5, λ)
    @show (ℓ, ₛλₗₘ)
end

Alternatively, you could use Iterates.take(λ, 6), for example.

Note that the iteration always begins with ℓ = Δ = max(abs(s), abs(m)).

Simple iterator to count down to 0, with alternating signs

julia> collect(AlternatingCountdown(5))
11-element Vector{Int64}:
  5
 -5
  4
 -4
  3
 -3
  2
 -2
  1
 -1
  0

ₛ𝐘(s, ℓₘₐₓ, [T=Float64], [Rθϕ=golden_ratio_spiral_rotors(s, ℓₘₐₓ, T)])

Construct a matrix of $ₛYₗₘ(Rθϕ)$ values for the input s and all nontrivial $(\ell, m)$ up to ℓₘₐₓ.

This is a fast and accurate method for mapping between the vector of spin-weighted spherical-harmonic mode weights $ₛ𝐟ₗₘ$ and the vector of function values on the sphere $ₛ𝐟ⱼₖ$, as

\[ₛ𝐟ⱼₖ = ₛ𝐘\, ₛ𝐟ₗₘ,\]

where the right-hand side represents the matrix-vector product. As usual, we assume that the $ₛ𝐟ₗₘ$ modes are ordered by increasing $m ∈ [-ℓ:ℓ]$, and $ℓ ∈ [|s|:ℓₘₐₓ]$. The ordering of the $ₛ𝐟ⱼₖ$ values will be determined by the ordering of the argument Rθϕ.

Note that the number of modes need not be the same as the number of points on which the function is evaluated, which would imply that the output matrix is not square. To be able to invert the relationship, however, we need the number of points $ₛ𝐟ⱼₖ$ to be at least as large as the number of modes $ₛ𝐟ₗₘ$.

Note that the usefulness of this approach is limited by the fact that the size of this matrix scales as ℓₘₐₓ⁴. As such, it is mostly useful only for ℓₘₐₓ of order dozens, rather than — say — the tens of thousands that CMB astronomy or lensing require, for example.

Direct application and inversion of this matrix are used in the "direct" methods of $s$-SHT transformations. See SSHTDirect for details about the implementation.

map2salm!(salm, map, spin, ℓmax)
map2salm!(salm, map, plan)

Transform map values sampled on the sphere to ${}_sa_{\ell, m}$ modes in place.

For details, see map2salm.

map2salm(map, spin, ℓmax)
map2salm(map, plan)

Transform map values sampled on the sphere to ${}_sa_{\ell, m}$ modes.

The map array should have size Nφ along its first dimension and Nϑ along its second; any number of dimensions may follow. The spin must be entered explicitly, and ℓmax is the highest ℓ value you want in the output.

For repeated applications of this function with different values of map, it is more efficient to pre-compute plan using plan_map2salm. These functions will create a new salm array on each call. To operate in place on a pre-allocated salm array, use map2salm!.

The core of this function follows the method described by Reinecke and Seljebotn.

plan_map2salm(map, spin, ℓmax)

Precompute values to use in executing map2salm or map2salm!.

The arguments to this function exactly mirror those of the first form of map2salm, and all but the first argument in the first form of map2salm!. The plan returned by this function can be passed to the second forms of those functions to avoid some computation and allocation costs.

Note that the plan object is not thread safe; a separate plan should be created for each thread that will use one, or locks should be used to ensure that a single plan is not used at the same time on different threads.