Maps and Spherical Harmonics#

The litebird_sim.maps_and_harmonics module provides robust data containers for HEALPix maps and Spherical Harmonic coefficients ($a_{ell m}$).

While libraries like healpy provide low-level array manipulation, this module introduces an object-oriented layer designed to prevent common errors in cosmological simulation pipelines, such as mixing coordinate systems or misinterpreting array shapes.

Design Philosophy#

In standard pipelines, a map is often just a NumPy array. This leads to ambiguity: is an array of shape (3, N) a T/Q/U map, or three T-only maps? Is the map in Galactic or Ecliptic coordinates?

The classes HealpixMap and SphericalHarmonics solve this by:

Strict Shape Enforcement: Data is always stored as (nstokes, size) for single-frequency, or (nfreqs, nstokes, size) for multi-frequency. A temperature-only map is (1, npix), never (npix,).
Metadata Propagation: Every object carries physical units and coordinates. Algebraic operations (like map_a + map_b) automatically check that these match, raising a ValueError if you try to add a Galactic map to an Ecliptic one.
Multi-Frequency Support: Both classes support storing data for multiple frequencies simultaneously, with automatic validation that frequency arrays match during operations.
Backend Agnosticism: While the API feels like Healpy, the heavy lifting (SHTs) is performed by ducc0, offering significant performance improvements and correct handling of spin-weighted transforms.

Data Structures#

Healpix Maps#

The HealpixMap class wraps a dense HEALPix map.

Note

This class does not perform reordering (RING vs NESTED) internally. It stores the ordering flag nest as metadata. Users must ensure they are using the correct geometry for their specific analysis tools.

Key features:

Geometry checks: The number of pixels is validated against the nside upon initialization.
Algebra: Supports +, -, * (scalar), and Stokes vector multiplication.
Multi-frequency support: Can store maps for multiple frequencies with shape (nfreqs, nstokes, npix). The frequencies_ghz attribute holds the frequency array, and nfreqs indicates the number of frequencies.
Static Helpers: Access HEALPix geometry math without importing healpy (e.g., nside_to_resolution_rad()).

Spherical Harmonics#

The SphericalHarmonics class wraps $a_{ell m}$ coefficients. It solves the “triangular array ambiguity” by explicitly storing lmax and mmax alongside the coefficients.

Storage: Standard HEALPix/ducc triangular layout.
Convolution: The convolve() method allows easy application of beams or transfer functions in harmonic space. Supports flexible multi-frequency convolution with 1D (broadcast to all frequencies), 2D (TEB filters), or 3D (per-frequency-stokes) filter arrays.
Resizing: The resize_alm() method allows you to truncate or zero-pad coefficients to a new $ell_{max}$.
Multi-frequency support: Can store coefficients for multiple frequencies with shape (nfreqs, nstokes, nalms). The frequencies_ghz attribute holds the frequency array, and nfreqs indicates the number of frequencies.

Spherical Harmonics Ordering#

A frequent source of confusion in CMB analysis is the ordering of $a_{ell m}$ coefficients in the 1D array.

The SphericalHarmonics class strictly enforces the Standard Healpix Ordering (also known as m-major).

Memory Layout

The coefficients are stored sequentially by iterating over $m$ (outer loop) and then $ell$ (inner loop). This allows efficient recursion for Legendre polynomial calculations.

Block $m=0$: contains $ell = 0, 1, dots, ell_{max}$
Block $m=1$: contains $ell = 1, 2, dots, ell_{max}$
…
Block $m=m_{max}$: contains $ell = m_{max}, dots, ell_{max}$

Indexing Formula

To find the index of a specific mode $(ell, m)$, the class provides a static method SphericalHarmonics.get_index(). It implements the standard formula:

\[\text{index} = m \cdot \frac{2 \ell_{max} + 1 - m}{2} + \ell\]

Example: Converting coordinates to indices

from litebird_sim import SphericalHarmonics
import numpy as np

lmax = 10

# Get index for a single mode (l=2, m=2)
# Note: Arguments are (lmax, l, m) to match healpy signature
idx = SphericalHarmonics.get_index(lmax, 2, 2)

# Vectorized usage
ls = np.array([2, 3, 4])
ms = np.array([2, 2, 2])
indices = SphericalHarmonics.get_index(lmax, ls, ms)

Transforms (SHT)#

We provide high-level wrappers around ducc0 for spherical harmonic transforms. These functions handle the complexity of spin-0 (Temperature) vs spin-2 (Polarization) transforms automatically.

estimate_alm(): Map $rightarrow$ $a_{ell m}$ (Analysis)
pixelize_alm(): $a_{ell m}$ $rightarrow$ Map (Synthesis)
interpolate_alm(): $a_{ell m}$ $rightarrow$ Values at arbitrary $(theta, phi)$
compute_cl(): Compute auto- or cross-power spectra from $a_{ell m}$
compute_dl(): Compute $D_{ell} = ell(ell+1)C_{ell}/(2pi)$ spectra

Tip

All transform functions accept a nthreads argument. Setting nthreads=0 (default) uses all available hardware threads, which is optimal for standalone scripts but should be adjusted when running inside an MPI environment.

Note

All transform functions support multi-frequency data. When operating on multi-frequency objects, transforms are applied independently to each frequency, and the output maintains the multi-frequency structure.

Cookbook#

Basic I/O and Algebra#

from litebird_sim.maps_and_harmonics import HealpixMap, SphericalHarmonics
from litebird_sim.constants import Units
from litebird_sim.coordinates import CoordinateSystem

# Loading a map from disk (wrapper around hp.read_alm usually not needed for maps,
# but classes support direct instantiation)
# Here we create a dummy map for demonstration
import numpy as np

nside = 128
npix = HealpixMap.nside_to_npix(nside)

# Create a Temperature-only map in Galactic coordinates
m_gal = HealpixMap(
    values=np.zeros((1, npix)),
    nside=nside,
    units=Units.K_CMB,
    coordinates=CoordinateSystem.Galactic
)

# Create a mask (unitless)
mask = HealpixMap(
    values=np.ones((1, npix)),
    nside=nside,
    units=Units.None,
    coordinates=CoordinateSystem.Galactic
)

# Multiplication works (Units.K_CMB * Units.None = Units.K_CMB)
masked_map = m_gal * mask

Smoothing a Polarization Map#

This example demonstrates how to smooth a T, Q, U map using a Gaussian beam.

from litebird_sim.maps_and_harmonics import estimate_alm, pixelize_alm
from litebird_sim.beam_synthesis import gauss_bl

# 1. Start with a 3-component map (T, Q, U)
# shape must be (3, npix)
m_pol = HealpixMap(..., nside=64, units=Units.uK_CMB)

# 2. Analysis: Convert to alms
# This automatically handles spin-0 for T and spin-2 for Q,U
alms = estimate_alm(m_pol, lmax=128)

# 3. Create a beam window function B_ell
# gauss_bl returns an array of size lmax+1
fwhm_rad = np.radians(1.0)
b_ell = gauss_bl(fwhm_rad, lmax=128)

# 4. Convolve (apply beam)
# The 'convolve' method applies the filter to all Stokes components
# if a single array is passed.
alms_smoothed = alms.convolve(b_ell)

# 5. Synthesis: Convert back to map
m_smoothed = pixelize_alm(alms_smoothed, nside=64)

Interpolating Alms#

If you need to evaluate a set of alms at specific locations (e.g., catalog positions):

from litebird_sim.maps_and_harmonics import interpolate_alm

# Define locations: (theta, phi) in radians
# Shape must be (N, 2)
locations = np.array([
    [1.57, 0.0],  # Equator, phi=0
    [0.0, 0.0]    # North Pole
])

# Interpolate
# If alms.nstokes == 3, this returns tuple (T, Q, U)
# where each is an array of length N (single-freq) or (nfreqs, N) (multi-freq)
t_vals, q_vals, u_vals = interpolate_alm(alms, locations)

Multi-Frequency Operations#

Both HealpixMap and SphericalHarmonics support multi-frequency data:

from litebird_sim.maps_and_harmonics import SphericalHarmonics, compute_cl
import numpy as np

# Create multi-frequency alms
lmax = 128
nfreqs = 3
frequencies = np.array([100.0, 143.0, 217.0])  # GHz

# Shape: (nfreqs, nstokes, nalms)
nalms = SphericalHarmonics.alm_array_size(lmax, lmax, 3)[1]
values = np.random.randn(nfreqs, 3, nalms) + 1j * np.random.randn(nfreqs, 3, nalms)

alms = SphericalHarmonics(
    values=values,
    lmax=lmax,
    frequencies_ghz=frequencies
)

# Multi-frequency smoothing: apply different FWHM to each frequency
fwhm_arcmin = np.array([30.0, 23.0, 14.0])
fwhm_rad = np.radians(fwhm_arcmin / 60.0)
alms_smooth = alms.apply_gaussian_smoothing(fwhm_rad=fwhm_rad)

# Compute power spectra: returns 2D arrays (nfreqs, lmax+1)
cls = compute_cl(alms_smooth, alms_smooth)
# cls["TT"] has shape (3, 129) for 3 frequencies
# cls["frequency_ghz"] contains the frequency array

# Arithmetic operations check frequency compatibility
alms2 = alms.copy()
alms_sum = alms + alms2  # Works: same frequencies

# This would raise ValueError:
# alms_bad = SphericalHarmonics(..., frequencies_ghz=np.array([95, 150, 220]))
# alms + alms_bad  # ValueError: frequencies not compatible!

API Reference#

class litebird_sim.maps_and_harmonics.SphericalHarmonics(values: ndarray, lmax: int, mmax: int | None = None, nfreqs: int | None = None, frequencies_ghz: ndarray | None = None, units: Units | None = None, coordinates: CoordinateSystem | None = None)#

Bases: object

A class that wraps spherical harmonics coefficients

The convention used in libraries like HealPy is to keep the a_ℓm coefficients of a spherical harmonic expansion in a plain NumPy array. However, this is ambiguous because it is not possible to uniquely determine the value of ℓ_max and m_max from the size of the array (unless you assume that ℓ_max == m_max).

This small data class keeps the array and the values l_max and m_max together. The shape of values is (nstokes, ncoeff) if nfreqs == 1, even if nstokes == 1. For multi-frequency data (nfreqs > 1), the shape is (nfreqs, nstokes, ncoeff).

It also provides algebraic operations and I/O utilities compatible with Healpy.

Ordering Convention (Healpix Scheme): The coefficients are stored in m-major order. 1. The outer loop iterates over m from 0 to mmax. 2. The inner loop iterates over l from m to lmax.

Sequence of coefficients in the array: Index 0: (l=0, m=0) Index 1: (l=1, m=0) … Index lmax: (l=lmax, m=0) Index lmax+1: (l=1, m=1) <– Note: l starts at m Index lmax+2: (l=2, m=1) …

This specific ordering is required for operations involving healpy.alm2map or ducc0.sht.

values#

The spherical harmonics coefficients, stored in a NumPy array of shape (nstokes, ncoeff) if nfreqs == 1, or (nfreqs, nstokes, ncoeff) if nfreqs > 1.

Type:: np.ndarray

lmax#

The maximum degree ℓ_max of the expansion.

Type:: int

mmax#

The maximum order m_max of the expansion. If None, it is set equal to lmax.

Type:: int, optional

frequencies_ghz#

Array of frequencies in GHz. Must be None if nfreqs == 1, otherwise must be an array of length nfreqs.

Type:: np.ndarray or None

units#

Physical units of the coefficients. If set to :data:None, the object is treated as unitless / unspecified.

Type:: Units or None

coordinates#

Sky coordinate system of these coefficients (e.g. Galactic or Ecliptic). If None, coordinates are unspecified.

Type:: CoordinateSystem or None

nstokes#

The number of Stokes parameters (1 for intensity-only, 3 for TEB).

Type:: int

nfreqs#

The number of frequency channels (1 for single frequency, > 1 for multi-frequency).

Type:: int

Arithmetic#

----------

The following operations are supported

- `+`, `-` between two SphericalHarmonics (same `lmax`, `mmax`, `nstokes`: and compatible units / coordinates)

- `*` with scalar or Stokes-vector (array of shape `(nstokes,)`): with optional units override: sh * a or sh.__mul__(a, units=...)

- `.convolve(f_ell)` applies a filter f_ell(ℓ) or f_ell^i(ℓ) per Stokes: (units are preserved)

allclose(other: SphericalHarmonics, rtol: float = 1e-05, atol: float = 1e-08) → bool#

Check if two SphericalHarmonics objects are equal within a tolerance.

Units must be compatible (same, or one/both None/Units.None). Coordinates must also be compatible (same, or one/both None).

Returns:: True if the objects are structurally consistent and their values are close within the specified tolerance. If the objects are inconsistent (e.g., different lmax, units, etc.), a warning is logged with the specific reason, and False is returned.
Return type:: bool

static alm_array_size(lmax: int, mmax: int | None = None, nstokes: int = 3) → tuple[int, int]#

Computes the expected shape of the a_ℓm array.

Returns:: The expected shape (nstokes, ncoeff).
Return type:: tuple[int, int]

static alm_l_array(lmax: int, mmax: int | None = None) → ndarray#: Return the ℓ values corresponding to each a_{ℓm} coefficient in Healpy’s flattened alm format.

apply_gaussian_smoothing(fwhm_rad: float | ndarray | list, inplace: bool = True) → SphericalHarmonics#

Apply Gaussian smoothing to the spherical harmonics coefficients.

Parameters:

fwhm_rad (float, np.ndarray, or list) –
Full Width at Half Maximum (FWHM) of the Gaussian beam in radians. Can be:
- A single float: Same FWHM applied to all frequencies.
- An array or list of floats: Per-frequency FWHM. Must have length equal to nfreqs.
inplace (bool, optional) – If True, modifies the object in place. Default is True.

Returns:

The smoothed object.

Return type:

SphericalHarmonics

apply_pixel_window(nside: int, inplace: bool = True) → SphericalHarmonics#

Apply the HEALPix pixel window function.

Parameters:

nside (int) – The HEALPix Nside resolution parameter.
inplace (bool, optional) – If True, modifies the object in place. Default is True.

Returns:

The object with pixel window applied.

Return type:

SphericalHarmonics

convolve(f_ell: ndarray | list[ndarray], inplace: bool = True) → SphericalHarmonics#

Apply a beam or filter to the SH coefficients.

Parameters:

f_ell (np.ndarray or list[np.ndarray]) –
The ℓ-dependent filter(s). Can be:
- A 1D array: Single filter applied to all alms regardless of nfreqs.
- A list of 3 filters or 2D array with shape (3, >=lmax+1): TEB filters. If nfreqs==1, applied as-is. If nfreqs>1, looped over frequencies and applied to each frequency’s TEB alms.
- A 3D array with shape (nfreqs, nstokes, >=lmax+1): Per-frequency windows. First two dimensions must match values shape.
inplace (bool, optional) – If True, modifies the coefficients of the current object. If False, returns a new SphericalHarmonics object. Default is True.

Returns:

The object itself (if inplace=True) or a new object (if inplace=False).

Return type:

SphericalHarmonics

copy() → SphericalHarmonics#: Returns a deep copy of this SphericalHarmonics object.

Calculate the 1D linear index in the standard Healpix (m-major) layout for the given (l, m) pairs.

Formula: index = m * (2 * lmax + 1 - m) / 2 + l

Parameters:

lmax (int) – Maximum degree of the expansion.
l_val (int or np.ndarray) – Degree l.
m_val (int or np.ndarray) – Order m.

Returns:

The 1D index. Returns an integer if inputs are scalars, otherwise a numpy array.

Return type:

int or np.ndarray

get_lm_arrays() → tuple[ndarray, ndarray]#

Generate arrays of l and m for every coefficient stored in this object.

This reconstructs the explicit (l, m) coordinates for the flattened data array, following the m-major ordering used internally.

Returns:

l_arr (np.ndarray) – Array containing the l (degree) for each element.
m_arr (np.ndarray) – Array containing the m (order) for each element.

is_consistent(other: SphericalHarmonics) → bool#: Check if two SphericalHarmonics objects are compatible for algebraic operations.

static lmax_from_num_of_alm(nalm: int, mmax: int | None = None) → int#: Returns the lmax corresponding to a given array size.

property mmax_int: int#: Returns the maximum order m_max of the expansion.

static num_of_alm_from_lmax(lmax: int, mmax: int | None = None) → int#: Computes the number of a_ℓm coefficients for given ℓ_max and m_max. If mmax is not provided, it is set equal to lmax.

property num_of_alm_per_stokes: int#: Returns the number of a_ℓm coefficients per Stokes component.

print_ordering_example()#: Prints the first few (l, m) pairs to demonstrate ordering.

static read_fits(filename: str) → SphericalHarmonics#: Load a SphericalHarmonics object from a FITS file. Compatible with both Explicit Index format and Standard Healpix format.

resize_alm(lmax_out: int, mmax_out: int | None = None, inplace: bool = False)#

Resizes the spherical harmonics coefficients, either truncating or padding them with zeros.

Units, coordinates, and frequencies are preserved.

write_fits(filename: str, overwrite: bool = True, idx_frequency: int | None = None)#

Save the SphericalHarmonics object to a FITS file using the Explicit Index scheme.

Format: - Separate HDUs for TEMPERATURE, E_MODE, B_MODE. - Each HDU has 3 columns: INDEX, REAL, IMAG. - INDEX = l^2 + l + m + 1 (FITS standard for sparse alm).

Parameters:

filename (str) – The FITS file path.
overwrite (bool, optional) – If True, overwrite existing file. Default is True.
idx_frequency (int or None, optional) – The frequency index to write (for multi-frequency objects). - If nfreqs is None (single-frequency): ignored; raises warning if not None. - If nfreqs is not None (multi-frequency): must be specified and < nfreqs. Default is None.

Raises:

ValueError – If nfreqs > 1 and idx_frequency is None. If idx_frequency is not in valid range [0, nfreqs).
UserWarning – If nfreqs is None and idx_frequency is not None.

Create a SphericalHarmonics object filled with zeros.

Parameters:

lmax (int) – Maximum degree.
mmax (int, optional) – Maximum order. Defaults to lmax.
nstokes (int, default=3) – Number of Stokes parameters (1 or 3).
nfreqs (int or None, default=None) – Number of frequency channels. If None (default), creates a single-frequency object. If int > 0, creates a multi-frequency object. If frequencies_ghz is provided, nfreqs is automatically computed and overrides this parameter.
frequencies_ghz (float, Sequence[float], np.ndarray, or None, optional) – Array of frequencies in GHz. Can be a single float or an array of floats. If provided, nfreqs is automatically set based on the size of frequencies_ghz. If not provided, nfreqs must be explicitly set or defaults to None (single-frequency mode).
dtype (type, default=np.complex128) – Data type for the array. Must be np.complex64 or np.complex128.
units (Units, optional) – Physical units.
coordinates (CoordinateSystem, optional) – Sky coordinate system.

Returns:

Instance initialized with zeros.

Return type:

SphericalHarmonics

Raises:

ValueError – If dtype is not np.complex64 or np.complex128.

class litebird_sim.maps_and_harmonics.HealpixMap(values: ndarray, nside: int | None = None, nfreqs: int | None = None, frequencies_ghz: ndarray | None = None, units: Units | None = None, coordinates: CoordinateSystem | None = None, nest: bool = False)#

Bases: object

A small container class for HEALPix maps, with shape checks and algebra.

This class stores one- or three-component HEALPix maps and their NSIDE consistently. Maps are always stored as a 2D NumPy array of shape (nstokes, npix), even if nstokes == 1.

If nside is not provided, it is inferred from the size of values.

values#

Map values, stored as a NumPy array of shape (nstokes, npix) if nfreqs == 1, or (nfreqs, nstokes, npix) if nfreqs > 1. If a 1D array is passed, it is promoted to shape (1, npix). If a tuple of arrays is passed (e.g. (I, Q, U)), it is stacked along the first axis.

Type:: np.ndarray

nside#

HEALPix NSIDE resolution parameter. If None, it is automatically inferred from the last dimension of values. If provided, it is checked against the size of values.

Type:: int | None

frequencies_ghz#

Array of frequencies in GHz. Must be None if nfreqs == 1, otherwise must be an array of length nfreqs.

Type:: np.ndarray or None

units#

Physical units of the map. If set to Units.None or None, the map is treated as unitless / unspecified.

Type:: Units or None

coordinates#

Sky coordinate system of the map (e.g. Galactic or Ecliptic). If None, coordinates are unspecified.

Type:: CoordinateSystem or None

nest#

If True, the map is in NESTED ordering; if False, in RING (default). This is just metadata here, no indexing operations are performed.

Type:: bool, optional

nstokes#

Number of Stokes parameters (1 for intensity-only, 3 for IQU).

Type:: int

nfreqs#

The number of frequency channels (1 for single frequency, > 1 for multi-frequency).

Type:: int

allclose(other: HealpixMap, rtol: float = 1e-05, atol: float = 1e-08) → bool#

Check if two HealpixMap objects are equal within a tolerance.

Checks for consistency in geometry (nside, nest, nstokes), metadata (units, coordinates), and numerical values.

Returns:: True if the objects are structurally consistent and their values are close within the specified tolerance. If the objects are inconsistent, a warning is logged with the specific reason, and False is returned.
Return type:: bool

copy() → HealpixMap#: Return a deep copy of this HealpixMap object.

is_consistent(other: HealpixMap) → bool#

Check if two HealpixMap objects are compatible for algebraic operations.

Two maps are considered consistent if they share the same nside, nest, nstokes, and nfreqs.

static is_npix_ok(num_of_pixels) → bool#

Return True or False whenever num_of_pixels is a valid number.

The function checks if the number of pixels provided as an argument conforms to the Healpix standard, which means that the number must be in the form 12 * NSIDE^2, where NSIDE is a power of 2.

>>> from litebird_sim import HealpixMap
>>> HealpixMap.is_npix_ok(48)
True
>>> HealpixMap.is_npix_ok(108) # 12 * 3^2, but 3 not power of 2
False

property npix: int#: Return the number of pixels in the map.

static npix_to_nside(num_of_pixels: int) → int#

Return NSIDE for a Healpix map containing num_of_pixels pixels.

If the number of pixels does not conform to the Healpix standard, an AssertionError exception is raised.

>>> from litebird_sim import HealpixMap
>>> HealpixMap.npix_to_nside(48)
2

static nside_to_npix(nside: int) → int#

Return the number of pixels in a Healpix map with the specified NSIDE.

If the value of nside is not valid (power of two), an AssertionError exception is raised.

>>> from litebird_sim import HealpixMap
>>> HealpixMap.nside_to_npix(1)
12

static nside_to_pixel_solid_angle_sterad(nside: int) → float#

Return the value of the solid angle of a pixel.

The result is exact, as all pixels in a Healpix map have the same area.

The result is in steradians.

static nside_to_resolution_rad(nside: int) → float#

Return an approximated resolution of a Healpix map, given its NSIDE.

The value is the square root of the pixel area (which is measured in steradians); see HealpixMap.nside_to_pixel_solid_angle_sterad().

The result is an angle in radians.

classmethod zeros(nside: int, nstokes: int = 3, nfreqs: int | None = None, frequencies_ghz: float | Sequence[float] | ndarray | None = None, dtype: type = <class 'numpy.float64'>, units: Units | None = None, coordinates: CoordinateSystem | None = None, nest: bool = False) → HealpixMap#

Create a HealpixMap object filled with zeros.

Parameters:

nside (int) – HEALPix resolution.
nstokes (int, default=3) – Number of Stokes parameters.
nfreqs (int | None, default=None) – Number of frequency channels. If None, creates single-frequency map (2D); if int > 0, creates multi-frequency map (3D). If frequencies_ghz is provided, nfreqs is automatically computed and overrides this parameter.
frequencies_ghz (float, Sequence[float], np.ndarray, or None, optional) – Array of frequencies in GHz. Can be a single float or an array of floats. If provided, nfreqs is automatically set based on the size of frequencies_ghz. If not provided, nfreqs must be explicitly set or defaults to None (single-frequency mode).
dtype (type, default=np.float64) – Data type of the map.
units (Units, optional) – Physical units.
coordinates (CoordinateSystem, optional) – Sky coordinate system.
nest (bool, default=False) – Ordering (Ring vs Nested).

Returns:

Instance initialized with zeros.

Return type:

HealpixMap

Raises:

ValueError – If both nfreqs and frequencies_ghz are inconsistent.

litebird_sim.maps_and_harmonics.estimate_alm(map: HealpixMap, *, lmax: int | None = None, mmax: int | None = None, maxiter: int | None = None, epsilon: float | None = None, nthreads: int = 0) → SphericalHarmonics#

Estimate spherical harmonic coefficients ($a_{ell m}$) from a HEALPix map.

This function performs a spherical harmonic analysis to transform the input map into harmonic space. Two backends are available depending on maxiter:

Default (``maxiter=None``) – uses ducc0.sht.adjoint_synthesis() to compute the pixel summation and scales the result by the pixel area, approximating the integral over the sphere:

\[a_{\ell m} \approx \Omega_{pix} \sum_{p} f(p) Y_{\ell m}^*(p)\]

where ducc0.sht.adjoint_synthesis() computes $sum f Y^*$ and the result is multiplied by $Omega_{pix}$ (base.pix_area()).
Iterative (``maxiter`` a positive integer) – uses ducc0.sht.pseudo_analysis(), an iterative least-squares solver that converges to more accurate coefficients at the cost of additional computation.

Both backends support scalar and polarized fields:

map.nstokes == 1 → spin-0 transform, returns a 1-component alm
map.nstokes == 3 → spin-0 for T and spin-2 for (Q, U), returns a 3-component alm (T, Q, U interpreted as spin-2)

Geometry is obtained from ducc0.healpix.Healpix_Base via its sht_info() method, which provides the arrays theta, nphi, phi0 and (optionally) ringstart needed by ducc’s SHT routines.

The units and coordinates metadata are propagated from map to the returned SphericalHarmonics.

Parameters:

map (HealpixMap) –
Input HEALPix map. Must satisfy:
- map.values.shape == (nstokes, npix)
- map.nstokes in (1, 3)
- map.nside is a power of two.
lmax (int, optional) –
Maximum multipoles to compute. Defaults are chosen in a Healpy-like way:
- if lmax is None → lmax_eff = 3 * nside - 1
- if mmax is None → mmax_eff = lmax_eff
mmax_eff must be ≤ lmax_eff.
mmax (int, optional) –
Maximum multipoles to compute. Defaults are chosen in a Healpy-like way:
- if lmax is None → lmax_eff = 3 * nside - 1
- if mmax is None → mmax_eff = lmax_eff
mmax_eff must be ≤ lmax_eff.
maxiter (int or None, optional) – If None (default), uses the non-iterative ducc0.sht.adjoint_synthesis() backend followed by pixel-area scaling. If a positive integer, uses ducc0.sht.pseudo_analysis() instead and performs at most maxiter iterations.
epsilon (float or None, optional) –
Desired accuracy for ducc0.sht.pseudo_analysis(). This parameter is used only when maxiter is a positive integer. If None, a dtype-dependent default is chosen:
- complex64 output -> 1e-6
- complex128 output -> 1e-13
nthreads (int, optional) – Number of threads passed to ducc. If zero (default), ducc chooses the number of threads.

Returns:

A spherical harmonics object with

values of shape (nstokes, nalm)
lmax = lmax_eff
mmax = mmax_eff
units = map.units
coordinates = map.coordinates

Return type:

SphericalHarmonics

Raises:

TypeError – If maxiter is not an integer or None.
ValueError – If maxiter is an integer but not positive (i.e. ≤ 0). If the HealpixMap object has unsupported nstokes, or if requested lmax/mmax are inconsistent.

litebird_sim.maps_and_harmonics.pixelize_alm(alms: SphericalHarmonics, nside: int, *, nest: bool = False, lmax: int | None = None, mmax: int | None = None, nthreads: int = 0) → HealpixMap#

Convert spherical harmonics coefficients to a HEALPix map using ducc0.sht.synthesis() on the HEALPix geometry.

This is essentially an pixelize_alm implemented on top of ducc0, with support for both scalar and polarized fields:

alms.nstokes == 1 → spin-0 synthesis, returns a 1-component map
alms.nstokes == 3 → spin-0 for T and spin-2 for (Q, U), returns a 3-component map

Geometry is obtained from ducc0.healpix.Healpix_Base via its sht_info() method, which provides the arrays theta, nphi, phi0 and (optionally) ringstart needed by ducc’s SHT routines.

The units and coordinates metadata are propagated from alms to the returned HealpixMap.

Parameters:

alms (SphericalHarmonics) –
Harmonic coefficients. Must satisfy:
- alms.values.shape == (nstokes, nalm)
- alms.nstokes in (1, 3)
- triangular Healpy/ducc layout in the last dimension (i.e. increasing m, l ≥ m).
nside (int) – Target HEALPix resolution (power of two). The output map will have npix = 12 * nside**2 pixels.
nest (bool, optional) – If False (default), the map is in RING ordering. If True, the map is in NESTED ordering. The underlying SHT is always performed in RING geometry; for NESTED, the geometry is built from a NESTED Healpix_Base but pixel ordering in the returned HealpixMap is marked accordingly.
lmax (int, optional) –
Maximum multipoles to use in the synthesis. If omitted, the values from alms are used:
- lmax_eff = alms.lmax if lmax is None
- mmax_eff = alms.mmax if mmax is None
Both must be ≤ the corresponding values in alms. If a smaller lmax/mmax is requested than what is stored in alms, the coefficient array is truncated accordingly.
mmax (int, optional) –
Maximum multipoles to use in the synthesis. If omitted, the values from alms are used:
- lmax_eff = alms.lmax if lmax is None
- mmax_eff = alms.mmax if mmax is None
Both must be ≤ the corresponding values in alms. If a smaller lmax/mmax is requested than what is stored in alms, the coefficient array is truncated accordingly.
nthreads (int, optional) – Number of threads passed to ducc. If zero (default), ducc chooses the number of threads.

Returns:

A HEALPix map with shape (nstokes, npix) and metadata

nside as specified
nest as specified
units = alms.units
coordinates = alms.coordinates

Return type:

HealpixMap

Raises:

ValueError – If the SphericalHarmonics object has unsupported nstokes, or if requested lmax/mmax are inconsistent with the available coefficients.

litebird_sim.maps_and_harmonics.interpolate_alm(alms: SphericalHarmonics, locations: ndarray, *, epsilon: float | None = None, nthreads: int = 0) → ndarray | tuple[ndarray, ndarray, ndarray]#

Interpolate spherical-harmonic coefficients at arbitrary positions using ducc0.sht.synthesis_general().

This wrapper supports both scalar and polarized fields:

If alms.nstokes == 1: a single spin-0 field is synthesized and the function returns a 1D array T of length N.
If alms.nstokes == 3: the first component is interpreted as a spin-0 field (T), and the remaining two as a spin-2 field, returning T, Q, U.

Parameters:

alms (SphericalHarmonics) –
Object exposing at least the attributes:
- values: ndarray, shape (nstokes, nalm), complex64/complex128
- lmax: int
- mmax: int
- nstokes: int (1 or 3)
The last dimension of values must follow the standard triangular (healpy/ducc) ordering.
locations (ndarray, shape (N, 2)) – Target positions on the sphere, in radians. locations[:, 0] = colatitude theta (0 .. π), locations[:, 1] = longitude phi (0 .. 2π).
epsilon (float, optional) –
Desired accuracy passed to synthesis_general(). If None, a safe default is chosen depending on the dtype:
- complex64 → 1e-6
- complex128 → 1e-13
nthreads (int, optional) – Number of threads for ducc. If 0 (default), ducc uses the number of hardware threads.

Returns:

T (ndarray) – If alms.nstokes == 1, a 1D array of length N with the interpolated scalar field.
T, Q, U (ndarray) – If alms.nstokes == 3, three 1D arrays of length N with the interpolated Stokes parameters.

Raises:

ValueError – If the input shapes are inconsistent or nstokes is not 1 or 3.