lenspyx.wigners

Wigner small-d functions and correlation functions for spin-weighted spherical harmonics.

This module provides efficient implementations of Wigner small-d functions and their associated correlation functions, which are fundamental for analyzing spin-weighted fields on the sphere (e.g., CMB polarization, gravitational lensing).

Overview

Wigner small-d functions \(d^\ell_{s_1 s_2}(\theta)\) are the angular parts of the Wigner D-matrices, describing the rotation of spin-weighted spherical harmonics. They appear in:

Correlation functions of spin-weighted fields (e.g., ξ± for CMB polarization or galaxy surveys)
Angular power spectrum estimators for CMB polarization
Lensing reconstruction and delensing operations
General spin transformations on the sphere

API Reference

Forward Transforms

lenspyx.wigners.wignerpos(cl: ndarray[float], theta: ndarray[float], s1: int, s2: int)[source]

Produces Wigner small-d transform defined by

\[\sum_\ell \frac{2\ell + 1}{4\pi} C_\ell d^\ell_{s_1 s_2}(\theta)\]

Parameters:

cl (array_like) – Spectrum of Wigner small-d transform (power spectrum \(C_\ell\))
theta (array_like) – Co-latitude in radians (in [0, π])
s1 (int) – First spin weight
s2 (int) – Second spin weight

Returns:

Real array of same size as theta containing the correlation function

Return type:

array_like

Notes

You can use wigner4pos() instead if you also need the result for -s2 (e.g., for computing \(\xi_{\pm}\) simultaneously).

See also

wigner4pos: Compute 4 correlation functions in one go
wignercoeff: Adjoint transform (correlation function to spectrum)

lenspyx.wigners.wigner4pos(gl: ndarray[float], cl: ndarray[float], theta: ndarray[float], s1: int, s2: int)[source]

Compute 4 Wigner correlation functions in one go.

Parameters:

gl (array_like) – First spectrum (power spectrum \(g_\ell\))
cl (array_like or None) – Second spectrum (power spectrum \(c_\ell\)). Can be set to None if irrelevant, and is ignored if s2 is zero.
theta (array_like) – Co-latitude in radians (in [0, π])
s1 (int) – First spin weight
s2 (int) – Second spin weight

Returns:

In the most general case, an array of shape (ncomp, ntheta) with:

Component 0: \(\sum_\ell g_\ell \frac{2\ell + 1}{4\pi} d^\ell_{s_1, |s_2|}(\theta)\)
Component 1: \(\sum_\ell g_\ell \frac{2\ell + 1}{4\pi} d^\ell_{s_1,-|s_2|}(\theta)\)
Component 2: \(\sum_\ell c_\ell \frac{2\ell + 1}{4\pi} d^\ell_{s_1, |s_2|}(\theta)\) (if cl is not None)
Component 3: \(\sum_\ell c_\ell \frac{2\ell + 1}{4\pi} d^\ell_{s_1,-|s_2|}(\theta)\) (if cl is not None)

The number of components ncomp in the output is:

4 if (s2 ≠ 0 and cl is not None)
2 if (s2 ≠ 0 and cl is None)
1 if s2 = 0

Return type:

array_like

Notes

This function is more efficient than calling wignerpos() multiple times when you need correlation functions for both +s2 and -s2 (e.g., ξ+ and ξ-).

See also

wignerpos: Compute single correlation function

Adjoint Transforms

lenspyx.wigners.wignercoeff(xi: ndarray[float], theta: ndarray[float], s1: int, s2: int, lmax: int)[source]

Computes spectrum of Wigner small-d correlation function (adjoint to wignerpos).

This is the adjoint transform that converts a correlation function to its power spectrum:

\[C_\ell = 2\pi \sum_\theta \xi(\theta) d^\ell_{s_1 s_2}(\theta)\]

Parameters:

xi (array_like) – Wigner correlation function (real array, one value per co-latitude)
theta (array_like) – Co-latitude in radians (in [0, π])
s1 (int) – First spin weight
s2 (int) – Second spin weight
lmax (int) – Maximum multipole (inclusive). Spectrum is calculated from 0 to lmax.

Returns:

Power spectrum \(C_\ell\) from ℓ=0 to ℓ=lmax

Return type:

array_like

Notes

This function implements the adjoint operation to wignerpos(). It is used to estimate power spectra from measured correlation functions.

See also

wignerpos: Forward transform (spectrum to correlation function)

Convolutions

lenspyx.wigners.wignerc(cl1: ndarray[float], cl2: ndarray[float], s1: int, t1: int, s2: int, t2: int, lmax_out: int = -1)[source]

Convolution of two Wigner small-d correlation functions.

Computes the power spectrum of the product of two correlation functions:

\[\xi_{s_1 t_1}(\mu) \times \xi_{s_2 t_2}(\mu)\]

where

\[\xi_{s_1 t_1}(\mu) = \sum_{\ell} C_{1,\ell} \frac{2\ell + 1}{4\pi} d^\ell_{s_1 t_1}(\mu)\]

\[\xi_{s_2 t_2}(\mu) = \sum_{\ell} C_{2,\ell} \frac{2\ell + 1}{4\pi} d^\ell_{s_2 t_2}(\mu)\]

Gauss-Legendre quadrature is used to solve this exactly (up to numerical precision).

Parameters:

cl1 (array_like) – Spectrum of first Wigner small-d function (\(C_{1,\ell}\))
cl2 (array_like) – Spectrum of second Wigner small-d function (\(C_{2,\ell}\))
s1 (int) – First spin of first function
t1 (int) – Second spin of first function
s2 (int) – First spin of second function
t2 (int) – Second spin of second function
lmax_out (int, optional) – Maximum multipole of output spectrum. Defaults to len(cl1) + len(cl2) - 2.

Returns:

Power spectrum of the product correlation function, from ℓ=0 to ℓ=lmax_out

Return type:

array_like

Notes

This function is useful for computing non-Gaussian covariances and higher-order statistics. The output has spins (s1+s2, t1+t2).

Wigner d-functions

lenspyx.wigners.wignerd(l: int, s1: int, s2: int, theta: ndarray)[source]

Returns Wigner small-d functions for a specific multipole.

Computes the normalized Wigner d-functions:

\[\frac{2\ell + 1}{4\pi} d^\ell_{s_1, |s_2|}(\theta)\]

and

\[\frac{2\ell + 1}{4\pi} d^\ell_{s_1,-|s_2|}(\theta)\]

for all θ. If s2 is zero, only the first component is returned.

Parameters:

l (int) – Multipole moment
s1 (int) – First spin weight
s2 (int) – Second spin weight
theta (array_like) – Co-latitude angles in radians

Returns:

Array of shape (2, ntheta) if s2 ≠ 0, or (1, ntheta) if s2 = 0

Return type:

array_like

See also

wignerdl: Returns d^l_{s1,s2}(θ) without normalization for all ℓ

lenspyx.wigners.wignerdl(s1: int, s2: int, theta: float, lmax: int)[source]

Returns the Wigner d-function for all multipoles.

Computes the Wigner d-function:

\[d^\ell_{s_1 s_2}(\theta)\]

for all ℓ from 0 to lmax.

Parameters:

s1 (int) – First spin weight
s2 (int) – Second spin weight
theta (float) – Co-latitude angle in radians (scalar)
lmax (int) – Maximum multipole moment

Returns:

Array of size lmax+1 containing \(d^\ell_{s_1 s_2}(\theta)\) for ℓ=0 to lmax

Return type:

array_like

Notes

This returns the unnormalized Wigner d-function, without the (2ℓ+1)/(4π) factor. For the normalized version, use wignerd().

See also

wignerd: Normalized Wigner d-function for a specific ℓ

Quadrature Utilities

lenspyx.wigners.get_thgwg(npts: int)[source]

Gauss-Legendre integration points and weights from DUCC0.

Provides quadrature points and weights for integration over the interval [0, π], optimized for integrating functions on the sphere.

Parameters:

npts (int) – Number of quadrature points

Returns:

tht (array_like) – Co-latitude points in radians (array of size npts)
wg (array_like) – Quadrature weights (array of size npts)

Notes

This uses DUCC0’s highly optimized Gauss-Legendre quadrature implementation. For a band-limited function up to ℓ_max, use npts ≥ (ℓ_max + 2) / 2 for exact integration.

Examples

>>> from lenspyx.wigners import get_thgwg
>>> theta, weights = get_thgwg(100)
>>> # Integrate a function f(θ) over the sphere:
>>> integral = np.sum(f(theta) * weights)

lenspyx.wigners.get_xgwg(a: float, b: float, npts: int)[source]

Gauss-Legendre points and weights for integration over an arbitrary interval.

Provides quadrature points and weights for integration over the interval [a, b].

Parameters:

a (float) – Lower bound of integration interval
b (float) – Upper bound of integration interval
npts (int) – Number of quadrature points

Returns:

xg (array_like) – Quadrature points within (a, b) (array of size npts)
wg (array_like) – Quadrature weights (array of size npts)

Notes

This function transforms the standard Gauss-Legendre quadrature on [-1, 1] to an arbitrary interval [a, b].

Examples

>>> from lenspyx.wigners import get_xgwg
>>> x, weights = get_xgwg(0, 10, 50)
>>> # Integrate a function f(x) over [0, 10]:
>>> integral = np.sum(f(x) * weights)

Examples

Computing CMB E-mode correlation function

import numpy as np
from lenspyx.wigners import wignerpos

# E-mode power spectrum
cl_ee = np.loadtxt('cl_ee.txt')
theta = np.linspace(0, np.pi, 100)

# Compute ξ_EE(θ) = Σ_ℓ (2ℓ+1)/(4π) C_ℓ^EE d^ℓ_{2,2}(θ)
xi_ee = wignerpos(cl_ee, theta, s1=2, s2=2)

Computing ξ± correlation functions

from lenspyx.wigners import wigner4pos

# Returns [ξ+, ξ-] where:
# ξ+ = Σ_ℓ (2ℓ+1)/(4π) C_ℓ d^ℓ_{2,+2}(θ)
# ξ- = Σ_ℓ (2ℓ+1)/(4π) C_ℓ d^ℓ_{2,-2}(θ)
xi_plus, xi_minus = wigner4pos(cl_ee, None, theta, s1=2, s2=2)

Inverting correlation function to power spectrum

from lenspyx.wigners import wignercoeff, get_thgwg

lmax = 2000
npts = (lmax + 2) // 2
theta, weights = get_thgwg(npts)

# Measure ξ(θ) from data
xi_measured = measure_correlation(theta)

# Compute C_ℓ = 2π Σ_θ ξ(θ) d^ℓ_{s1,s2}(θ)
cl_recovered = wignercoeff(xi_measured, theta, s1=2, s2=2, lmax=lmax)

Convolving correlation functions

from lenspyx.wigners import wignerc

cl1 = cl_ee  # E-mode spectrum
cl2 = cl_bb  # B-mode spectrum

# Compute spectrum of ξ_EE(θ) * ξ_BB(θ)
cl_product = wignerc(cl1, cl2, s1=2, t1=2, s2=2, t2=2)

Notes

Spin values s1, s2 can be positive, negative, or zero
The normalization includes (2ℓ+1)/(4π) factors
Gauss-Legendre quadrature ensures exact integration for band-limited functions

The implementation uses DUCC0’s alm2leg and leg2alm functions, which achieve excellent performance through compiler optimization.