Internals

_backend_nthreads(exec) -> Int

Thread count for a library plan built under exec: Threads.nthreads() on a threaded backend, 1 otherwise.

Every provider is handed this explicit count, so the thread budget follows the -t Julia was started with. Several of them (FINUFFT, and NUFSHT's inner FINUFFT) read 0 as "every core on the machine", which is a different quantity.

_to_host(a) -> Array

a as a host Array, returning it unchanged when it already is one. copyto! covers every KA / GPUArrays device array, so this is device-generic.

DirectSumCartesianPlan{FT}

Reusable Cartesian direct-sum plan: the per-axis dense DFT matrices, the working arrays the contraction walks through, the wavenumber axes, the Nyquist-twin storage, and the quadrature-scaled field buffer, all built once for a fixed grid, ms, iflag and batch shape.

The matrices depend only on the grid and ms, and every buffer only on the shapes, so a time loop rebuilds none of them. Execute with calculate_spectrum!(coeffs, plan, field).

DirectSumSphericalPlan{FT}

Reusable spherical direct-sum plan: the grid's materialized nodes (or its ring table, or its longitude DFT matrix and latitude quadrature — whichever its layout reads), the Legendre recurrence tables, and the (lmax+1, nrings, B) longitude buffer, all built once for a fixed grid, lmax and batch shape. A complex field additionally needs the exp(+imλ) longitude transform, so the plan holds its matrix and buffer too; the element type passed to plan_spectrum is what decides.

Building these is a minority of a transform's time and about three quarters of its allocations (the latitude quadrature is a Gauss–Legendre root solve, and materialize is documented as an explicit allocating call), so reuse across a time loop removes most of the heap traffic. Execute with calculate_spectrum!(coeffs, plan, field).

DirectSumSynthesisPlan{FT}

Reusable direct-sum inverse: the grid, resolution, sign and batch shape the inverse is fixed by, so synthesize! writes into the caller's array and allocates nothing of its own.

The Cartesian inverse reads the grid's axes per call, and the spherical one its nodes and Legendre tables; a Nyquist twin arrives with the coefficients on ks, never held here.

HybridPlan{T,D,R}

Reusable plan for a Cartesian grid uniform in some directions and stretched in others: the FFTW plan over the uniform axes, one 1-D type-1 NUFFT per stretched axis (each with its own strengths/spectrum buffers), and the D - length(sdims) + 1 working arrays the passes write through.

R marks a real field, so the publish branch folds at compile time. Execute with calculate_spectrum!(coeffs, plan, field).

_hybrid_plan(nufft, exec, grid, ::Type{T}, ms, umask; batch, iflag, eps, kwargs...)

Build the reusable hybrid composite. The FFTW extension's plan_spectrum calls this where the grid is uniform in some directions and stretched in others, the same split its forward makes.

_hybrid_derive(grid, ::Type{Tr}, ms, umask, R, iflag, eps, batch) -> NamedTuple

Everything the hybrid composite needs before any transform runs: the per-axis coordinates, origins and Fourier lengths, which dims the FFT takes and which the NUFFT takes, whether a Nyquist twin is required (and so whether the FFT pass halves axis 1), the publish extents and phase, and the wavenumber axes with their twins.

The one-shot and the reusable plan both read this, so the two cannot disagree about the derivation.

cart_setup(layout, grid, T, ms, batch, iflag) -> NamedTuple

Everything a Cartesian direct sum reads from grid: the wavenumber axes, and for a tensor layout the per-axis DFT matrices and the working arrays their contraction walks through, plus the Nyquist-twin storage. T is the field's element type, whose realness fixes the packed layout.

Run it against a field with cart_run!. The result depends on the grid, ms, iflag and the batch shape alone, so DirectSumCartesianPlan holds it across executions.

cart_run!(coeffs, layout, setup, field) -> ks

Fill coeffs with the Cartesian spectrum of field against a cart_setup result, and return the wavenumber axes, the halved axis carrying a Nyquist twin where one is required.

field arrives already scaled by the grid's quadrature factor. A tensor layout contracts axis by axis, a point cloud sums directly, and both write only into coeffs and the setup's own buffers.

AxisPass{FT, ND, M, A}

One axis's contraction: its dense DFT matrix, the two matrix aliases the matmul reads and writes, and — for an axis past the first — the staging buffers and the permutation bringing it to the front and back.

Every alias and permutation is built once, so a pass reshapes nothing and inverts no permutation per call. Run it with _axis_pass!.

sph_setup(layout, grid, T, lmax, B; sampling, weights) -> NamedTuple

Everything a spherical direct sum reads from grid: its materialized nodes, or its ring table, or its longitude DFT matrix and latitude quadrature — whichever its layout uses — plus the Legendre recurrence tables and the longitude buffer. T is the coefficient element type; a complex one adds the exp(+imλ) partner transform.

Run it against a field with sph_run!. DirectSumSphericalPlan holds it across executions.

sph_run!(coeffs, layout, setup, field, lmax, B) -> ks

Fill coeffs (lmax+1, 2lmax+1, B) with the spherical spectrum of field against a sph_setup result, and return (0:lmax, -lmax:lmax).

A tensor layout runs one longitude DFT for every latitude, a ring layout one per ring, and a point cloud projects per node; all three then contract over θ against the Legendre tables.

sph_synth_setup(grid, ::Type{FT}, lmax) -> NamedTuple

Everything the spherical inverse reads from grid: its nodes, the Legendre recurrence tables and the per-node Plm scratch. Synthesis carries no quadrature, so this is the node list only.

Run it against coefficients with sph_synth_run!; DirectSumSynthesisPlan holds it.

sph_synth_run!(out, setup, coeffs, lmax) -> out

Evaluate the spherical harmonic sum at the setup's nodes, writing out (N, batch…).

physical_wavenumbers(Ls, ms, ::Val{real}) -> NTuple
physical_wavenumbers(grid, ms, ::Val{real}) -> NTuple

Packed native-order wavenumber axes matching the packed coefficient layout: for a real transform axis 1 is a nonnegative rfft axis (Packing.RFFTAxis, length ms[1]÷2+1) and axes 2:D are full fft axes (Packing.FFTAxis); for a complex transform every axis is a full fft axis. ms[d] is the full transform length along d.

PointwiseCartesian

Cartesian grids whose coordinates are stored per point: a node cloud, and a curvilinear grid, whose coordinate arrays hold one value per cell.

Both index linearly in the field's own order, so a transform reads them identically — the direct sum, its Nyquist twin, both inverses, and a NUFFT provider's point list. A tensor grid stores axes instead and has its own factorized methods, which are more specific than this.

SphericalHarmonicGeometry

Geometries whose points are addressed by a DIRECTION, so a spherical-harmonic transform applies to a field sampled on them: a sphere, and a spheroid.

A spherical harmonic is orthonormal over directions, and a spheroid node's direction is its geocentric one, which _colatitude reads through the geometry's own embedding. What such a transform returns is the SURFACE expansion of the field: the spheroid's radius varies with latitude, and that radial dependence belongs to the solid-harmonic expansion, which is a different transform.

field_batch_shape(grid, field) -> Tuple

Trailing batch sizes of field on grid: the dims after the leading ndims(grid) spatial ones, per the (spatial…, batch…) contract Problem states. ndims(grid) is 1 for a node cloud and N for a grid that indexes as an N-D array, so a path serving both reads its batch from here.

point_coordinates(::Type{FT}, grid, D) -> (coords::NTuple{D}, spatial::Tuple)

The D per-point coordinate vectors of grid in the field's own column-major order, and the spatial shape a field on it carries.

FlowGeometries.Grids.materialize answers this for every architecture: a node cloud returns its nodes, a tensor grid the expansion of its axes, a curvilinear grid its per-cell values, and a pixelization its pixel centres. It allocates D·n numbers on a tensor grid (its docstring says so), so a repeated caller holds the result in a plan.

axis_geometry(::Type{FT}, grid, D) -> (offsets, ranges)

Per-axis origin and Fourier length: offsets[d] is the smallest coordinate along direction d, and ranges[d] its wrap length where the direction wraps and 1 where it does not.

offsets comes from FlowGeometries.Grids.bounds, which is O(1) on a structured grid (an axis is monotone, so its extremes are its endpoints) against an O(N) scan of the coordinates.

ranges reads isperiodic before period, whose own contract states it is meaningful only where the direction wraps: a grid may declare a direction non-periodic while carrying a nonzero period entry, so the flag decides. A non-periodic direction gets 1, making its wavenumbers raw per-sample.

axis_range(FT, grid, d) — direction d's Fourier length: its wrap period, or 1 where it does not wrap.

_sph_points(grid) -> (θ, φ)

Per-point colatitude/longitude lists for any spherical grid, in the transform's (θ = colatitude, φ = longitude) convention.

The nodes come from FlowGeometries.Grids.materialize, which every grid architecture answers: on an unstructured grid it returns the nodes themselves, on a structured (nlon, nlat) grid the expansion in column-major order (longitude fastest, matching size and the field layout), and on a pixelization (HEALPix, cubed-sphere, icosahedral, ring, Yin-Yang) the pixel centres. Those pixelizations answer neither coordinates(grid, d) nor sampling(grid), so materialize is the accessor that spans them.

_colatitude(geometry, λ, φ) -> FT

Colatitude of the direction the node (λ, φ) lies in.

On a sphere that is π/2 − φ. On a spheroid the stored φ is the geodetic latitude, whose difference from the geocentric one reaches ~0.19° at mid-latitudes on Earth's ellipsoid, so the node is embedded through the geometry's own geodetic_to_cartesian and the direction read back from it. The result is independent of λ on both (a spheroid is a surface of revolution), so a latitude row remains one iso-latitude ring and the ring-factorized transform still applies.

Tensor-product (nlon, nlat) sphere: one longitude axis shared by every latitude.

Iso-latitude rings whose longitude count varies by ring (HEALPix, reduced Gaussian).

No iso-latitude structure: a point cloud (Fibonacci, cubed-sphere, icosahedral, Yin-Yang).

_sph_sampling(grid) -> AbstractSphericalSampling or nothing

The grid's own node-set recipe. A structured grid records one and carries it in its TYPE, so the traits read from it fold at compile time. A HEALPixGrid answers from its resolution parameter, nside fixing its recipe entirely. Anything else records none.

_sph_layout(grid) -> TensorSphere | RingSphere | ScatteredSphere

Which spherical algorithm grid admits, taken from its sampling's declared traits.

_ring_table(grid, ::Type{FT}; lmax=nothing) -> (; ranges, θ, w, λ)

Per iso-latitude ring: its slice of the flattened point vector, its colatitude, and the per-point quadrature weight inside it; plus every point's longitude.

Colatitude and the weight are read at the ring's first point, both being constant along a ring. The weight comes from Grids.measure, rescaled so the weights total 4π — the transform's convention, the same total _sph_node_weights carries. A ring layout's measure is a RingwiseVector and an equal-area one's an Axes.ConstantVector, whose sum is O(nrings) and O(1), so no dense ∏N-sized measure is built.

_sht_weights(grid, nlat; sampling=nothing, weights=nothing) -> Vector

Per-colatitude quadrature weights for a structured spherical grid, from its own node-set recipe FlowGeometries.Grids.sampling(grid) via SphericalSampling.latitude_weights (Σw = 2, carrying the sinθ Jacobian — the transform's convention). weights or sampling override it.

The recipe is the only thing that determines the quadrature: Gauss–Legendre and an arbitrary lat–lon grid are the same coordinates to round-off, and only the recipe says which rule is exact on them. A grid built from raw axes records none, so it raises here and the caller states the quadrature.

_sph_node_weights(grid, ::Type{FT}, N, weights) -> Vector{FT}

Per-node quadrature weights for a scattered spherical grid, from its own per-node measure rescaled to sum to 4π — the transform's convention, the same total the structured path's latitude_weights ⊗ 2π/nlon carries. weights overrides it. Only the relative measures matter, so a grid of true dual-cell areas integrates exactly and a grid of equal-area nodes gives 4π/N.

quadrature_scale(grid, ::Type{FT}, N) -> Union{Nothing, AbstractVector}

Per-node factor α_j = N·w_j/Σw taking a node-count-normalized transform to the measure-weighted quadrature over grid, or nothing where the measure is constant and α_j is exactly one. N is the number of spatial samples the field carries. Flat and of length N, indexed by the same column-major point order the field and the coordinates use, so it reads against a strengths vector and stages to a device as one.

Read from Grids.measure, never measure_array: a rectilinear grid's measure is a SeparableMeasure and a ring layout's a RingwiseVector, whose sum and extrema cost O(Σ Nᵈ) and O(nrings) against a dense O(∏ Nᵈ), and whose scaling by a constant stays in that representation. The result is therefore lazy on such a grid, carrying O(Σ Nᵈ) numbers for N entries.

quadrature_weighted(grid, field) -> field

field scaled by the grid's per-node quadrature factor (see Grids.quadrature_scale), so a backend's node-count normalization becomes the measure-weighted quadrature over grid. Returns field itself where the measure is constant, which is every uniform grid and every scattered grid whose nodes carry equal weight. The scaled copy is what both the coefficients and the Nyquist twin are built from, so they stay consistent.

_warn_bandlimit(sampling, nlat, lmax)

Warn when degree lmax exceeds what nlat latitude samples of sampling support, so the coefficients above it are not resolved by that grid's quadrature.

Warned, never raised: an over-resolved request still returns coefficients, and a caller may want them. The relation is the sampling's own SphericalSampling.bandlimit — notably nlat÷2 − 1 for Driscoll–Healy against nlat − 1 for Gauss–Legendre and Clenshaw–Curtis. maxlog = 1 caps this per logger instance, so a fresh logger still observes it.

_zeroed_inactive(field, grid) -> AbstractArray

field with its inactive cells set to zero, or field itself where every cell is active.

A transform that consumes complete longitude rows — a library SHT on its own grid, or a longitude rfft — is linear in the field, so zeroing an inactive cell's datum removes it from the sum exactly, matching what the per-point paths get by skipping it. The weights stay the sampling's own, so the active cells still carry their own solid angles and total the covered one. Writing zero also clears a NaN, which masked data commonly holds.

covered_area(grid) -> T

Total measure of the cells that participate: the extent a transform over grid integrates over.

Equal to sum(measure(grid)) on an unmasked grid — 4πR² for a layout tiling the whole sphere, which a GridMeasure answers without visiting a cell — and the active cells' share of it under a mask.

sky_fraction(grid) -> T

Fraction of grid's cells' total measure that participates: covered_area(grid) / sum(measure(grid)). One on an unmasked grid. On a partial-sky spherical grid this is f_sky, so a caller wanting the full-sphere-normalized convention scales coefficients by inv(sky_fraction(grid)).

_sht_applicable(grid, ms) -> Bool

Whether the FastSphericalHarmonics analysis applies to grid at ms. Its transform is addressed by array position on FastTransforms' own Clenshaw–Curtis nodes, so only that grid at the matching size qualifies; the extension owns the check and answers false until it loads.

hermitian_request_size(ms) -> NTuple

Mode counts a transform must produce so that every mode of a real field's packed half is present. Axis 1 asks for one extra at even N₁, which makes its length odd and so puts both ±N₁/2 on the native axis; the packed half's axis-1 frequencies 0:N₁÷2 are then the leading N₁÷2+1 entries for even and odd N₁ alike. Negating a packed entry whose axis 1 sits at +N₁/2 and whose axis d sits at −N_d/2 lands on +N_d/2, off-axis, so a conjugate read cannot substitute.

offset_phase(Tr, ms, offsets, ranges, M, ::Val{R}, iflag=1) -> Array

Grid-offset correction exp(-iflag·i·k·x₀)/M over the packed half (R = true) or the full native spectrum (R = false), for points scaled from a domain starting at offsets with period ranges.

publish_packed!(coeffs, fk, phase, ns, pms, ntrans, neg) -> coeffs

Write a real field's packed half from the full native spectrum fk of requested size ns (see hermitian_request_size). The half is the leading pms[1] entries of axis 1, so each axis-1 run copies contiguously out of the longer ns[1] run and the axes above it match one for one. phase is the packed-layout offset correction; neg conjugates the published values.

This is the HOST write: a linear scalar sweep, allocation-free on a host array. A device path takes packed_half_view and broadcasts, which reads no element from the host.

packed_half_view(fk, ns::NTuple{D,Int}, pms::NTuple{D,Int}, ntrans) -> SubArray

The packed half of a full native spectrum as a lazy view: the leading pms[1] entries of axis 1 with every other axis whole, batch axis last.

A device path publishes with out .= view .* phase (or its conjugate); publish_packed! is the host form of the same write.

NyquistTwin(slices)

Coefficients at the modes a packed half cannot reach by index negation. Negating a mode sends an even axis sitting at −N_d/2 to +N_d/2, which is off the native axis; the two agree only when the grid samples periodically. slices[mask] covers the modes whose axes-at-−N_d/2 are exactly mask (bit d−2 per axis d ≥ 2): entry I holds C at k₁ kept, +N_d/2 on the masked axes, and −k_e elsewhere. Each slice collapses its masked axes to extent 1, so the whole object holds a hyperplane's worth of values. Index it as t[mask, I] with mask = nyquist_mask(ks, I).

twin_at(t, mask, I, b)

Entry of the mask slice at spectral index I and batch member b. getindex(t, mask, I) takes an index whose rank already spans the batch dims; this form takes a purely spectral I of rank D and the batch member separately, for a caller that walks the spectral dims and offsets the batch linearly.

twin_table(Tr, ms, mask, dims, offsets, ranges, M; phis, normfactor, conjugate) -> (shape, src, fac)

Gather table for the mask slice of a NyquistTwin read out of a spectrum with array dimensions dims: twin[i] = S(fk[src[i]]) * fac[i], where S is conj when conjugate and the identity otherwise. Each entry's target frequency is +N_d/2 on the masked axes, k₁ on axis 1, and −k_d elsewhere; fac carries the grid-offset phase at that target frequency together with 1/M and, when phis is given, normfactor / ∏ phis[d].

src indexes the target frequency itself (conjugate = false, for a spectrum that holds +N_d/2 outright, such as an oversampled NUFFT grid) or its negative (conjugate = true, for a Hermitian spectrum of a real field, where fk[-k] = conj(fk[k]) puts the value within a native mode set). Both indices come from the target frequency, so an axis already sitting at −N_d/2 still reads its +N_d/2 partner. phis holds one per-axis coefficient table sampled at the native wavenumbers, even in k.

twin_slice_shape(ms, mask) -> NTuple

Shape of the mask slice of a NyquistTwin: the masked axes collapse to their single +N_d/2 mode, and an odd axis has none, which empties the slice.

ConjTwinSlice(slice, src, fac)

One mask's NyquistTwin slice together with the table that fills it by conjugate reads of a Hermitian native spectrum (see twin_table with conjugate = true).

conj_twins(Tr, ks_phys, ms, ns, offsets, ranges, M, batch) -> (ks, slices)

ks_phys with a NyquistTwin attached to its halved axis, plus the per-mask ConjTwinSlices that fill it from a Hermitian native spectrum of requested size ns. A real field's transform is Hermitian, so each twin value is a conjugate read at the negated frequency, which that spectrum holds.

gather_conj_twins!(slices, fk, Pm, ntrans, neg)

Refill every ConjTwinSlice from the native spectrum fk, whose transforms are Pm modes apart. Entry i of transform t is conj(fk[(t-1)Pm + src[i]]) * fac[i], conjugated again where neg.

axis_layout(sz, d) -> (pre, N_d, post)

sz viewed as (pre, N_d, post) for a transform along axis d: pre is the product of the axes below d, post the product of those above. The lines to transform number pre · post, and line r starts at mod1(r, pre) + ((r - 1) ÷ pre) · pre · N_d.

axis_work_shape(Ns, ns, d, batch) -> NTuple

Shape a separable transform's working array carries entering pass d: the axes below d already hold their ns modes, the rest still hold their Ns grid points, and the batch rides above. d = 1 is the field's own shape and d = D+1 the finished spectrum's.

axis_chunk_offsets!(inoff, outoff, base, nvalid, pre, Nd, m)

Starting linear indices of lines base+1 … base+nvalid in the input (Nd per line) and in the output (m per line). Lines sharing a slab get consecutive offsets, so the gather reads contiguously.

gather_axis_block!(cjs, A, inoff, nvalid, pre, Nd)

Fill cjs[t][i] = A[inoff[t] + (i-1)·pre] for t ≤ nvalid, zeroing the strength vectors past nvalid so a partial final chunk still presents the plan's full transform count.

scatter_axis_block!(out, fks, outoff, nvalid, pre, m)

Write out[outoff[t] + (j-1)·pre] = fks[t][j] for t ≤ nvalid, the inverse of gather_axis_block! on the transformed axis length m.

spatial_shape(prob) — the leading spatial array sizes (N_1, …).

batch_shape(prob) — the trailing batch sizes (() if none).

n_batch(prob) — number of trailing batch dims.

batch_length(prob) — total batch slices, prod(batch) (1 if no batch axes).

coeff_output_size(spectral::Tuple, prob) -> NTuple

Shape of the coefficient array: (spectral…, batch…).

stack_fields(fields::Tuple) -> AbstractArray

Convenience for the multi-field call form: stack equal-shaped field arrays along a new trailing batch axis ((spatial…, batch…, NU)). This materializes a combined array — pass a single (spatial…, batch…) array to avoid the copy.

window_function(win::AbstractWindow, n::Integer, ::Type{T}=Float64) -> Vector{T}

Allocate and return the length-n taper win.

window_function!(w::AbstractVector, win::AbstractWindow) -> w

Fill w (length n) in place with the taper win.

window_correction(w::AbstractVector) -> (S1, S2)

Coherent-gain factor S1 = (Σ w)/n and power factor S2 = (Σ w²)/n. Amplitude spectra divide by S1; power/energy spectra divide by S2 to preserve variance.

detrend!(x::AbstractVector, d::AbstractDetrend) -> x

Detrend x in place according to d.

is_identity(spec) — whether spec leaves a field untouched, so the transform reads it directly.

axis_taper(::Type{FT}, window, n) -> Vector{FT}

Length-n taper for one axis, scaled to unit mean square (Σw²/n == 1).

Under that scaling the tapered field carries the variance of the original, so Parseval holds on the coefficients with no correction applied afterwards: with C_w(k) = (1/N) Σ w f e^{-ikx}, the folded sum Σ_k |C_w|² is mean|w·f|², which returns mean|f|² exactly when Σw²/n = 1 on every axis (the tensor product's mean square is the product of the axes'). NoWindow already satisfies it, so a rectangular window is an exact no-op.

detrend_spatial!(field, detrend, nsp::Int)

Detrend the leading nsp spatial dims of field in place, per batch slice.

Demean subtracts each slice's spatial mean, which needs no axes and so applies to any grid. LinearDetrend subtracts the least-squares linear trend along each spatial axis in turn — a separable plane removal — and so reads an axis order; a caller whose grid has none gets an ArgumentError from the entry point.

apply_window!(field, tapers::Tuple, nsp::Int)

Multiply the leading nsp spatial dims of field by the tensor product of tapers, in place, over every batch slice. One taper per spatial axis, each already scaled by axis_taper.

sided_factor(s::AbstractSidedness, k, kmax) -> Real

Folding multiplier. TwoSided → 1 everywhere. OneSided → 2 for interior wavenumbers, 1 at DC (k≈0) and Nyquist (k≈kmax) which have no negative-frequency partner.

LegendreTables{FT}

Precomputed, point-independent recurrence coefficients for the normalized associated Legendre functions $\bar P_\ell^m$ up to degree lmax. Built once per transform via legendre_tables; the per-point table is then filled in O(lmax²) with no sqrt calls by fill_legendre!.

legendre_tables(::Type{FT}, lmax::Int) -> LegendreTables{FT}

Precompute the recurrence coefficients up to degree lmax.

fill_legendre!(Plm::AbstractMatrix, t::LegendreTables, x, s, lmax)

Fill Plm[l+1, m+1] = \bar P_\ell^m(x) for m = 0:lmax, l = m:lmax at a single point with x = cosθ, s = sinθ, reusing the precomputed coefficients in t. Entries with l < m are left untouched (the projection never reads them).

normalized_legendre(l, m, x, s) -> FT

Single normalized associated Legendre value $\bar P_\ell^m(x)$ (m ≥ 0) computed by on-the-fly recurrence. Reference implementation used for validation; the hot path uses fill_legendre!.