Simple Usage

from serval.containers import SkyModel, TIRSPowerBeam
from serval.core import SingleVisibilityGenerator
from serval import sht
from serval.utils import (
    plane_wave_mag_from_bandlimit,
    harmonic_point_source,
    power_law_sky,
)
import numpy as np
import matplotlib.pyplot as plt
from astropy import units
import matplotlib as mpl

plt.style.use("./example_style.mplstyle")
mpl.rcParams["text.usetex"] = False
_ = np.seterr(divide="ignore")

Set up simulation settings

Here we set the spherical harmonic bandlimits of the power primary beam, baseline term and the resulting maximum sky resolution we need. We also define a single-channel frequency axis and a fixed CIRS epoch for the sky model.

beam_lmax = 40  # Maximum harmonic degree for power primary beam
baseline_lmax = 50  # Maximum harmonic degree for baseline term
sky_lmax = beam_lmax + baseline_lmax  # Resulting maximum harmonic degree for sky term

# Telescope position in TIRS lat, lon in radians.
telescope_latitude = 0.0
telescope_longitude = 0.0

# Observed frequency in MHz
frequency_MHz = 600.0
frequencies_MHz = np.array([frequency_MHz])
sky_epoch = "J2026.0"

Simulated Sky Containers

We'll use two different simulated skies. A power-law Gaussian random field simulates diffuse emission, and a second container holds a single point source. In both cases the container keeps the frequency axis and CIRS metadata together with the harmonic coefficients.

Power-law sky

Generate a power-law random sky, a Gaussian random field with power spectrum \(\propto \ell^{-2}\).

pl_sky_grid = power_law_sky(sky_lmax, power_law_index=2, seed=123)
pl_sky = SkyModel(
    lmax=sky_lmax,
    frequencies_MHz=frequencies_MHz,
    alm=sht.analysis(pl_sky_grid, return_sparse=False)[None, ...],
    map_unit="K",
    polarisation="I",
    coordinate_basis="CIRS",
    epoch=sky_epoch,
)
_ = pl_sky.as_grid(freq_inds=0, return_shgrid=True).plot(show=False)

Point Source Sky

Generate a sky container with a single point source at CIRS R.A. = \(180^\circ\) and Dec. = \(0^\circ\).

ps_sky = SkyModel(
    lmax=sky_lmax,
    frequencies_MHz=frequencies_MHz,
    alm=harmonic_point_source(180, 0, sky_lmax, return_sparse=False)[None, ...],
    map_unit="K",
    polarisation="I",
    coordinate_basis="CIRS",
    epoch=sky_epoch,
)
_ = ps_sky.as_grid(freq_inds=0, return_shgrid=True).plot(show=False)

Power Primary Beam

Generate a simulated Gaussian primary beam for a telescope at the given latitude pointing at zenith, then wrap it in a TIRSPowerBeam container. As with the sky containers, plotting stays container-native and we only unwrap the coefficients when the core visibility API needs them.

# D_eff gives FWHM = λ/D_eff ≈ 20° at 600 MHz
D_eff = 0.5 / np.radians(20)
power_beam = TIRSPowerBeam.from_gaussian(
    D_eff=D_eff,
    lmax=beam_lmax,
    frequencies_MHz=frequencies_MHz,
    pol_product="XX",
    sky_pol="I",
    latitude=telescope_latitude,
    longitude=telescope_longitude,
    altitude=np.pi / 2,
    azimuth=0,
)
_ = power_beam.as_grid(freq_inds=0, return_shgrid=True).plot(show=False)

Baseline Term

Set the baseline length such that it has equal East and North projections and a length that is band-limited with the set baseline_lmax

bs_pw_mag = plane_wave_mag_from_bandlimit(baseline_lmax) / 2 / np.pi
baseline_length = (frequency_MHz * units.MHz).to(
    "m", equivalencies=units.spectral()
) * bs_pw_mag  # Compute band-limited baseline length
baseline_unit_enu = np.array(
    [1.0, 1.0, 0]
)  # Baseline direction in East-North-Up coordinates.
baseline_unit_enu /= np.linalg.norm(baseline_unit_enu)
baseline_enu = tuple(baseline_unit_enu * baseline_length.value)  # Rescale to length
baseline_enu

(np.float64(1.2782342874738888),
 np.float64(1.2782342874738888),
 np.float64(0.0))

Generate Visibilities and m-modes

Now we'll generate m-modes and visibilities for both methods and compare them. This is the only place we leave the container interface: SingleVisibilityGenerator currently expects raw harmonic coefficient arrays, so we pass sky_model.as_coeff(...) and power_beam.as_coeff(...) explicitly.

For Power-Law Sky

single_visgen = SingleVisibilityGenerator(
    latitude=telescope_latitude,
    longitude=telescope_longitude,
    baseline_enu=baseline_enu,
    frequency_MHz=frequency_MHz,
    power_beam_alm=power_beam.as_coeff(freq_inds=0),
    sky_alm=pl_sky.as_coeff(freq_inds=0),
    method="gaunt",  # Use the Gaunt coefficient method
)

mmodes_gaunt = single_visgen.mmodes()
vis_gaunt = single_visgen.visibilities()

single_visgen.method = "grid"  # Switch to grid based method

mmodes_grid = single_visgen.mmodes()
vis_grid = single_visgen.visibilities()

eras = np.linspace(0, 360, len(vis_gaunt), endpoint=False)
ms = single_visgen.sky_m_values

fig, axs = plt.subplots(2, 2, figsize=(8, 6))

axs[0, 0].plot(eras, vis_gaunt.real)
axs[0, 0].plot(eras, vis_grid.real, ls="--")
axs[0, 0].set_ylabel(r"Re[$\mathcal{V}(\theta_{\mathrm{ERA}})$]")
axs[0, 0].set_xlabel(r"$\theta_{\mathrm{ERA}}$ [deg]")

axs[1, 0].plot(eras, vis_gaunt.imag)
axs[1, 0].plot(eras, vis_grid.imag, ls="--")
axs[1, 0].set_ylabel(r"Im[$\mathcal{V}(\theta_{\mathrm{ERA}})$]")
axs[1, 0].set_xlabel(r"$\theta_{\mathrm{ERA}}$ [deg.]")

axs[0, 1].plot(ms, mmodes_gaunt.real, label="Gaunt")
axs[0, 1].plot(ms, mmodes_grid.real, ls="--", label="Grid")
axs[0, 1].set_ylabel(r"Re[$\mathcal{V}(m)$]")
axs[0, 1].set_xlabel(r"$m$")
axs[0, 1].legend(loc="upper right")

axs[1, 1].plot(ms, mmodes_gaunt.imag)
axs[1, 1].plot(ms, mmodes_grid.imag, ls="--")
axs[1, 1].set_ylabel(r"Im[$\mathcal{V}(m)$]")
axs[1, 1].set_xlabel(r"$m$")

fig.tight_layout()

For Point Source Sky

single_visgen = SingleVisibilityGenerator(
    latitude=telescope_latitude,
    longitude=telescope_longitude,
    baseline_enu=baseline_enu,
    frequency_MHz=frequency_MHz,
    power_beam_alm=power_beam.as_coeff(freq_inds=0),
    sky_alm=ps_sky.as_coeff(freq_inds=0),
    method="gaunt",
)

mmodes_gaunt = single_visgen.mmodes()
vis_gaunt = single_visgen.visibilities()

single_visgen.method = "grid"

mmodes_grid = single_visgen.mmodes()
vis_grid = single_visgen.visibilities()

eras = np.linspace(0, 360, len(vis_gaunt), endpoint=False)
ms = single_visgen.sky_m_values

fig, axs = plt.subplots(2, 2, figsize=(8, 6))

axs[0, 0].plot(eras, vis_gaunt.real)
axs[0, 0].plot(eras, vis_grid.real, ls="--")
axs[0, 0].set_ylabel(r"Re[$\mathcal{V}(\theta_{\mathrm{ERA}})$]")
axs[0, 0].set_xlabel(r"$\theta_{\mathrm{ERA}}$ [deg]")

axs[1, 0].plot(eras, vis_gaunt.imag)
axs[1, 0].plot(eras, vis_grid.imag, ls="--")
axs[1, 0].set_ylabel(r"Im[$\mathcal{V}(\theta_{\mathrm{ERA}})$]")
axs[1, 0].set_xlabel(r"$\theta_{\mathrm{ERA}}$ [deg.]")

axs[0, 1].plot(ms, mmodes_gaunt.real, label="Gaunt")
axs[0, 1].plot(ms, mmodes_grid.real, ls="--", label="Grid")
axs[0, 1].set_ylabel(r"Re[$\mathcal{V}(m)$]")
axs[0, 1].set_xlabel(r"$m$")
axs[0, 1].legend(loc="upper right")

axs[1, 1].plot(ms, mmodes_gaunt.imag)
axs[1, 1].plot(ms, mmodes_grid.imag, ls="--")
axs[1, 1].set_ylabel(r"Im[$\mathcal{V}(m)$]")
axs[1, 1].set_xlabel(r"$m$")

fig.tight_layout()

Split by m

The Gaunt based mode can also compute separately for ranges of sky |m|. Here we compute only the low sky |m| modes.

single_visgen = SingleVisibilityGenerator(
    latitude=telescope_latitude,
    longitude=telescope_longitude,
    baseline_enu=baseline_enu,
    frequency_MHz=frequency_MHz,
    power_beam_alm=power_beam.as_coeff(freq_inds=0),
    sky_alm=ps_sky.as_coeff(freq_inds=0),
    method="gaunt",
    sky_absm_limits=(0, 50),  # Only |m| < 50
)

mmodes_gaunt = single_visgen.mmodes()
vis_gaunt = single_visgen.visibilities()

eras = np.linspace(0, 360, 2 * single_visgen.sky_lmax + 1, endpoint=False)
cut_ms = single_visgen.sky_m_values

fig, axs = plt.subplots(2, 2, figsize=(8, 6))

axs[0, 0].plot(eras, vis_gaunt.real)
axs[0, 0].plot(eras, vis_grid.real, ls="--")
axs[0, 0].set_ylabel(r"Re[$\mathcal{V}(\theta_{\mathrm{ERA}})$]")
axs[0, 0].set_xlabel(r"$\theta_{\mathrm{ERA}}$ [deg]")

axs[1, 0].plot(eras, vis_gaunt.imag)
axs[1, 0].plot(eras, vis_grid.imag, ls="--")
axs[1, 0].set_ylabel(r"Im[$\mathcal{V}(\theta_{\mathrm{ERA}})$]")
axs[1, 0].set_xlabel(r"$\theta_{\mathrm{ERA}}$ [deg.]")

axs[0, 1].plot(cut_ms, mmodes_gaunt.real, label="Gaunt")
axs[0, 1].plot(ms, mmodes_grid.real, ls="--", label="Grid")
axs[0, 1].set_ylabel(r"Re[$\mathcal{V}(m)$]")
axs[0, 1].set_xlabel(r"$m$")
axs[0, 1].legend(loc="upper right")

axs[1, 1].plot(cut_ms, mmodes_gaunt.imag)
axs[1, 1].plot(ms, mmodes_grid.imag, ls="--")
axs[1, 1].set_ylabel(r"Im[$\mathcal{V}(m)$]")
axs[1, 1].set_xlabel(r"$m$")

fig.tight_layout()