NUFSHT.jl

Non-Uniform Fast Spherical Harmonic Transforms — native Julia implementation of the Double Fourier Sphere (DFS) + NUFFT algorithm for spherical harmonic transforms at arbitrary scattered (colatitude, longitude) points.

Results

Synthesis + Round-Trip Accuracy (10⁻¹¹ relative error)

Synthesis and Accuracy

Inversion at Arbitrary Scattered Points

Inversion

Spectral Filtering

Ocean Mask + Renormalization

Mask Renorm

What it does

Given a field sampled at M arbitrary points on the sphere, NUFSHT.jl can:

Synthesise (Type 2): Evaluate a bandlimited field (given as SH coefficients) at any scattered point set in O(K log K + M) time.
Analyse (Type 1): Project scattered values back to SH coefficients; exact on the Clenshaw-Curtis (CC) quadrature grid.
Solve (LSMR): Exactly invert the synthesis operator at any scattered point set.
Filter: Apply isotropic spectral filters (Gaussian, top-hat, custom) entirely in harmonic space.

Quick start

using Pkg
Pkg.add(url="https://github.com/jbphyswx/NUFSHT.jl")

using NUFSHT, FastSphericalHarmonics

lmax = 30
θ = rand(5000) .* π       # colatitudes ∈ (0,π)
φ = rand(5000) .* 2π      # longitudes ∈ [0,2π)
plan = make_plan(θ, φ, lmax; tol=1e-8)

# Synthesise at scattered points
C = zeros(lmax+1, 2lmax+1)
C[sph_mode(2, 0)] = 1.0
f = zeros(length(θ))
nusht_type2!(f, C, plan)

# Exact inversion (for non-CC scattered points)
C_rec = similar(plan.C)
C_rec, iters, rel_res = nusht_solve!(C_rec, f, plan; rtol=1e-6)

Which function should I use?

Scenario	Function
Evaluate SH expansion at scattered points	`nusht_type2!`
Apply the adjoint `A†`	`nusht_type1!`
Invert at arbitrary scattered points	`nusht_solve!`
Apply spectral filter at scattered points	`nusht_filter!`
Filter with land/ocean mask	`nusht_filter!` + `nusht_filter_renorm!`

Module

NUFSHT.NUFSHT — Module

NUFSHT.jl — Non-Uniform Fast Spherical Harmonic Transform (native Julia)

Double Fourier Sphere (DFS) + nuFFT spherical harmonic transforms at arbitrary scattered (colatitude, longitude) points. The synthesis operator factors as A = N·F·S:

S (plan_sph2fourier): the Legendre step, taking SH coefficients to the Double-Fourier-Sphere bivariate Fourier series. No equiangular grid is formed.
F (_assemble_modes!): the cos/sin bivariate Fourier basis S produces, rewritten as the complex exponentials a NUFFT evaluates. Both mode axes carry wavenumbers -lmax…lmax.
N (FINUFFT guru type 2 / type 1): non-uniform FFT evaluating the 2D Fourier series at the scattered points.

A NUSHTplan owns persistent FINUFFT guru plans (built once, points set once) and every work buffer, so repeated transforms — filtering, or the hundreds of matvecs in nusht_solve! — allocate nothing and never re-plan. All calls transform a batch of B = plan.B co-located fields (ntrans); B = 1 methods accept plain vectors/matrices.

References

Merilees (1973); Townsend & Olver (2015); Reinecke & Seljebotn (2013, A&A 554 A112); Keiner, Kunis & Potts (2009); Belkner et al. (2024, arXiv:2406.14542).
FastSphericalHarmonics.jl, FINUFFT.jl, FastTransforms.jl.

source