$s$-SHT Transformations

SSHT{T}

Supertype of the spin-weighted spherical-harmonic transforms. See SSHT (the constructor), SSHTMatrix, SSHTRS and SSHTMinimal. T is the real type the transform works in.

SSHT(s, ℓₘₐₓ, [T=Float64]; [method="RS"], [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 (synthesis) as

f = 𝒯 * f̃

or solve for the mode weights from the function values (analysis) as

f̃ = 𝒯 \ f

The first dimension of f̃ must index the mode weights in the canonical ordering [f̃(ℓ, m) for ℓ ∈ abs(s):ℓₘₐₓ for m ∈ -ℓ:ℓ] (see Yindex), and the first dimension of f must index the locations at which the function is evaluated, in the order given by rotors or pixels. Any following dimensions will be broadcast over. The transform works in the floating-point type T, and its results are complex numbers of that type; the data given to it may be real or complex.

The available methods, each given as a string or a Symbol, are

• "RS" (default): SSHTRS, the ring-based algorithm of Reinecke and Seljebotn, which scales as $ℓₘₐₓ^3$ and is the choice for large $ℓₘₐₓ$;
• "Matrix": SSHTMatrix, the direct dense-matrix method, which is the fastest for small $ℓₘₐₓ$ and lets the sample points be chosen freely, with round-trip errors of about $10^{-14}$ for $ℓₘₐₓ ≲ 24$, somewhat larger than those of "RS";
• "Minimal": SSHTMinimal, the optimal-dimensionality algorithm of Elahi et al., which uses exactly as many samples as there are modes; mostly experimental, not very accurate.

The remaining keyword arguments are passed to the constructor of the chosen type.

Mode weights

A ModeWeights is accepted wherever mode weights are, and its labels are checked against the transform: weights of another spin weight are refused, even when they have the right number of entries. Synthesis, 𝒯 * w, accepts weights of the transform's spin weight with any range of $ℓ$ up to the transform's $ℓₘₐₓ$: the modes that w lacks are zero, and its entries with $ℓ < |s|$, which belong to no function of spin weight $s$, are ignored. The three-argument mul!(f, 𝒯, w) and the outputs of analysis need exactly the range $|s| ≤ ℓ ≤ ℓₘₐₓ$; ModeWeights(w; ℓₘᵢₙ, ℓₘₐₓ) copies mode weights into another range. Analysis of one-dimensional data returns a ModeWeights, and of more dimensions a plain array.

Acting in place

The "RS" method never acts in place: 𝒯 * f̃ and 𝒯 \ f allocate their results and leave their arguments alone. The other two methods have an inplace keyword argument, which is part of the type of the resulting object, and which makes them reuse the storage of their argument. It is on by default for both — for "Matrix" whenever there are exactly as many sample points as modes, as there are with its default points — so that whether 𝒯 \ f preserves f depends on the method; pass inplace=false for the behavior of "RS". With inplace=true,

• "Matrix" acts in place for analysis only: 𝒯 \ f overwrites f with the mode weights, while 𝒯 * f̃ always allocates its result, since a matrix product cannot write over its own input;
• "Minimal" acts in place in both directions: 𝒯 * f̃ overwrites f̃ with the function values, and 𝒯 \ f overwrites f with the mode weights. Analysis in place of one-dimensional data returns a ModeWeights wrapping the argument's own storage, now holding the mode weights, and synthesis in place returns the storage itself, now holding the function values, as a plain array — not the ModeWeights that may have held the input, whose labels no longer describe it. Mode weights whose range of $ℓ$ differs from the transform's are first copied into its range, and it is the copy that synthesis overwrites and returns. Acting in place needs storage that can hold the complex results, and real storage is refused with an ArgumentError; real data are transformed with inplace=false.

Whatever the option, the two-argument LinearAlgebra.ldiv!(𝒯, x) analyzes in place for "Matrix" and "Minimal", returning the mode weights as 𝒯 \ x would, and the two-argument LinearAlgebra.mul!(𝒯, x) synthesizes in place for "Minimal", returning the function values as a plain array. "RS" has neither, because its number of sample points differs from the number of modes, and "Matrix" has no two-argument mul!. The three-argument mul!(f, 𝒯, f̃) and ldiv!(f̃, 𝒯, f) write into the given output for every type, and never touch their input. The output of ldiv!(f̃, 𝒯, f) for one-dimensional data may be a bare vector at least nmodes(𝒯) long, and the mode weights then come back as a ModeWeights over its first entries, rather than as the vector.

Sample data

Sample data given as keyword arguments — the rotors Rθϕ of "Matrix", and the colatitudes θ and quadrature_weights of "RS" and "Minimal" — must already be in the type T that the transform works in: rotors as Rotor{T}, and real numbers of type T, or integers, which convert exactly. Anything else is rejected, rather than converted. In particular a general Quaternion or a QuatVec is not taken for a rotor, and BigFloat data is not rounded to the default T=Float64; pass BigFloat as T to work in that type.

Use from several tasks

Each transform runs on the thread of the task that calls it, including a transform of many columns of data at once. Parallelism is therefore a matter of running transforms in several tasks at the same time, and the three types differ in what that needs.

The "RS" and "Minimal" objects hold workspace, which every transform overwrites, so one object must never be used by two tasks at the same time. Give each task its own: copy(𝒯) returns an independent transform that shares the read-only tables and FFT plans of 𝒯 and allocates new workspace, at a small fraction of the cost of constructing one. A vector of objects indexed by Threads.threadid() is no substitute: a task may move to another thread whenever it yields, and a second task, then running on the thread it left, is given the same object. A pattern that works is to divide the data into chunks and to spawn one task per chunk, each with its own copy:

chunks = Iterators.partition(maps, cld(length(maps), Threads.nthreads()))tasks = [Threads.@spawn(let 𝒯ₖ = copy(𝒯); [𝒯ₖ \ f for f ∈ chunk]; end) for chunk ∈ chunks]mode_weights = reduce(vcat, fetch.(tasks))

The "Matrix" object holds no workspace — its transforms only read the matrix and its decomposition, through BLAS and LAPACK — so one object may be used by any number of tasks at once, and its copy is the object itself.

A deepcopy of a transform is also independent, and shares the FFT plans of the original as a copy does, since a plan is never modified once it is made. A transform may be serialized, as it is when it is sent to another process with Distributed; the FFT plans, which belong to the process that made them, are then made again where the transform is deserialized.

Half-integer indices

The spin weight and $ℓₘₐₓ$ may be half-integers, passed as Rationals with denominator 2 — SSHT(1//2, 7//2) — for the "RS" and "Matrix" methods; the "Minimal" method is defined only for integer spin weights, and says so. The mode weights are then indexed by half-odd $ℓ$ and $m$ in the same canonical ordering, and the function values are, as before, the values at the points that rotors returns. Those points deserve more attention than they need for an integer spin weight. A function of half-integer spin weight is a function on the double cover of the rotation group rather than on the sphere, so its value "at a point of the sphere" depends on which of the two rotors above that point is meant; the rotors used here are from_spherical_coordinates(θ, ϕ), and the function values are therefore antiperiodic in the azimuth — a full circuit in $ϕ$ arrives at the antipodal rotor, and changes the sign. Each ring needs $N_ϕ ≥ 2ℓₘₐₓ+1$ points, as for an integer $ℓₘₐₓ$, and the default $2ℓₘₐₓ+1$ rings of $2ℓₘₐₓ+1$ points are even numbers when $ℓₘₐₓ$ is a half-odd-integer. A quadrature rule symmetric about the equator, such as Fejér's or Clenshaw–Curtis, then needs only $2ℓₘₐₓ$ rings, one fewer than it would for an integer $ℓₘₐₓ$.

pixels(𝒯)

Return the spherical coordinates (θ, ϕ) at which the transform 𝒯 evaluates functions, in the order used by the first dimension of the function values. See also rotors.

rotors(𝒯)

Return the Rotors at which the transform 𝒯 evaluates functions, in the order used by the first dimension of the function values. The vector is a new one, which may be modified without affecting the transform. See also pixels.

SSHTMatrix(s, ℓₘₐₓ, [T=Float64]; Rθϕ=leja_rotors(s, ℓₘₐₓ, T), decomposition, inplace)

Construct an $s$-SHT object that uses the "Matrix" method. The dense matrix of harmonics sYlm_matrix evaluated at the rotors Rθϕ is stored, so that synthesis is a matrix-vector product and analysis is a linear solve using the given decomposition of that matrix. Also see SSHT for general information about how to use these objects.

The parameters of the type SSHTMatrix{T, Inplace, Tdecomp, IT} are as follows:

• T is the real type the transform works in.
• Inplace is true when analysis acts in place, as set by the inplace keyword.
• Tdecomp is the type of the decomposition of the matrix.
• IT is the index type, Int or HalfOddInteger.

By default, this uses precisely optimal sampling — meaning that the number of points on which the function is evaluated (the length of Rθϕ), is equal to the number of modes $(ℓₘₐₓ+1)²$. However, it is equally possible to evaluate on more points than there are modes, in which case the analysis is a least-squares solve. 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. The rotors must be Rotor{T}s, in the type T that the transform works in.

With precisely optimal sampling, the accuracy is limited by the conditioning of the points. The default points, those of leja_rotors, are chosen to keep the matrix well conditioned: with the default LU decomposition, a round trip of random mode weights loses about 2.5 digits (in the norm of the error relative to that of the weights) at ℓₘₐₓ = 32, and 3 at ℓₘₐₓ = 64. Choosing them costs about 2–3 times as much as the decomposition itself. Points that merely spread evenly over the sphere are not enough: on the golden-ratio spiral of golden_ratio_spiral_rotors, a round trip loses about 4 digits at ℓₘₐₓ = 16, 5 at 32, and all of them at 64. The constructor measures the error of one round trip, and warns when fewer than half the digits of T would survive.

LU decomposition (LinearAlgebra.lu) is used by default for a square matrix; any function that decomposes the matrix into something capable of solving the linear problem may be passed as decomposition. In particular, QR decomposition (LinearAlgebra.qr) will typically be about twice as slow, with errors a few times smaller; it is also the default when there are more points than modes, where an LU decomposition, which cannot solve a least-squares problem, is rejected.

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. It applies to analysis only: with inplace=true, 𝒯 \ f overwrites the storage of f with the mode weights, and returns them as a ModeWeights wrapping it, while 𝒯 * f̃ always allocates a new array for the function values, since a matrix product cannot write over its own input. The two-argument LinearAlgebra.ldiv!(𝒯, f) analyzes in place whatever the option, including for more points than modes, when the mode weights are left in the first nmodes(𝒯) entries of f; there is no two-argument mul!. Real function values are analyzed with inplace=false, since in place the complex mode weights would have to be written into real storage. See SSHT for a description of in-place operation.

The object holds no workspace: its transforms only read the matrix and its decomposition, through BLAS and LAPACK, so one object may be used by any number of tasks at the same time, and copy(𝒯) is 𝒯 itself.

This method is fastest for $ℓₘₐₓ ≲ 24$, where its round-trip errors are about $10^{-14}$, somewhat larger than those of the "RS" method. However, this advantage quickly falls away. A warning is issued if ℓₘₐₓ is greater than about 64, because this method is not likely to be the most efficient or most accurate choice.

SSHTMinimal(s, ℓₘₐₓ, [T=Float64]; θ=minimal_rings(s, ℓₘₐₓ, T).θ, plan_fft_flags=FFTW.ESTIMATE, plan_fft_timelimit=Inf, inplace=true)

Construct an $s$-SHT object that uses the optimal-dimensionality algorithm described by Elahi et al., which samples the function at exactly as many points as there are modes. This may also be achieved by calling the main SSHT function with the same keywords, along with method="Minimal".

The parameters of the type SSHTMinimal{T, Inplace, P, BP} are as follows:

• T is the real type the transform works in.
• Inplace is true when the transforms act in place, as set by the inplace keyword.
• P and BP are the types of the forward and backward FFT plans.

The function is sampled on $ℓₘₐₓ-|s|+1$ "rings" at constant colatitude, each holding an odd number of equally spaced points starting at $ϕ = 0$. Their sizes, and the default colatitudes, are given by minimal_rings; the θ keyword argument may give other colatitudes, which must be distinct, one for each ring in the order listed there, which is also the order of the sample points, given in the type T that the transform works in, or as integers. See pixels and rotors for the sample points themselves.

For $s = 0$ the rings have $1, 3, …, 2ℓₘₐₓ+1$ points. For any other spin weight that choice is badly conditioned: near the north pole a function of spin weight $s$ is dominated by the modes with $m$ near $-s$, and near the south pole by those near $+s$, so that the smallest rings — which are placed nearest the poles — see some of the frequencies they are responsible for more weakly than the higher frequencies that alias onto them. The error of a round trip then grows by more than an order of magnitude with each unit of ℓₘₐₓ, and no choice of colatitudes cures it. For $s ≠ 0$ the rings are therefore arranged so that each polar ring is centered, in frequency, on the modes that dominate near its pole; see minimal_rings. The analysis is then no longer a sequence of solves for one $m$ at a time, but for small groups of $m$ values that alias into one another — in every case measured, at most $4|s|-1$ of them, however large ℓₘₐₓ is.

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, and are passed to AbstractFFTs.plan_fft!. One pair of plans is made for each distinct number of points on a ring, and each plan runs on a single thread.

Because this algorithm achieves optimal dimensionality, the transformation is performed in place by default, in both directions: 𝒯 * f̃ overwrites the storage of f̃ with the function values (and returns that storage, as a plain array), and 𝒯 \ f overwrites f with the mode weights (and returns them, for one-dimensional f, as a ModeWeights wrapping that storage). That storage must hold complex floating-point numbers. If this is not desired, pass the keyword argument inplace=false, which makes those operations work on a copy of the input, in Complex{T}, and so accept real data. Whatever the option, the two-argument LinearAlgebra.mul!(𝒯, x) and LinearAlgebra.ldiv!(𝒯, x) act in place, the first returning the function values as a plain array and the second the mode weights as 𝒯 \ x would. See SSHT.

The values ${}_sλ_{ℓ,m}(θ_r)$ of every mode on every ring are precomputed at construction (with one batched sλlmCalculator) and stored, which takes $O(ℓₘₐₓ^3)$ memory, as are the LU decompositions of the matrices for the groups of $m$ values. The object also holds workspace for the transforms, so two tasks must not use it at the same time; copy(𝒯) returns another transform, which shares those tables and the FFT plans of 𝒯 and has workspace of its own. See SSHT.

Even so, the sample points become increasingly badly conditioned as ℓₘₐₓ grows, for every spin weight: in Float64 a round trip loses about 3.5 digits by ℓₘₐₓ = 32, 5.5 by ℓₘₐₓ = 48 and 8 by ℓₘₐₓ = 64 at $s = 0$, and more at larger $|s|$ — about 9 by ℓₘₐₓ = 32 at $s = 2$. The constructor therefore measures the error of one round trip, and warns when fewer than half the digits of T would survive; the "RS" method has no such limitation.

This method is defined only for integer spin weights. Its bookkeeping — rings of an odd number of points, and the aliasing of $m$ into rings too small to hold it — is written for integer indices throughout, and has not been extended; a half-integer spin weight is rejected with a message naming the two methods, "RS" and "Matrix", that do accept one.

SSHTRS(s, ℓₘₐₓ, [T=Float64]; [θ, quadrature_weights], Nϕ=2ℓₘₐₓ+1, plan_fft_flags=FFTW.ESTIMATE, plan_fft_timelimit=Inf)

Construct an $s$-SHT object that uses the ring-based algorithm described by Reinecke and Seljebotn. This may also be achieved by calling the main SSHT function with the same keywords, along with method="RS" (the default).

The parameters of the type SSHTRS{T, ST, P, BP, B, IT} are as follows:

• T is the real type the transform works in.
• ST is the storage type of the $H$ wedge in the transform's sλlmCalculator.
• P and BP are the types of the forward and backward FFT plans.
• B is true when that calculator is batched (see isbatched).
• IT is the index type, Int or HalfOddInteger.

The spin-weighted spherical harmonics are evaluated on a series of "rings" at constant colatitude, whose locations are given by the θ keyword argument, and the analysis integrates over $θ$ with the quadrature_weights of the rule that placed those rings. When both are omitted they are the nodes and weights of Fejér's first rule, fejer1_rings(2ℓₘₐₓ+1, T) and fejer1(2ℓₘₐₓ+1, T). Only the caller knows which rule placed a given set of rings, so the two must be given together — for example θ=clenshaw_curtis_rings(N, T) with quadrature_weights=clenshaw_curtis(N, T) — and either one without the other is refused. The analysis is exact for band-limited functions when the quadrature rule integrates polynomials of degree $2ℓₘₐₓ$ in $\cos θ$ exactly, as the Fejér and Clenshaw–Curtis rules with at least $2ℓₘₐₓ+1$ nodes do. The constructor checks this, and warns when the rule falls short; synthesis does not use the weights, and is exact on any rings.

On each ring, an FFT is performed. To reach the band limit of $m = ±ℓₘₐₓ$, the number of points along each ring must 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. 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 used for all rings.

The sample points are ordered ring by ring, with the azimuth $ϕ_k = 2πk/N_ϕ$ for $k = 0, …, N_ϕ-1$ varying fastest; see pixels and rotors. The transform works in the floating-point type T, and the colatitudes and weights must be given in that type, or as integers. It never acts in place: 𝒯 * f̃ and 𝒯 \ f allocate their results, and real data are accepted; see SSHT.

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, and are passed to AbstractFFTs.plan_fft!. One pair of plans is made for each distinct number of points on a ring, and each plan runs on a single thread, since a ring is too short for FFTW's threads to be worth their cost. For other element types the ring FFTs are computed by generic code, and for BigFloat they are most of the cost of a transform; one transform still runs on a single thread, so the way to use several threads is to give each task its own copy(𝒯).

The harmonics on the rings are computed with one batched sλlmCalculator, one $ℓ$ at a time, so the cost is $O(N_θ ℓₘₐₓ^2)$ and the memory $O(N_θ ℓₘₐₓ)$. Data of several columns are transformed up to eight columns at a time, with the harmonics of each $ℓ$ computed once for all of them, which makes each column about 1.5 to 2 times faster to transform; the results of each column are exactly those of a transform of that column alone. The object holds workspace for these computations, so two tasks must not use it at the same time; copy(𝒯) returns another transform, which shares the rings, weights and FFT plans of 𝒯 and has workspace of its own. See SSHT.

Half-integer s and ℓₘₐₓ are accepted as Rationals with denominator 2; see SSHT for what the function values then mean. The algorithm is the same one: the harmonics on a ring are $i^{2s}$ times real functions of $θ$, so the $θ$ stage works with those real functions and restores the phase once per ring, and $e^{imϕ}$ with half-odd $m$ is $e^{iϕ/2}$ times a Fourier mode of integer frequency $m - 1/2$, so the FFT runs at integer frequencies and each sample of a ring is multiplied by $e^{±iϕ/2}$. The defaults are $2ℓₘₐₓ+1$ rings of $2ℓₘₐₓ+1$ points, as for integer indices, and these are then even numbers.

map2salm(map, s, ℓₘₐₓ)
map2salm(map, 𝒯::SSHTRS)

Transform function values map sampled on an equiangular grid to spin-weighted spherical-harmonic mode weights ${}_sa_{ℓ,m}$.

The map array, of real or complex numbers, must have size $N_ϕ$ along its first dimension and $N_θ ≥ 2$ along its second, with any number of dimensions following, sampled at the Clenshaw–Curtis nodes clenshaw_curtis_rings(Nθ) in $θ$ (which include both poles) and at $ϕ_k = 2πk/N_ϕ$; pixels returns that grid for a constructed SSHTRS. For the analysis to be exact for a band-limited function, one needs $N_ϕ ≥ 2ℓₘₐₓ+1$ points on each ring and $N_θ ≥ 2⌊ℓₘₐₓ⌋+1$ rings — that is, $2ℓₘₐₓ+1$ rings for an integer $ℓₘₐₓ$, and $2ℓₘₐₓ$ for a half-odd one; with fewer, the result is not exact, and a warning is issued. The transform works in the floating-point type of the map's elements, float(real(eltype(map))).

The result is a ModeWeights for a one-dimensional map ($N_ϕ × N_θ$), or an array whose first dimension indexes the modes in the canonical ordering ℓ ∈ abs(s):ℓₘₐₓ, m ∈ -ℓ:ℓ otherwise. Note that the output starts at $ℓ = |s|$, not at $ℓ = 0$.

For repeated use with different maps of the same shape, construct the transform once with SphericalFunctions.map2salm_plan(map, s, ℓₘₐₓ) (an SSHTRS on the Clenshaw–Curtis rings) and pass it as the second argument. See also salm2map.

The spin weight and $ℓₘₐₓ$ may be half-integers, passed as Rationals with denominator 2. The map is then antiperiodic in $ϕ$ — its values are those of the function at the rotors from_spherical_coordinates(θ, ϕ), and a full circuit of the azimuth reaches the antipodal rotor — and each ring needs $N_ϕ ≥ 2ℓₘₐₓ+1$ points, as for an integer $ℓₘₐₓ$; see SSHT.

map2salm_plan(map, s, ℓₘₐₓ)

Construct the SSHTRS transform used by map2salm and salm2map for maps of the shape of map ($N_ϕ × N_θ × …$) on the Clenshaw–Curtis grid, in the floating-point type of the map's elements.

salm2map(salm, s, ℓₘₐₓ, Nϕ, Nθ)
salm2map(w::ModeWeights, Nϕ, Nθ)
salm2map(salm, 𝒯::SSHTRS)

Evaluate the spin-weighted function with mode weights salm on the equiangular $N_ϕ × N_θ$ grid used by map2salm. The mode weights must be given in the canonical ordering ℓ ∈ abs(s):ℓₘₐₓ, m ∈ -ℓ:ℓ along the first dimension, as real or complex numbers, or as a ModeWeights, which may then cover any range of $ℓ$ up to $ℓₘₐₓ$ (see SSHT). Given a ModeWeights w, the spin weight and $ℓₘₐₓ$ may be omitted, and are then those of w. The result has size $N_ϕ × N_θ$ followed by the trailing dimensions of salm, and is computed in the floating-point type of the elements of salm.

Synthesis is exact on any grid of $N_θ ≥ 2$ rings, so no condition on $N_ϕ$ or $N_θ$ is needed beyond that; but a map with fewer than $2ℓₘₐₓ+1$ points on each ring, or fewer rings than map2salm needs, cannot be analyzed back to the same mode weights.

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

Cover the sphere 𝕊² with pixels given by the Driscoll and Healy [12] equiangular grid. They write:

Let $f(θ, ϕ)$ be a band-limited function such that $\hat{f}(l, m) = 0$ for $l ≥ b$. We sample the function at the equiangular grid of points $(θ_i, ϕ_j)$, $i = 0, \ldots, 2b-1$, $j = 0, \ldots, 2b-1$, where $θ_i = π i/2b$ and $ϕ_j = π j/b$.

The band limit here is $b = ℓₘₐₓ + 1$, and ℓₘₐₓ must be a non-negative integer. The grid is defined for integer indices only; functions of half-integer spin are sampled with golden_ratio_spiral_pixels, leja_pixels or sorted_ring_pixels instead.

These are the $(2b)^2$ points returned, ordered with $ϕ$ varying fastest. The first ring, $i = 0$, lies at the north pole $θ = 0$, and all $2b$ of its points are included, as distinct rotors that differ by a rotation about the pole. The grid does not reach the south pole, since the largest colatitude is $π(2b-1)/2b$.

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

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

Cover the sphere 𝕊² with pixels given by the Driscoll and Healy [12] equiangular grid.

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

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 a common angular pattern of phyllotaxis in plants, found in sunflowers, pine cones, and artichokes, for example. Wrapped around 𝕊², 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 number of pixels is Ysize(abs(s), ℓₘₐₓ), the number of modes of spin weight s with $|s| ≤ ℓ ≤ ℓₘₐₓ$, which is $(ℓₘₐₓ+1)^2 - s^2$ for either kind of index. The spin weight and ℓₘₐₓ may be integers or half-odd-integers, the latter passed as Rationals with denominator 2, as in golden_ratio_spiral_pixels(1//2, 7//2). The two must be of one kind; a call that mixes them, such as golden_ratio_spiral_pixels(1//2, 3), is an error, and so is a spin weight with $|s| > ℓₘₐₓ$, for which there are no modes.

Pixel $k = 0, 1, …, N-1$ of the $N$ pixels has $\cos θ = 1 - (2k+1)/N$, the center of the $k$th of $N$ bands of equal width in $\cos θ$, and so of equal area. Its azimuth is $kΔϕ$ reduced modulo $2π$, so that every $ϕ$ lies in $[0, 2π)$. Both coordinates are formed from exact integer quantities — the reduction modulo $2π$ is made in exact arithmetic — so that their only errors are those of the final few roundings, of about an ulp of T, and the pixels in different types T are roundings of the same points.

This is also known as the "Fibonacci sphere" or "Fibonacci lattice" — though Fibonacci had nothing to do with it; later authors pointed out relationships between his sequence and the golden ratio.

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.

Note that these points can be constructed analytically, as is done here, but do not lead to very well conditioned matrices of spin-weighted spherical harmonics. A somewhat better choice is given by leja_pixels, which chooses a subset of the points of a denser golden-ratio spiral with better conditioning, and is the default for SSHTMatrix.

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, including the half-odd-integer values that s and ℓₘₐₓ may take. The quantity returned by this function is a vector of Rotors providing each pixel.

leja_pixels(s, ℓₘₐₓ, [T=Float64]; oversampling=2)

Choose Ysize(abs(s), ℓₘₐₓ) pixels — as many as there are modes of spin weight s with $|s| ≤ ℓ ≤ ℓₘₐₓ$ — so that the matrix of spin-weighted spherical harmonics at those points is well conditioned. This is needed to solve for the mode weights from function values at exactly that many points.

The golden-ratio spiral of golden_ratio_spiral_pixels spreads its points evenly, but that is not what such a solve needs, and with exactly as many points as modes the matrix becomes badly conditioned as ℓₘₐₓ grows: its condition number is about $10^4$ at ℓₘₐₓ = 16, $10^5$ at 32 and $10^{14}$ at 64. The points chosen here are "discrete Leja points" [13], which approximate the Fekete points that maximize the determinant of the matrix: the harmonics are evaluated on a golden-ratio spiral of oversampling times as many candidate points, and an LU decomposition with partial pivoting of that tall matrix picks, one row at a time, the candidate that is least well represented by those already picked. With the default oversampling=2, the condition number is about 30 at ℓₘₐₓ = 16, 70–100 at 32 and 200–400 at 64, for either $s = 0$ or $s = 2$, and a round trip through SSHTMatrix on these points loses about 3 digits at ℓₘₐₓ = 64 rather than all of them. (Up to ℓₘₐₓ ≈ 12 the two are comparable, with condition numbers of order 10, and the spiral is sometimes slightly the better; beyond that it degrades and these points do not.) Larger values of oversampling, which must be finite and at least 1, did not help consistently in tests, and cost more.

The points depend on s, not just on their number, because the matrix does. They are a subset of that spiral, returned in its order (from the north pole to the south), and are deterministic. The candidates are chosen in Float64, by LAPACK's LU decomposition of the harmonics at the Float64 spiral, whatever T is, and only the chosen candidates are then evaluated in T; so every T gets the same points, each rounded to T, and a function sampled in one precision can be analysed in another. The cost is that of the LU decomposition of an $N_c × N$ matrix, with $N$ the number of modes and $N_c$ the number of candidates — about 2–3 times that of the decomposition SSHTMatrix itself performs in Float64, or 2–3 seconds at ℓₘₐₓ = 64 — and its $O(N^2)$ storage. These are the default sample points of SSHTMatrix.

The spin weight and ℓₘₐₓ may be integers or half-odd-integers, as for golden_ratio_spiral_pixels, on the same terms. The returned quantity is a vector of 2-SVectors providing the spherical coordinates of each pixel. See also leja_rotors for the corresponding Rotors.

leja_rotors(s, ℓₘₐₓ, [T=Float64]; oversampling=2)

Choose pixels on which the spin-weighted spherical harmonics are well conditioned, as discrete Leja points drawn from a golden-ratio spiral.

See leja_pixels for more detailed explanation, including the half-odd-integer values that s and ℓₘₐₓ may take. The quantity returned by this function is a vector of Rotors providing each pixel.

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

Cover the sphere 𝕊² with pixels given by the McEwen and Wiaux [14] equiangular grid. They consider "signals on the sphere bandlimited at $L$, that is signals such that ${}_sf_{ℓm} = 0, ∀ℓ ≥ L$", so that $L = ℓₘₐₓ + 1$ here, and write:

We adopt an equiangular sampling of the sphere with sample positions given by $θ_t = \frac{π(2t+1)}{2L-1}$, where $t ∈ \{0, 1, …, L-1\}$ […] and $ϕ_p = \frac{2πp}{2L-1}$, where $p ∈ \{0, 1, …, 2L-2\}$.

The band limit ℓₘₐₓ must be a non-negative integer. The grid is defined for integer indices only; functions of half-integer spin are sampled with golden_ratio_spiral_pixels, leja_pixels or sorted_ring_pixels instead.

These are the $L(2L-1)$ points returned, ordered with $ϕ$ varying fastest. The last ring, $t = L-1$, lies at the south pole $θ = π$, and all $2L-1$ of its points are included, as distinct rotors that differ by a rotation about the pole; McEwen and Wiaux count the pole once, for $(L-1)(2L-1)+1$ samples. Their algorithm also extends the $θ$ domain to $[0, 2π)$, reflecting the samples through the south pole to $t ∈ \{L, …, 2L-2\}$, but those are not further sample positions on the sphere, and are not returned.

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

Cover the sphere 𝕊² with pixels given by the McEwen and Wiaux [14] equiangular grid.

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

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

The rings on which SSHTMinimal samples a function of spin weight s band-limited at ℓₘₐₓ, as a named tuple (; Nϕ, centers, θ). Its entries are the number of points on each ring, the center of each ring's window of frequencies, and the default colatitude of each ring. The rings are listed in order of size (for rings of equal size, the northern first), and there are $ℓₘₐₓ-|s|+1$ of them, holding $(ℓₘₐₓ+1)^2 - s^2$ points in all — exactly the number of modes.

A ring of $N = 2k+1$ points cannot distinguish frequencies $m$ that differ by a multiple of $N$; the analysis treats its Fourier coefficients as measuring the $N$ consecutive frequencies $|m - c| ≤ k$ of its window, centered on $c$, and removes the aliases of all other frequencies from them. For each $m$ there must be as many rings whose windows include $m$ as there are modes with that $m$, namely $ℓₘₐₓ - \max(|m|, |s|) + 1$. The windows $|m| ≤ a$ for $a ∈ |s|:ℓₘₐₓ$, one ring for each, satisfy this, and are what is used for $s = 0$, with the colatitudes of sorted_rings.

For $s ≠ 0$ that arrangement is badly conditioned (see SSHTMinimal). Instead, pairs of those windows are recentered. The two windows $|m| ≤ a$ and $|m| ≤ a+2d$ cover every $m$ exactly as often as the two windows $|m + d| ≤ a+d$ and $|m - d| ≤ a+d$ do, and the latter become a ring of $2(a+d)+1$ points in the northern hemisphere, whose window is centered on $-d\,\mathrm{sign}(s)$, and one of the same size in the southern hemisphere, centered on $+d\,\mathrm{sign}(s)$ — toward the frequencies $∓s$ that dominate near each pole. The pairs are chosen greedily, first with $d = |s|$ and $a$ in increasing order, whenever both windows are still available, and then with successively smaller $d$ down to 1, which matters when $ℓₘₐₓ < 3|s|$ and no window has a partner $2|s|$ larger. The windows left unpaired — the largest ones — remain centered on 0, and become the rings nearest the equator. The default colatitudes are equally spaced, $θ = iπ/(n+1)$ for $i ∈ 1:n$ with $n$ the number of rings; the northern rings take the slots nearest the north pole, smallest first, and likewise in the south, and the rings centered on 0 take the remaining slots in the order sorted_rings uses. For $s = 0$ this reproduces sorted_rings exactly. The order is decided by the positions of the slots, not by their rounded colatitudes, so the rings are the same for every T.

Like SSHTMinimal, this is defined for integer indices only.

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

Cover the sphere 𝕊² with Ysize(abs(s), ℓₘₐₓ) pixels, distributed in the rings provided by sorted_rings. That count is the number of modes of spin weight s with $|s| ≤ ℓ ≤ ℓₘₐₓ$, which is $(ℓₘₐₓ+1)²-s²$; see the documentation of sorted_rings for more description, including the half-odd-integer values that s and ℓₘₐₓ may take.

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 Ysize(abs(s), ℓₘₐₓ) pixels, distributed in the rings provided by sorted_rings. That count is the number of modes of spin weight s with $|s| ≤ ℓ ≤ ℓₘₐₓ$, which is $(ℓₘₐₓ+1)²-s²$; see the documentation of sorted_rings for more description, including the half-odd-integer values that s and ℓₘₐₓ may take.

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

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

Compute locations of a series of rings labelled by $j ∈ |s|:ℓₘₐₓ$ (analogous to $ℓ$). 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, and for $s = 0$ it gives the default rings of the Minimal algorithm. For $s ≠ 0$ these rings are badly conditioned, and the Minimal algorithm uses the different arrangement of minimal_rings instead; see SSHTMinimal.

The rings occupy the $n = ℓₘₐₓ - |s| + 1$ interior points of an equally spaced grid on $[0, π]$, and are assigned to them in order of size, from the slots farthest from the equator to the nearest. Two slots that are mirror images of each other across the equator are equally far from it, and of the two rings assigned to such a pair, the larger goes to the northern slot for $s ≥ 0$ and to the southern for $s < 0$. The order is decided by the positions of the slots in the grid, not by their rounded colatitudes, so it is the same for every T, and for $s ≠ 0$ the rings for $-s$ are exactly the mirror images of those for $s$. (For $s = 0$ the larger ring of each mirror pair is the northern one, so the rings are not symmetric about the equator.)

The spin weight and ℓₘₐₓ may be integers or half-odd-integers, the latter passed as Rationals with denominator 2, as in sorted_rings(1//2, 7//2); the two must be of one kind. There are then $ℓₘₐₓ - |s| + 1$ rings, and ring $j$ has $2j+1$ pixels, both of which are whole numbers for either kind of index. A spin weight with $|s| > ℓₘₐₓ$ describes no modes, and is an error.

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

\[θ_k = k \frac{π}{n-1} \quad k=0, \ldots, n-1,\]

as given by clenshaw_curtis_rings.

The weights approximate $∫_0^π f(θ) \sin θ\, dθ = ∫_{-1}^{1} f(x)\, dx$, where $x = \cos θ$, so that they sum to 2. This function uses Waldvogel's method.

The number of nodes n must be at least 2, since the nodes include both poles. The type T may be any AbstractFloat, but defaults to Float64; the weights are returned as a Vector{T}.

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. The number of rings N must be at least 2, since the rings include both poles.

fejer1(n, [T])

Compute n weights for Fejér's first rule, corresponding to n evenly spaced nodes between 0 and π exclusive, offset by half a spacing from each end. That is, the nodes are located at

\[θ_k = (2k+1) \frac{π}{2n} \quad k=0, \ldots, n-1,\]

as given by fejer1_rings.

The weights approximate $∫_0^π f(θ) \sin θ\, dθ = ∫_{-1}^{1} f(x)\, dx$, where $x = \cos θ$, so that they sum to 2. This function uses Waldvogel's method.

The number of nodes n must be at least 1. The type T may be any AbstractFloat, but defaults to Float64; the weights are returned as a Vector{T}.

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. The number of rings N must be at least 1.

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

\[θ_k = k \frac{π}{n+1} \quad k=1, \ldots, n,\]

as given by fejer2_rings.

The weights approximate $∫_0^π f(θ) \sin θ\, dθ = ∫_{-1}^{1} f(x)\, dx$, where $x = \cos θ$, so that they sum to 2. 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 number of nodes n must be at least 1. The type T may be any AbstractFloat, but defaults to Float64; the weights are returned as a Vector{T}.

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. The number of rings N must be at least 1.