healpix_analyse.healpix_sht =========================== .. py:module:: healpix_analyse.healpix_sht .. autoapi-nested-parse:: healpix_sht.py ============== 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 ---------- .. autoapisummary:: healpix_analyse.healpix_sht.ArrayLike Classes ------- .. autoapisummary:: healpix_analyse.healpix_sht.HEALPixSHT Module Contents --------------- .. py:data:: ArrayLike .. py:class:: HEALPixSHT(level: int, lmax: Optional[int] = None, dtype: torch.dtype = torch.float32, device: Optional[Union[str, torch.device]] = None, ellipsoid: str = 'sphere') Ring-based Spherical Harmonic Transform for full-sky HEALPix maps. :Parameters: * **level** (:py:class:`int Grid4Earth/HEALPix level (integer >= 0).`) -- Internally, nside = 2**level. * **lmax** (:py:class:`int` or :py:class:`None Maximum multipole. Default`: :py:class:`3·2**level - 1.`) * **dtype** (:py:class:`torch.dtype torch.float32` or :py:class:`torch.float64.`) * **device** (:py:class:`str` or :py:class:`device Defaults` to :py:class:`CUDA when available.`) * **ellipsoid** (str ``"sphere"`` (default) or ``"WGS84"``:py:class:`.`) .. py:attribute:: level .. py:attribute:: nside .. py:attribute:: lmax .. py:attribute:: dtype :value: Ellipsis .. py:attribute:: cdtype :value: Ellipsis .. py:attribute:: ellipsoid :value: '' .. py:attribute:: device .. py:attribute:: n_pix .. py:attribute:: n_rings .. py:attribute:: n_alm .. py:method:: map2alm(im: ArrayLike, nest: bool = False) -> torch.Tensor Spin-0 analysis: HEALPix map → spherical harmonic coefficients. a_lm = ∫ f(θ,φ) · Y_lm*(θ,φ) dΩ Fully differentiable with respect to ``im``. :Parameters: * **im** (:py:class:`(...`, :py:class:`N) real array-like / Tensor. RING order by default.`) * **nest** (:py:class:`bool Set True for NESTED input.`) :returns: **alm** (:py:class:`(...`, :py:class:`K) complex Tensor.`) -- K = (lmax+1)·(lmax+2)//2. Flat layout: ``[m=0: l=0..lmax | m=1: l=1..lmax | … | m=lmax]``. .. py:method:: alm2map(alm: ArrayLike, nest: bool = False) -> torch.Tensor Spin-0 synthesis: a_lm → HEALPix map. f(θ,φ) = Σ_{lm} a_lm · Y_lm(θ,φ) Fully differentiable with respect to ``alm``. :Parameters: * **alm** (:py:class:`(...`, :py:class:`K) complex array-like / Tensor.`) * **nest** (:py:class:`bool Set True for NESTED output.`) :returns: **im** (:py:class:`(...`, :py:class:`N) real Tensor.`) .. py:method:: map2alm_spin(Q: ArrayLike, U: ArrayLike, spin: int = 2, nest: bool = False) -> Tuple[torch.Tensor, torch.Tensor] 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** (:py:class:`(...`, :py:class:`N) real Tensors. RING order by default.`) * **spin** (:py:class:`int Spin weight (1` or :py:class:`2; use 0 for scalar pair).`) * **nest** (:py:class:`bool Set True for NESTED input.`) :returns: **almE, almB** (:py:class:`(...`, :py:class:`K) complex Tensors.`) .. py:method:: alm2map_spin(almE: ArrayLike, almB: ArrayLike, spin: int = 2, nest: bool = False) -> Tuple[torch.Tensor, torch.Tensor] 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** (:py:class:`(...`, :py:class:`K) complex Tensors.`) * **spin** (:py:class:`int Spin weight.`) * **nest** (:py:class:`bool Set True for NESTED output.`) :returns: **Q, U** (:py:class:`(...`, :py:class:`N) real Tensors.`) .. py:method:: anafast(im: ArrayLike, map2: Optional[ArrayLike] = None, spin: int = 0, nest: bool = False) -> torch.Tensor 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** (:py:class:`(...`, :py:class:`N)` or :py:class:`(...`, ``2``, :py:class:`N) array-like.`) -- Scalar map, or [Q, U] stacked on axis -2 for spin > 0. * **map2** (:py:class:`same shape as im` or :py:class:`None. If given`, :py:class:`compute cross-spectrum.`) * **spin** (:py:class:`int 0 for scalar`, ``1`` or :py:class:`2 for polarisation.`) * **nest** (:py:class:`bool`) :returns: **cl** (:py:class:`(...`, :py:class:`lmax+1) float Tensor for spin=0.`) -- (..., 3, lmax+1) float Tensor for spin>0 — [EE, BB, EB]. .. py:method:: uv_to_curl_div(u: ArrayLike, v: ArrayLike, nest: bool = False) -> Tuple[torch.Tensor, torch.Tensor] 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** (:py:class:`(...`, :py:class:`N) real arrays — east` and :py:class:`north components.`) :returns: **div, curl** (:py:class:`(...`, :py:class:`N) real Tensors.`)