healpix_analyse.healpix_sht#

Ring-based Spherical Harmonic Transform (SHT) for full-sky HEALPix maps.

Algorithm#

Adapts the FOSCAT ring-FFT algorithm to the healpix-analyse package. By exploiting the HEALPix ring structure, the number of Legendre summations is reduced from 12·nside² (all-pixel, alm_latlon.py) to 4·nside-1 (one per iso-latitude ring), giving a ~3·nside speed-up at fixed lmax.

For each iso-latitude ring r with N_r uniformly-spaced pixels at colatitude θ_r and starting longitude φ_0^r:

  1. F_r(m) = FFT_m(f_ring_r) × exp(-i·m·φ_0^r) [phase-corrected ring FFT]

  2. a_lm = sqrt(2l+1)/N × Σ_r √4π·P̃_lm(cos θ_r) × F_r(m)

where N = 12·nside². This follows the standard orthonormal convention (same as healpy and alm_latlon.py):

Y_lm(θ,φ) = sqrt((2l+1)/(4π)) · P̃_lm(cos θ) · exp(i·m·φ) a_lm = ∫ f(θ,φ) · Y_lm*(θ,φ) dΩ

Differentiability#

All hot-path operations use torch.fft.fft/ifft and torch.einsum, which are supported by torch.autograd. Geometry (Legendre tables, phase matrices, permutation indices) is precomputed once and stored as non-grad constants.

Spin support (requires quaternionic + spherical)#

spin=0 : scalar (temperature) field spin=1 : vector field (curl/divergence decomposition) spin=2 : CMB polarisation Q/U → E/B modes

The relation between (Q, U) and the spin-weighted harmonics sY_lm is:

(Q ± iU)(n) = Σ_lm [±almE ∓ i·almB] · ±sY_lm(n)

which gives:

almE = -(Σ+ + Σ-) / 2 almB = (Σ+ - Σ-) / (2i)

where Σ+(m) = Σ_r +sY_lm*(θ_r)·F_r^+(m), F_r^+ = FFT(Q + iU), and Σ-(m) = Σ_r -sY_lm*(θ_r)·F_r^-(m), F_r^- = FFT(Q - iU).

Public API#

HEALPixSHT(level, lmax, dtype, device, ellipsoid)

.map2alm(im, nest=False) → Tensor […, K] .alm2map(alm, nest=False) → Tensor […, N] .map2alm_spin(Q, U, spin, nest) → (almE, almB) each […, K] .alm2map_spin(almE, almB, spin, nest) → (Q, U) each […, N] .anafast(im, map2, spin, nest) → Tensor […, lmax+1]

Flat alm layout (same as healpy / alm_latlon.py):

[m=0: l=0..lmax | m=1: l=1..lmax | … | m=lmax: l=lmax] K = (lmax+1)·(lmax+2)//2

Attributes#

ArrayLike

Classes#

HEALPixSHT

Ring-based Spherical Harmonic Transform for full-sky HEALPix maps.

Module Contents#

healpix_analyse.healpix_sht.ArrayLike#
class healpix_analyse.healpix_sht.HEALPixSHT(level: int, lmax: int | None = None, dtype: torch.dtype = torch.float32, device: str | torch.device | None = None, ellipsoid: str = 'sphere')[source]#

Ring-based Spherical Harmonic Transform for full-sky HEALPix maps.

Parameters:

  • level (int            Grid4Earth/HEALPix level (integer >= 0).) – Internally, nside = 2**level.

  • lmax (int or None    Maximum multipole.  Default: 3·2**level - 1.)

  • dtype (torch.dtype    torch.float32 or torch.float64.)

  • device (str or device  Defaults to CUDA when available.)

  • ellipsoid (str "sphere" (default) or "WGS84")

level#
nside#
lmax#
dtype = Ellipsis#
cdtype = Ellipsis#
ellipsoid = ''#
device#
n_pix#
n_rings#
n_alm#
map2alm(im: ArrayLike, nest: bool = False) → torch.Tensor[source]#

Spin-0 analysis: HEALPix map → spherical harmonic coefficients.

a_lm = ∫ f(θ,φ) · Y_lm*(θ,φ) dΩ

Fully differentiable with respect to im.

Parameters:

  • im ((..., N) real array-like / Tensor.  RING order by default.)

  • nest (bool  Set True for NESTED input.)

Returns:

alm ((..., K) complex Tensor.) – K = (lmax+1)·(lmax+2)//2. Flat layout: [m=0: l=0..lmax | m=1: l=1..lmax | … | m=lmax].

alm2map(alm: ArrayLike, nest: bool = False) → torch.Tensor[source]#

Spin-0 synthesis: a_lm → HEALPix map.

f(θ,φ) = Σ_{lm} a_lm · Y_lm(θ,φ)

Fully differentiable with respect to alm.

Parameters:

  • alm ((..., K) complex array-like / Tensor.)

  • nest (bool  Set True for NESTED output.)

Returns:

im ((..., N) real Tensor.)

map2alm_spin(Q: ArrayLike, U: ArrayLike, spin: int = 2, nest: bool = False) → Tuple[torch.Tensor, torch.Tensor][source]#

Spin-s analysis: (Q, U) → E-mode and B-mode alm tensors.

For spin=2:

Q = +P cos(2χ), U = +P sin(2χ) (IAU CMB convention)

For spin=1:

Q = longitudinal (div) component, U = transverse (curl) component

For spin=0:

Delegates to map2alm(Q) and map2alm(U) independently.

The E/B decomposition follows:

Σ+(m) = Σ_r +sY_lm*(θ_r) · F_r^+(m), F^+ = FFT(Q + iU) Σ-(m) = Σ_r -sY_lm*(θ_r) · F_r^-(m), F^- = FFT(Q - iU) almE = -(Σ+ + Σ-) / 2 almB = (Σ+ - Σ-) / (2i)

Fully differentiable with respect to Q and U.

Parameters:

  • Q, U ((..., N) real Tensors.  RING order by default.)

  • spin (int  Spin weight (1 or 2; use 0 for scalar pair).)

  • nest (bool  Set True for NESTED input.)

Returns:

almE, almB ((..., K) complex Tensors.)

alm2map_spin(almE: ArrayLike, almB: ArrayLike, spin: int = 2, nest: bool = False) → Tuple[torch.Tensor, torch.Tensor][source]#

Spin-s synthesis: (almE, almB) → (Q, U) maps.

This is the adjoint of map2alm_spin.

For spin=2, the reconstructed maps are:

Q(n) = Re[ Σ_{lm} (-almE + i·almB) · +2Y_lm(n) ] U(n) = Im[ Σ_{lm} (-almE + i·almB) · +2Y_lm(n) ]

For spin=0: delegates to alm2map(almE) and alm2map(almB).

Fully differentiable with respect to almE and almB.

Parameters:

  • almE, almB ((..., K) complex Tensors.)

  • spin (int  Spin weight.)

  • nest (bool  Set True for NESTED output.)

Returns:

Q, U ((..., N) real Tensors.)

anafast(im: ArrayLike, map2: ArrayLike | None = None, spin: int = 0, nest: bool = False) → torch.Tensor[source]#

Angular power spectrum C_l (or cross-spectrum C_l^{12}).

For spin=0 (scalar):

C_l = 1/(2l+1) × [ |a_l0|² + 2 Σ_{m=1}^{l} |a_lm|² ]

For spin>0 (E/B modes):

C_l^{EE}, C_l^{BB}, C_l^{EB} are returned as a (3, lmax+1) tensor.

Fully differentiable with respect to im (and map2).

Parameters:

  • im ((..., N) or (..., 2, N) array-like.) – Scalar map, or [Q, U] stacked on axis -2 for spin > 0.

  • map2 (same shape as im or None.  If given, compute cross-spectrum.)

  • spin (int  0 for scalar, 1 or 2 for polarisation.)

  • nest (bool)

Returns:

cl ((..., lmax+1) float Tensor for spin=0.) – (…, 3, lmax+1) float Tensor for spin>0 — [EE, BB, EB].

uv_to_curl_div(u: ArrayLike, v: ArrayLike, nest: bool = False) → Tuple[torch.Tensor, torch.Tensor][source]#

Decompose a tangent-plane vector field (u, v) into curl and divergence.

Uses spin=1 spherical harmonics:

almE (divergence-free part) → div map almB (curl-free part) → curl map

Parameters:

u, v ((..., N) real arrays — east and north components.)

Returns:

div, curl ((..., N) real Tensors.)