API Reference

make_plan([FE = Float64,] θ_nodes, φ_nodes, lmax; tol=1e-8, ntrans=1, tuning=NoTuning(), …)

FE is the field element type, positional as it is for zeros(T, …), and it selects the transform: Float64/Float32 build the real specialization (Hermitian torus spectrum → rfft/brfft, half the θ wavenumbers reach the NUFFT), ComplexF64/ComplexF32 the complex one (fft/bfft, full spectrum). Both are the same spherical harmonic transform; only the symmetry exploited differs. The T = FE keyword form forwards to the positional one.

Construct a NUSHTplan for M scattered points at colatitudes θ_nodes ∈ [0,π] and longitudes φ_nodes ∈ [0,2π), up to spherical harmonic degree lmax, transforming ntrans co-located fields per call. Builds the FINUFFT guru plans and sets the nonuniform points once; they are freed by a finalizer (or eagerly via close!).

Keyword arguments:

• tol: FINUFFT accuracy tolerance.
• T: floating-point type (Float64/Float32).
• ntrans: batch size B — transform B co-located fields (same nodes) per call.
• tuning: an AbstractPlanTuning — NoTuning (default), AutoTuning or ThoroughTuning. The default already pins the settings that matter most; the searching strategies trade O(100 ms)-O(1 s) of trial builds for roughly a further 1.1x per transform, so they are worth it only for a plan reused thousands of times. Outcomes are memoized per problem shape, so building many plans of one size pays the search once.
• ft_fftw_nthreads, ft_fftw_flags: override the FFTW planner thread count / flags used for the sphere synthesis and analysis plans.
• nthreads, upsampfac: override FINUFFT's thread count (0 is its "all cores" sentinel) and upsampling factor.

Each override keyword defaults to nothing, meaning "whatever tuning decides". An explicit value is honoured exactly and skips the search for that setting.

FINUFFT accepts coordinates in [-3π, 3π], so natural [0,π]/[0,2π) coordinates are passed directly.

NUSHTplan{T}

Pre-computed plan for non-uniform spherical harmonic transforms at M scattered points, up to degree lmax, transforming B co-located fields per call (ntrans = B).

Fields:

• lmax, Nθ = lmax+1, Nφ = 2lmax+1, B (batch size / FINUFFT ntrans)
• tol: FINUFFT accuracy tolerance
• nodes: the AbstractNodeSet holding θ_nodes, φ_nodes (colatitudes ∈ [0,π] and longitudes ∈ [0,2π) of the M points), the (M, B) strengths buffer fbuf, and the two FINUFFT guru plans. set_nodes! re-points it.
• C, F: coefficient scratch and the bivariate Fourier coefficients P·C, both (Nθ, Nφ, B) with eltype FE
• Fhat: the (2lmax+1, 2lmax+1, B) complex mode array the NUFFT evaluates, assembled from F by _assemble_modes!. Both axes carry wavenumbers -lmax…lmax.
• Fslice: (Nθ, Nφ) scratch Matrix each batch slice is copyto!-ed through for the S-step. FastTransforms has no Float32 sphere plan, so it stays double precision and the copy converts.
• sph_plan, sph_plan_adj: FastTransforms plan_sph2fourier (P) and its adjoint on a (Nθ, Nφ) slice. P alone is the whole S-step — its output is already the bivariate Fourier series, so synthesis is P·C → assemble → one NUFFT, with no equiangular grid in between.

nodes.nufft_type2 is the guru type-2 plan (iflag = +1, synthesis N) and nodes.nufft_type1 the type-1 plan (iflag = -1, adjoint N†). iflag = +1 supplies the reconstruction sign directly, so the modes need no conjugate-transpose. The axis convention is x = θ, y = φ, and both mode axes are in centered order (modeord = 0), which is what -lmax…lmax already is.

close!(plan::NUSHTplan)

Eagerly free the FINUFFT guru plans owned by plan (otherwise freed by their finalizers). Safe to call more than once (finufft_destroy! is idempotent).

set_nodes!(plan, θ_nodes, φ_nodes) -> plan

Move a plan's nodes to new positions, reusing every structure fixed by (lmax, B, T) — the FastTransforms sphere plans, the FFTW plans and all coefficient buffers. Only FINUFFT's point tables are rebuilt, so this costs 0.01–0.10 ms against 3–7 ms to build an equivalent plan from scratch (measured, lmax 45–128).

The points may move anywhere. Changing how many there are additionally needs the plan's point-indexed buffers to be replaced, which requires make_plan(…; variable_npts = true); a plan built with the default FixedCountNodes throws instead of silently reallocating.

With the count unchanged this rewrites array contents only and allocates nothing, on either node-set form.

AbstractNodeSet

The point-dependent half of a plan: the M scattered nodes, the (M, B) strengths buffer, and the two FINUFFT guru plans — which own the loaded point tables, so they belong with the points rather than with the bandlimit machinery.

Both concrete forms let the nodes move freely (that only rewrites array contents). They differ in whether the count may change, which is the only thing that requires rebinding a field: FixedCountNodes is immutable, VariableCountNodes is not. make_plan's variable_npts keyword picks one, so mutability is opt-in rather than imposed.

FixedCountNodes <: AbstractNodeSet

Immutable node set — the default. The nodes may move anywhere via set_nodes!; only their count is fixed, because it sizes the point-indexed buffers.

VariableCountNodes <: AbstractNodeSet

Node set whose three point-sized fields are assignable, so set_nodes! also accepts a different number of points. The guru plans stay const — finufft_setpts! updates their point count in place. Request one with make_plan(…; variable_npts = true).

AbstractPlanTuning

How hard make_plan / make_spin_plan searches for the library settings that the problem does not fix: the FFTW planner thread count and effort for the sphere plans, and FINUFFT's nthreads and upsampfac. Subtype this and add _tune_sph / _tune_nufft methods to define a custom strategy.

NoTuning <: AbstractPlanTuning

Build with fixed, measured-good settings and time nothing — the default. Costs nothing beyond an ordinary plan build, and captures the large win of pinning the FFTW planner thread count instead of inheriting whatever the last foreign-library call left in that process global.

AutoTuning <: AbstractPlanTuning

Time the candidate FFTW planner thread counts and FINUFFT nthreads/upsampfac pairs under FFTW.ESTIMATE planning, keeping the fastest. Adds roughly 1.1x over NoTuning and costs O(100 ms) of trial builds, so it needs on the order of a thousand transforms on one plan to pay for itself. Outcomes are memoized per problem shape.

ThoroughTuning <: AbstractPlanTuning

AutoTuning plus a search over FFTW.MEASURE planner effort, which can beat ESTIMATE severalfold on the sphere synthesis/analysis step at large lmax but costs FFTW up to seconds to plan. For a plan that will live for a very long time.

nusht_type2!(f, C, plan)

Type 2 (synthesis): evaluate the field with spherical harmonic coefficients C at the M scattered points, writing values into f. Batched: C is (Nθ, Nφ) / (Nθ, Nφ, B) and f is length-M / (M, B).

Algorithm A = N·F·S: plan_sph2fourier to the bivariate Fourier series (S) → assemble the complex mode array (F) → FINUFFT type 2 (N).

nusht_type2!(fs, Cs, plans, backend = AutoBackend()) -> fs

Synthesize a collection of independent problems: for each i, nusht_type2!(fs[i], Cs[i], plans[i]). ThreadedBackend requires using OhMyThreads.

nusht_type1!(C, f, plan)

Type 1 (adjoint): given field values f at the M scattered points, apply A† and write the result to C. This is the exact Euclidean adjoint of nusht_type2! — the transpose, not an inverse — so A†A is symmetric positive definite and usable by an iterative solver. Use nusht_solve! to invert.

For exact analysis of a field already sampled on the Clenshaw-Curtis grid, use FastSphericalHarmonics.sph_transform: that is quadrature, needs no scattered-point machinery, and is a different operation from the adjoint.

nusht_type1!(Cs, fs, plans, backend = AutoBackend()) -> Cs

Adjoint analysis over a collection of independent problems; see nusht_type2!.

nusht_solve!(C, f, plan; ws=LSMRWorkspace(plan), maxiter=500, rtol=1e-6, conlim=0, verbose=false)

Exact inversion: solve min ‖A c − f‖ for coefficients C by LSMR on the Golub–Kahan bidiagonalization of A. Batched (B > 1) runs the columns as independent single-column solves.

Returns (C, iters, rel_res, converged) with rel_res = max_k ‖A†r_k‖/‖A†f_k‖ and converged = rel_res < rtol; ws.colres carries the same residual per column. rel_res is LSMR's own recurrence value for ‖A†r‖, floored at eps(T) since a relative residual is not resolvable below that — so an rtol under machine precision never reports convergence.

A column also stops when LSMR's condition estimate exceeds conlim (default 1/eps(T)), or when the bidiagonalization terminates exactly. That matters when the points do not determine the coefficients — M below (lmax+1)², or clustered so that they effectively do not — where A is rank deficient and the iteration has nothing left to resolve. converged == false is the signal that the point set, not the budget, was the limit.

nusht_solve!(Cs, fs, plans, backend = AutoBackend(); kwargs...) -> Cs

Exact inversion over a collection of independent problems; kwargs go to the single-problem method.

make_spin_plan([FE = ComplexF64,] θ_nodes, φ_nodes, lmax, s; tol=1e-10, ntrans=1, …)

Build a SpinNUSHTplan. Colatitudes θ ∈ [0,π], longitudes φ ∈ [0,2π).

FE is the field element type, positional as in make_plan; a spin field is complex in general, hence the default. A real FE asserts the field VALUES are real, which makes the mode array conjugate-symmetric and halves both the Δ-contraction and the NUFFT's θ axis — correct only if the coefficients satisfy the reality condition.

tuning (AbstractPlanTuning) and the nthreads / upsampfac overrides behave as in make_plan; there are no FastTransforms plans here, so only the FINUFFT settings are searched.

SpinNUSHTplan{T}

Plan for spin-s non-uniform spherical harmonic transforms at M scattered points up to degree lmax, transforming B co-located fields per call. Owns persistent FINUFFT guru plans (type 2 iflag=+1, type 1 iflag=−1) whose points are set to −θ once. The Wigner Δ^ℓ = d^ℓ(π/2) planes are generated on the fly by the Trapani–Navaza recurrence into two reused (2lmax+1)² buffers (dl_curr/dl_prev), so the plan is O(lmax²) memory, and the recurrence is numerically stable to ℓ ≈ 1024 (the explicit-factorial wigner_d sum it replaced loses all accuracy above ℓ ≈ 40).

Pass a WignerTable as wigner_table to precompute those planes instead: faster per transform, at O(lmax³) memory.

nusht_type2_spin!(f, sf, plan) -> f

Spin-weighted synthesis: evaluate the spin-s field with coefficients sf at the M scattered points, writing complex values into f. Batched: sf is (lmax+1, 2lmax+1[, B]), f is length-M / (M, B).

nusht_type2_spin!(fs, sfs, plans, backend = AutoBackend()) -> fs

Spin-weighted synthesis over a collection of independent problems. This path touches no FastTransforms state, so it carries none of the in-task hazard the scalar path works around.

nusht_type1_spin!(sf, f, plan) -> sf

Exact Euclidean adjoint of nusht_type2_spin!: scattered values f → spin-s coefficients sf.

nusht_type1_spin!(sfs, fs, plans, backend = AutoBackend()) -> sfs

Spin-weighted adjoint analysis over a collection; see nusht_type2_spin!.

nusht_solve_spin!(sf, f, plan; ws=LSMRWorkspace(plan), maxiter=500, rtol=1e-8, conlim=0, verbose=false)

Exact inversion of the spin-weighted synthesis at arbitrary scattered points: solve min ‖A sf − f‖ by LSMR on the Golub–Kahan bidiagonalization of A. Batched (B > 1) runs the columns as independent single-column solves. Same contract and return as nusht_solve!.

nusht_solve_spin!(sfs, fs, plans, backend = AutoBackend(); kwargs...) -> sfs

Spin-weighted exact inversion over a collection; see nusht_type2_spin!.

WignerTable(lmax, s; T = Float64)

Precomputed table of Q^ℓ_{m'm} = Δ^ℓ_{m'm} · Δ^ℓ_{m',−s} for every degree ℓ ≤ lmax — the only combination of the Wigner-d(π/2) planes that the spin contraction reads. It depends solely on (lmax, s, T), so one table is read-only and shareable: hand the same one to any number of plans, on any number of threads.

Without a table the whole O(lmax³) Trapani–Navaza sweep is regenerated on every transform — twice per solver iteration in nusht_solve_spin!, where the planes are identical every time. With one, that sweep is paid once and the contraction additionally runs ℓ-innermost, keeping each output column of G hot across the degree sum instead of re-streaming G once per degree.

The cost is O(lmax³) memory — roughly 2.8 MiB at lmax = 64, 22 MiB at 128, 175 MiB at 256 and 1.4 GiB at 512 in Float64. These are the same two modes s2fft exposes as "precompute" and "on the fly"; choose per problem rather than globally.

tbl = WignerTable(64, 1)
plans = [make_spin_plan(θs[i], φs[i], 64, 1; wigner_table = tbl) for i in eachindex(θs)]

sYlm(s, ℓ, m, θ, φ) -> Complex

Spin-weighted spherical harmonic ₛY_{ℓm}(θ,φ) = N_ℓ d^ℓ_{m,−s}(θ) e^{imφ} in the Goldberg / Newman–Penrose convention, with d the rotation matrix of wigner_d. Provided for direct evaluation.

spin_coeff_index(ℓ, m, lmax) -> CartesianIndex

Index into the dense (lmax+1, 2lmax+1) spin coefficient array for degree ℓ, order m.

wigner_d(ℓ, m, n, β)

Wigner small-d matrix element d^ℓ_{mn}(β) = ⟨ℓm|exp(-iβ Jy)|ℓn⟩ in the Condon–Shortley basis, via the explicit alternating sum with log-gamma factorials for stability.

Accuracy degrades above ℓ ≈ 40; use _wigner_d_halfpi_step! for β = π/2 at large ℓ.

_wigner_d_halfpi_step!(dl, dlp, ℓ, off) -> dl

Overwrite dl with the degree-ℓ Wigner-d plane at β = π/2, from the degree-(ℓ-1) plane dlp, via the Trapani–Navaza recurrence (Trapani & Navaza 2006; ssht ssht_dl.c; s2fft trapani.py). Both dl/dl_prev are (2lmax+1)², off = lmax+1. O(ℓ²) work, one previous plane — so a full sweep ℓ = 0…lmax is O(lmax³) work in O(lmax²) memory, numerically stable to ℓ ≈ 1024 (double).

Layout is n-major: dl[n+off, m+off] = d^ℓ_{m,n}(π/2) (McEwen–Wiaux/ssht convention, transposed). The recurrence runs downward in m, so storing n first makes every inner loop contiguous and independent — Stage A, Stage B and the n → −n fill all vectorize — and it is also the order the contraction reads, since Δ^ℓ_{m',m} = dl[m'+off, m+off]. So the relation to wigner_d is transposed: wigner_d(ℓ, m, n, π/2) == dl[n+off, m+off].

Only the eighth 0 ≤ n ≤ m ≤ ℓ is recurred; the rest is filled by the d(π/2) symmetries (transpose, m→−m, n→−n).

nusht_filter!(f_out, f_in, filter, plan; ws=LSMRWorkspace(plan), kwargs...)

Apply a spectral filter to f_in at the scattered points, writing to f_out (both length-M / (M, B)): fit coefficients with nusht_solve! → apply_transfer! (× H(ℓ)) → nusht_type2!. kwargs (rtol, maxiter, …) go to the solve. Uses the plan's coefficient scratch, and is allocation-free when a ws is supplied.

The fit is what makes this a filter. A H A† — the adjoint in place of a fit — is a smoothing operator, not A H A⁺: at scattered points A† is the transpose of the synthesis, not its inverse, and only on a quadrature grid do the two coincide. So filtering scattered data is iterative; hold a ws and reuse it across calls.

nusht_filter!(outs, ins, filter, plans, backend = AutoBackend()) -> outs

Spectral filtering over a collection of independent problems; see nusht_type2!.

nusht_filter_renorm!(f_out, mask, filter, plan; mask_filt=similar(f_out), ws, C_mask=nothing)

Renormalise the output of nusht_filter! to correct for land/ocean masking: divide by the filtered mask (the fraction of kernel weight over ocean). f_out must have been produced by nusht_filter!(f_out, f .* mask, filter, plan). Points where the filtered mask is below 0.01 are set to 0. Pass a reusable mask_filt scratch (shaped like f_out) to run allocation-free.

GaussianTransfer

Gaussian spectral filter: H(ℓ) = exp(-ℓ(ℓ+1) σ²/2)

where σ = scalem / Rm is the dimensionless filter width. In physical space this corresponds to convolution with a Gaussian-like kernel on the sphere. At degree ℓ the spatial scale is approximately R/ℓ.

Reference: Eq. (3) of Aluie et al. (2018), analogous to Gaussian in Fourier space.

gaussian_from_scale(scale_m, R_m=6.371e6)

Construct a GaussianTransfer from a physical filter scale in meters. Use this instead of GaussianTransfer(scale_m) to avoid ambiguity with the struct constructor which takes σ² directly.

TopHatTransfer

Ideal low-pass (sharp spectral cutoff) filter: H(ℓ) = 1 for ℓ ≤ L, else 0. Physical-space equivalent is a convolution with a zonal kernel whose Legendre spectrum is a boxcar. Note: this is spectrally sharp but spatially oscillatory (Gibbs phenomenon).

SharpSpectralTransfer

Alias for TopHatTransfer — identical sharp low-pass cutoff at degree L.

kernel_transfer(filter, ℓ) -> H

Evaluate the spectral transfer function H(ℓ) for a given filter type and degree ℓ. Returns a real scalar in [0, 1].

cutoff_degree(scale_m, R_m)

Compute the spherical harmonic cutoff degree L corresponding to a physical filter scale in meters on a sphere of radius R_m.

The relationship is L ≈ π * Rm / scalem, analogous to Nyquist for a circle of circumference 2π R_m.

nusht_type2(θs, φs, Cs, lmax, backend = AutoBackend(); tol, ntrans, tuning) -> fs

Synthesize N independent problems given their node sets: for each i a plan is built from (θs[i], φs[i]), Cs[i] is evaluated, and the field is returned as fs[i].

Use this instead of the plan-collection nusht_type2! when the backend is a DistributedBackend. When you already hold plans and are on a local backend, prefer the in-place form — it reuses them.

nusht_solve(θs, φs, fs, lmax, backend = AutoBackend(); tol, rtol, maxiter, …) -> Cs

Exact inversion of N independent problems given their node sets; see nusht_type2.

plot_field(θ, φ, f; colormap=:RdBu, markersize=8, title="", colorbarlabel="Field value") -> Figure

Scatter-plot a scalar field f sampled at scattered colatitude/longitude points (θ ∈ [0,π], φ ∈ [0,2π)), coloured by real(f) (so a complex/spin field plots its real part). Longitude on x, colatitude on y (poles top/bottom). Method supplied by NUFSHTCairoMakieExt — load it with using CairoMakie.

apply_transfer!(C, filter, lmax)

Multiply a FastSphericalHarmonics coefficient array C — size (lmax+1, 2lmax+1) or batched (lmax+1, 2lmax+1, B) — in-place by the transfer function H(ℓ) for each degree ℓ (broadcast across the batch dimension).

_nusht_true_adjoint!(C, f, plan, k = plan.B, kdfn = k)

The exact Euclidean adjoint of nusht_type2!: type-1 NUFFT, the transpose of the mode assembly, then P' per column. k bounds the columns the sphere loop visits and kdfn the NUFFT plan width.

_assemble_modes!(Z, G, lmax)

G (bivariate Fourier coefficients, sph2fourier layout, (lmax+1, 2lmax+1, B)) → Z, the (2lmax+1, 2lmax+1, B) complex mode array with Z[kθ+lmax+1, kφ+lmax+1] the coefficient of exp(i(kθ·θ + kφ·φ)). Exact: no grid, no doubling, no FFT.

_assemble_modes_adjoint!(G, Z, lmax)

Exact Euclidean adjoint of _assemble_modes!: each G[i,j] gathers conj(w)·Z over the same (up to four) mode entries the forward map scattered into. A real G takes the real part, which is the adjoint of the real→complex embedding; a complex one does not.