active

Software for Radio Cosmology

Open-source Python projects for 21cm intensity mapping — beam modelling, data simulation, Bayesian calibration, and power spectrum estimation

Software for Radio Cosmology

Open-source Python projects for 21cm intensity mapping — beam modelling, data simulation, Bayesian calibration, and power spectrum estimation

The 21cm Challenge

Cosmic 21cm signal history across redshift

Detecting the cosmological 21cm signal from neutral hydrogen is one of the great observational challenges in modern astrophysics. The signal is buried beneath astrophysical foregrounds four to five orders of magnitude brighter, corrupted by instrumental systematics — beam chromaticity, gain drifts, correlated 1/f noise — and contaminated by radio frequency interference. Every stage of the data pipeline must be handled with care — motivating the projects presented here.

Projects at a Glance

Over the past several years I have developed a collection of Python packages that address different aspects of single-dish radio cosmology — from characterizing antenna beams to estimating power spectra. Each project is self-contained, though they naturally complement one another.

block-beta
  columns 4
  TIBEC["TIBEC\nBeam Transforms"]:1
  limTOD["limTOD\nTOD Simulator"]:1
  hydra["hydra-tod\nBayesian Calibration"]:1
  comat["comat\nFast Toeplitz Ops"]:1
  QPSE["weighted QPSE\nPower Spectrum"]:1
  MERS["MERS\nForeground Modelling"]:1
  G21["Global21cm\nSignal Emulation"]:1
  space:1

  style TIBEC fill:#4a86c8,color:#fff
  style limTOD fill:#4a86c8,color:#fff
  style hydra fill:#e8744f,color:#fff
  style comat fill:#7ab648,color:#fff
  style QPSE fill:#9b59b6,color:#fff
  style MERS fill:#f1c40f,color:#333
  style G21 fill:#95a5a6,color:#fff

Full Stokes Beam Integration — TIBEC

Full Stokes antenna beam patterns (I, Q, U, V)

TIBEC (Time Integrated Beam pattern in Equatorial Coordinates) transforms antenna far-field E-field patterns into full Stokes (I, Q, U, V) beam matrices in equatorial coordinates, with support for time integration over Local Sidereal Time.

Reads FEKO/CST E-field data formats
Computes Stokes beams via Pauli matrix algebra
Two APIs: E_field (LST-based, multi-frequency) and EFieldBeamAngle (beam-angle-based, HEALPix)

TOD Simulation — limTOD

Coordinate system transformations for beam pointing

limTOD simulates realistic time-ordered data for single-dish intensity mapping experiments such as MeerKAT/MeerKLASS. It handles asymmetric beam rotation in spherical harmonic space, sky convolution via HEALPix, and correlated 1/f gain noise generation.

Full Stokes (I, Q, U, V) beam convolution
Coordinate transforms: telescope (az/el) to equatorial via ZYZYZ Euler decomposition
Built-in high-pass + Wiener filter map-making (HPW_mapmaking)
MPI parallelization over frequencies

Bayesian Calibration — hydra-tod

Reconstructed sky maps from Bayesian Gibbs sampling

hydra-tod performs joint Bayesian inference of sky temperature maps, instrument gains, and correlated noise parameters from radio telescope time-ordered data. Using Gibbs sampling with MPI parallelization, it alternates between sampling sky, gain, and noise conditional distributions.

Data model: TOD = T_sys × (1 + noise) × gain
Three gain parameterizations: linear, log-linear, factorized
Flicker noise (1/f^α) covariance via incomplete gamma functions
Polynomial emulators for ~1700× speedup of covariance evaluation
Published in RAS Techniques and Instruments (arXiv:2509.10992)

Gibbs sampling for calibration, map-making, and noise estimation

The Gibbs sampler alternates between four conditional updates per iteration. Steps (a)–(c) run independently per TOD on each MPI rank; step (d) synchronises the shared sky parameters across all ranks.

Init

Initialise
p_sky⁽⁰⁾, {p_loc,j⁽⁰⁾}, {p_g,j⁽⁰⁾}, {p_n,j⁽⁰⁾}

→

Per TOD j (parallel over MPI ranks)

(a)

Sample gains p_g,j⁽ⁱ⁾

↓

(b)

Sample local temp p_loc,j⁽ⁱ⁾

↓

(c)

Sample noise p_n,j⁽ⁱ⁾

→
sync

Global (MPI allreduce)

(d)

Sample sky temp p_sky⁽ⁱ⁾

→

i < N_samp ?
Yes ↻ No ↓

Fast Toeplitz Operations — comat

The flicker noise covariance matrices in hydra-tod are Toeplitz-structured. comat provides a Cython + OpenMP implementation of the Levinson-Durbin algorithm for O(N²) computation of:

Log-determinant: log|C| for the noise covariance matrix
Quadratic form: x^T C⁻¹ x for likelihood evaluation

These two quantities are the computational bottleneck of the noise parameter sampler. By exploiting Toeplitz structure, comat avoids the O(N³) cost of general matrix factorization, making the Gibbs sampler tractable for long time streams.

Power Spectrum Estimation — weighted QPSE

Power spectrum mode selection and k-space coverage

The weighted Quadratic Power Spectrum Estimator extracts cosmological 21cm power spectra from interferometric visibility data, with principled handling of non-uniform frequency weighting caused by RFI flagging.

Half-cosine apodization around flagged channels to minimize spectral leakage
Comprehensive window function analysis
Automatic noise bias correction
Both theoretical (Gaussian) and Monte Carlo covariance estimation

Foreground Modelling — MERS

Synchrotron emission mechanism — the dominant radio foreground

MERS (Moment Expanded Radio Sky) characterizes radio foregrounds using a power-law ensemble moment expansion formalism. It enables systematic comparison of sky models (CNN-PL vs GSM), quantifies the impact of beam chromaticity and point sources on foreground subtraction, and provides foreground mitigation strategies for both global 21cm and interferometric experiments.

Companion emulators in Global21cm (21cmVAE, Poly21cm) provide fast neural-network-based predictions of the 21cm signal for use in forecasting and analysis pipelines.

Key Takeaways

Broad coverage: Python packages addressing key challenges in single-dish 21cm intensity mapping — beam transforms, TOD simulation, Bayesian calibration, map-making, power spectrum estimation, and foreground modelling.
Full Stokes support: Beam integration, TOD simulation, and map-making all handle the full polarization state (I, Q, U, V).
Production-ready: MPI-parallelized solvers, Cython/OpenMP routines (comat), and polynomial emulators for computationally intensive operations.
Open-source: All packages available as installable Python libraries with Jupyter notebook examples and comprehensive documentation.
Actively developed: Ongoing work on MeerKLASS and RHINO telescope applications.