Simple Usage
from serval.core import SingleVisibilityGenerator
from serval import sht
from serval.utils import (
pointed_gaussian,
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
mpl.rcParams["text.usetex"] = True
This notebook shows a simple example of computing mmodes and visibilities for a single baseline and frequency channel. Here we use either a direct grid based method or a Gaunt co-efficient based method. Serval in general is designed for re-using cached computations, allowable by the Gaunt methods. In the single visibility case, it is therefore quite inefficient compared to the grid bases method.
Set up simulation settings
Here we set the spherical harmonic bandlimits of the power primary beam, baseline term and the resulting maxmimum sky resolution we need. Additionally some arbitrary meta-data for the telescope being simulated is set here.
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 maxmimum harmonic degree for sky term
# Telescope position in TIRS lat, lon in radians.
telescope_latitude = 0
telescope_longitude = 0
# Observed frequency in MHz
frequency_MHz = 600
Simulated Sky Harmonics
We'll use two different simulated skies. A power-law Gaussian random field, simulating diffuse emission and a sky with a single point source.
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)
pl_sky_alm = sht.analysis(pl_sky_grid, return_sparse=False)
pl_sky_grid.plot()
/home/devin/Dropbox/arcturus/Projects/radio/hirax/systematics/dev/serval/src/serval/utils.py:82: RuntimeWarning: divide by zero encountered in power
power = ells ** (-power_law_index)
/home/devin/miniforge3/envs/serval_dev/lib/python3.11/site-packages/pyshtools/shclasses/shgrid.py:1464: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
fig.show()
(<Figure size 640x768 with 2 Axes>,
array([<Axes: title={'center': 'Real component'}, xlabel='Longitude', ylabel='Latitude'>,
<Axes: title={'center': 'Imaginary component'}, xlabel='Longitude', ylabel='Latitude'>],
dtype=object))
Point Source Sky
Generate a sky with a single point source at CIRS R.A. = \(180^\circ\) and Dec. = \(0^\circ\)
ps_sky_alm = harmonic_point_source(180, 0, sky_lmax).todense()
ps_sky_grid = sht.array_synthesis(ps_sky_alm)
ps_sky_grid.plot()
/home/devin/miniforge3/envs/serval_dev/lib/python3.11/site-packages/pyshtools/shclasses/shgrid.py:1464: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
fig.show()
(<Figure size 640x768 with 2 Axes>,
array([<Axes: title={'center': 'Real component'}, xlabel='Longitude', ylabel='Latitude'>,
<Axes: title={'center': 'Imaginary component'}, xlabel='Longitude', ylabel='Latitude'>],
dtype=object))
Power Primary Beam
Generate a simulated Gaussian primary beam for a telesope at the given latitude pointing at zenith.
beam = pointed_gaussian(
altitude=np.pi / 2, # Pointing at zenith
azimuth=0,
latitude=telescope_latitude,
longitude=telescope_longitude,
fwhm=np.radians(20), # 20 deg wide primary beam
template=sht.grid_template(beam_lmax),
)
beam_alm = sht.analysis(beam, return_sparse=False)
beam.plot()
/home/devin/miniforge3/envs/serval_dev/lib/python3.11/site-packages/pyshtools/shclasses/shgrid.py:1464: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
fig.show()
(<Figure size 640x768 with 2 Axes>,
array([<Axes: title={'center': 'Real component'}, xlabel='Longitude', ylabel='Latitude'>,
<Axes: title={'center': 'Imaginary component'}, xlabel='Longitude', ylabel='Latitude'>],
dtype=object))
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
(1.2782342874738888, 1.2782342874738888, 0.0)
Generate Visibilities and m-modes
Now we'll generate mmodes and visibilities for both methods and compare them.
For Power-Law Sky
single_visgen = SingleVisibilityGenerator(
latitude=telescope_latitude,
longitude=telescope_longitude,
baseline_enu=baseline_enu,
frequency_MHz=frequency_MHz,
power_beam_alm=beam_alm,
sky_alm=pl_sky_alm,
method="gaunt", # Use the Gaunt co-efficient 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_{\rm ERA})$]")
axs[0, 0].set_xlabel(r"$\theta_{\rm 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_{\rm ERA})$]")
axs[1, 0].set_xlabel(r"$\theta_{\rm 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=beam_alm,
sky_alm=ps_sky_alm,
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_{\rm ERA})$]")
axs[0, 0].set_xlabel(r"$\theta_{\rm 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_{\rm ERA})$]")
axs[1, 0].set_xlabel(r"$\theta_{\rm 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=beam_alm,
sky_alm=ps_sky_alm,
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_{\rm ERA})$]")
axs[0, 0].set_xlabel(r"$\theta_{\rm 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_{\rm ERA})$]")
axs[1, 0].set_xlabel(r"$\theta_{\rm 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()
