"""Angular power spectrum C_ℓ via healpy pseudo-Cl estimator."""
from __future__ import annotations
import healpy as hp
import numpy as np
from ..catalogue import GalaxyCatalogue
[docs]
def cl_angular(
gal: GalaxyCatalogue,
nside: int,
ell_bins: np.ndarray,
mask: np.ndarray | None = None,
use_master_correction: bool = False,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
"""Angular power spectrum C_ℓ of a galaxy overdensity map.
The galaxy catalogue is pixelised onto a HEALPix map at resolution
``nside``, converted to an overdensity field δ = n/n̄ − 1 (weighted
by galaxy weights), and the pseudo-Cl is estimated via ``healpy.anafast``.
Shot noise is subtracted analytically.
Parameters
----------
gal : GalaxyCatalogue
Galaxy catalogue (ra, dec, weight required).
nside : int
HEALPix resolution (power of 2, e.g. 512 or 1024).
ell_bins : ndarray, shape (n_bins+1,)
Multipole bin edges (integers).
mask : ndarray, shape (n_pix,), optional
HEALPix binary or apodised mask. Defaults to a mask derived from
occupied pixels.
use_master_correction : bool
Apply an analytical MASTER correction for the mask (approximate).
Returns
-------
ell_eff : ndarray, shape (n_bins,)
Effective multipole per bin (simple mean).
cl : ndarray, shape (n_bins,)
Angular power spectrum [sr].
nl : ndarray, shape (n_bins,)
Shot noise power spectrum [sr].
cov : ndarray, shape (n_bins, n_bins)
Gaussian covariance of C_ℓ (diagonal, approximate).
Performance
-----------
~500 ms/call (N_gal=1e5, nside=512, CPU x86-64)
"""
n_pix = hp.nside2npix(nside)
omega_pix = hp.nside2pixarea(nside) # sr per pixel
# Pixelise the weighted galaxy map
pix_ids = hp.ang2pix(nside, gal.ra, gal.dec, lonlat=True)
delta_map = np.zeros(n_pix)
weight_map = np.zeros(n_pix)
for i in range(gal.n):
delta_map[pix_ids[i]] += gal.weight[i]
weight_map[pix_ids[i]] += gal.weight[i]
# Build mask from occupied pixels if not provided
if mask is None:
mask = (weight_map > 0).astype(np.float64)
occupied = mask > 0
n_gal_eff = np.sum(gal.weight[occupied[pix_ids]])
n_pix_eff = np.sum(mask)
# Mean density per pixel
n_bar = n_gal_eff / n_pix_eff if n_pix_eff > 0 else 1.0
# Overdensity field
overdensity = np.zeros(n_pix)
overdensity[occupied] = delta_map[occupied] / n_bar - 1.0
overdensity *= mask
# Survey fraction
f_sky = n_pix_eff / n_pix
# Pseudo-Cl
lmax = int(ell_bins[-1]) + 1
cl_full = hp.anafast(overdensity, lmax=lmax)
# Mode-coupling correction (approximate MASTER)
if use_master_correction and f_sky > 0:
cl_full = cl_full / f_sky
# Shot noise per mode: N_ℓ = Ω_pix / N_gal (per steradian)
nl_val = omega_pix / n_gal_eff if n_gal_eff > 0 else 0.0
# Bin the power spectrum
ell_bins = np.asarray(ell_bins, dtype=int)
n_bins = len(ell_bins) - 1
ells = np.arange(len(cl_full))
ell_eff = np.zeros(n_bins)
cl_binned = np.zeros(n_bins)
nl_binned = np.zeros(n_bins)
cov_binned = np.zeros((n_bins, n_bins))
for i in range(n_bins):
lo, hi = ell_bins[i], ell_bins[i + 1]
sel = (ells >= lo) & (ells < hi)
n_modes = np.sum(sel)
if n_modes > 0:
ell_eff[i] = np.mean(ells[sel])
cl_binned[i] = np.mean(cl_full[sel]) - nl_val
nl_binned[i] = nl_val
# Gaussian variance: Var(C_ℓ) = 2 (C_ℓ + N_ℓ)^2 / (2ℓ+1) / f_sky
cl_tot = cl_binned[i] + nl_val
var_cl = 2.0 * cl_tot**2 / (n_modes * f_sky)
cov_binned[i, i] = var_cl
return ell_eff, cl_binned, nl_binned, cov_binned