"""Telescope sensitivity models and simulated interferometric UV coverage.
The module combines three related pieces:
* :class:`Telescope` subclasses describe an array layout, location, primary
beam, useful baseline range, and frequency-dependent SEFD;
* :class:`TelescopeSimu` projects the physical baselines through an observation;
* :class:`SimuGridded` stores the resulting UV weights and derives thermal
noise or a matching power-spectrum estimator.
Sensitivity models can be used independently of a UV simulation::
import numpy as np
from ps_eor import obssimu
freqs = np.arange(50, 80, .2) * 1e6
telescope = obssimu.Telescope.from_name('nenufar')
stokes_i_sefd = telescope.get_i_sefd(freqs)
For an array with an available station layout, configure the declination, hour
angle range, and time resolution before choosing a gridding strategy::
telescope = obssimu.DEx(n_antenna_side=16, sep_antenna=6)
observation = obssimu.TelescopeSimu(
telescope,
freqs,
dec_deg=-27,
hal=-1,
har=1,
timeres=60,
)
coverage = observation.image_gridding(
fov_deg=telescope.fov,
min_weight=1,
)
noise_std = coverage.get_noise_std_cube(total_time_sec=100 * 3600)
The same :class:`SimuGridded` object provides a power-spectrum estimator whose
frequency range, UV limits, primary beam, and default weights match the
simulation. The noise-standard-deviation cube can be passed directly to it to
obtain the expected thermal-noise power::
ps_gen = coverage.get_ps_gen(
filter_kpar_min=.05,
filter_wedge_theta=0,
)
kbins = np.logspace(np.log10(ps_gen.kmin), np.log10(.5), 8)
noise_spectra = ps_gen.get_all(kbins, noise_std)
noise_spectra.ps.plot()
noise_spectra.ps2d.plot()
noise_spectra.ps3d.plot()
The ``.data`` attribute of each result is the expected thermal-noise power,
not the sensitivity. The corresponding one-sigma sensitivity is stored in
``.err``; for example, use ``noise_spectra.ps3d.err`` for the spherical
sensitivity. A random noise realization is only needed when testing a
realization-dependent analysis and can be drawn with
``noise_std.generate_noise_cube()``.
Frequencies are in Hz, physical coordinates in metres, UV coordinates in
wavelengths, angles in degrees unless a parameter explicitly ends in
``_rad``, hour angles in hours, and times in seconds.
"""
import itertools
import os
from enum import Enum, auto
from importlib.resources import files as resource_files
from typing import ClassVar
import astropy.constants as const
import numpy as np
import scipy.interpolate
from astropy import units
from astropy.coordinates import EarthLocation
from fast_histogram import histogram2d
INSTRU_DIR = str(resource_files('ps_eor').joinpath('data', 'instruments'))
[docs]
class CoordinateType(Enum):
"""Enumeration of coordinate systems used to interpret station positions."""
ENU = auto() # East-North-Up (local topocentric)
ECEF = auto() # Earth-Centered, Earth-Fixed (global geodetic)
XYZ = auto() # Rotated Cartesian XYZ used for UVW conversion
[docs]
def enu_to_ecef(location, enu):
"""Convert local east-north-up coordinates to Earth-fixed coordinates.
Args:
location (astropy.coordinates.EarthLocation): Reference geodetic location.
enu (np.ndarray): ENU coordinates, shape (N, 3), in meters.
Returns:
np.ndarray: ECEF coordinates, shape (N, 3), in meters.
"""
e, n, u = np.hsplit(enu, 3)
lon = location.geodetic[0].to(units.rad).value
lat = location.geodetic[1].to(units.rad).value
alt = location.geodetic[2].to(units.m).value
x, y, z = lla_to_ecef(lat, lon, alt)
sin_lat, cos_lat = np.sin(lat), np.cos(lat)
sin_lon, cos_lon = np.sin(lon), np.cos(lon)
X = x - sin_lon * e - sin_lat * cos_lon * n + cos_lat * cos_lon * u
Y = y + cos_lon * e - sin_lat * sin_lon * n + cos_lat * sin_lon * u
Z = z + cos_lat * n + sin_lat * u
return np.hstack([X, Y, Z])
[docs]
def ecef_to_xyz_matrix(long_rad):
"""Return the rotation matrix from ECEF to interferometric XYZ.
Args:
long_rad (float): Longitude in radians.
Returns:
np.ndarray: 3×3 rotation matrix.
""" # noqa: RUF002
return np.array([
[ np.cos(long_rad), np.sin(long_rad), 0],
[-np.sin(long_rad), np.cos(long_rad), 0],
[0, 0, 1]
])
[docs]
def lla_to_ecef(lat, lon, alt):
"""Convert WGS84 geodetic coordinates to ECEF coordinates.
Args:
lat (float): Geodetic latitude in radians.
lon (float): Geodetic longitude in radians.
alt (float): Height above ellipsoid in meters.
Returns:
tuple of float: Cartesian ECEF coordinates (x, y, z) in meters.
"""
WGS84_a = 6378137.0
WGS84_b = 6356752.31424518
N = WGS84_a**2 / np.sqrt(
WGS84_a**2 * np.cos(lat)**2 + WGS84_b**2 * np.sin(lat)**2
)
x = (N + alt) * np.cos(lat) * np.cos(lon)
y = (N + alt) * np.cos(lat) * np.sin(lon)
z = ((WGS84_b**2 / WGS84_a**2) * N + alt) * np.sin(lat)
return x, y, z
[docs]
def xyz_to_uvw_matrix(ha_rad, dec_rad):
"""Return the rotation matrix from interferometric XYZ to UVW.
Args:
ha_rad (float): Hour angle in radians.
dec_rad (float): Declination in radians.
Returns:
np.ndarray: 3x3 transformation matrix from XYZ to UVW.
"""
return np.array([
[ np.sin(ha_rad), np.cos(ha_rad), 0.0],
[-np.sin(dec_rad) * np.cos(ha_rad), np.sin(dec_rad) * np.sin(ha_rad), np.cos(dec_rad)],
[ np.cos(dec_rad) * np.cos(ha_rad), -np.cos(dec_rad) * np.sin(ha_rad), np.sin(dec_rad)]
])
[docs]
class Telescope:
"""Base description of a radio telescope used by the simulators.
Subclasses provide their station layout and single-polarization SEFD.
Attributes:
name: Identifier accepted by :meth:`from_name`.
pb_name: Primary-beam identifier accepted by
:meth:`ps_eor.datacube.PrimaryBeam.from_name`.
n_elements_per_stations: Number of layout elements forming one station.
only_drift_mode: Whether the model is intended only for drift scans.
redundant_array: Whether to use the redundant-baseline gridding path.
redundant_baselines: Representative redundant baseline lengths in
metres.
umin, umax: Default baseline limits in wavelengths.
coord_type: Coordinate system of the station positions.
location: Array reference location.
"""
name = 'none'
pb_name = name
n_elements_per_stations = 1
only_drift_mode = False
redundant_array = False
redundant_baselines: ClassVar[list] = []
umin = 0
umax = 10000
coord_type = CoordinateType.ENU
location = EarthLocation(lon=6.8670 * units.deg, lat=52.9088 * units.deg, height=15.0 * units.m)
[docs]
def get_stat_pos_file(self):
"""The station-position file supplied by the telescope model.
Returns:
str: Path to an ``(n_elements, 3)`` text array in metres.
"""
[docs]
def get_stat_pos(self):
"""Load station positions from :meth:`get_stat_pos_file`.
Returns:
ndarray: Station positions shaped ``(n_elements, 3)`` in metres.
"""
return np.loadtxt(self.get_stat_pos_file())
[docs]
def get_sefd(self, freq):
"""Return the single-polarization SEFD in Jy.
Args:
freq (np.ndarray or float): Frequency in Hz.
Returns:
ndarray or float: Single-polarization SEFD in Jy.
"""
[docs]
def get_i_sefd(self, freq):
"""Return the Stokes-I SEFD in Jy.
Args:
freq (np.ndarray or float): Frequency in Hz.
Returns:
ndarray or float: Stokes-I SEFD in Jy.
"""
return 1 / np.sqrt(2) * self.get_sefd(freq)
[docs]
def sky_temperature(self, freq, tsys_sky=60, temp_power_law_index=2.55):
"""Estimate the sky temperature from a wavelength power law.
Args:
freq (np.ndarray or float): Frequency in Hz.
tsys_sky (float): Sky temperature at one-metre wavelength, in K.
temp_power_law_index (float): Wavelength power-law index.
Returns:
np.ndarray or float: Estimated sky temperature in Kelvin.
"""
lamb = const.c.value / freq
return tsys_sky * lamb ** temp_power_law_index
[docs]
def get_dipole_aeff(self, freq, distance_between_dipole):
"""Return the effective area of one dipole in square metres.
Args:
freq (np.ndarray or float): Frequency in Hz.
distance_between_dipole (float): Physical spacing between dipoles in meters.
Returns:
np.ndarray or float: Effective area in m².
"""
lamb = const.c.value / freq
return np.min(
[lamb ** 2 / 3, np.ones_like(lamb) * np.pi * distance_between_dipole ** 2 / 4.],
axis=0)
[docs]
def get_dish_aeff(self, freq, diameter, efficiency):
"""Return the effective area of one dish in square metres.
Args:
freq (np.ndarray or float): Frequency in Hz.
diameter (float): Diameter of the dish in meters.
efficiency (float): Aperture efficiency.
Returns:
np.ndarray or float: Effective area in m².
"""
lamb = const.c.value / freq
return lamb ** 2 / (4 * np.pi) * efficiency * (np.pi * diameter / lamb) ** 2
[docs]
@staticmethod
def from_name(name):
"""Create a registered telescope model from its name.
Args:
name (str): Telescope name, for example ``"ska_low"``.
Returns:
Telescope: Instance of the corresponding Telescope subclass.
Raises:
ValueError: If no matching telescope subclass is found.
"""
klasses = Telescope.__subclasses__()
[klasses.extend(k.__subclasses__()) for k in klasses[:]]
for klass in klasses:
if hasattr(klass, 'name') and klass.name == name:
return klass()
raise ValueError(f'No telescope with name: {name}')
[docs]
class DEx(Telescope):
"""Conceptual lunar dipole array with square grid layout."""
name = 'dex'
umin = 0.5
umax = 25
fov = 120
du = 1
only_drift_mode = True
coord_type = CoordinateType.ENU
location = EarthLocation(lon=116.67081524 * units.deg, lat=-26.70122102627586 * units.deg,
height=0.0 * units.m)
def __init__(self, n_antenna_side=32, sep_antenna=6):
"""Configure a square array.
Args:
n_antenna_side: Number of antennas along each side.
sep_antenna: Antenna separation in metres.
"""
Telescope.__init__(self)
self.sep_antenna = sep_antenna
self.n_antenna_side = n_antenna_side
@property
def pb_name(self):
"""Dipole beam name, following the configured antenna separation.
:meth:`get_sefd` models one element as an aperture of diameter
``sep_antenna`` (see :meth:`Telescope.get_dipole_aeff`), so the beam
uses the same aperture. The tapering matches the other dipole arrays
here (``a12_lba``, ``ovro_lwa``). Unlike every other telescope the
element geometry is a constructor argument, so this cannot be a fixed
class attribute.
"""
return f'ant_{self.sep_antenna}_1.1_gaussian'
[docs]
def get_stat_pos(self):
"""Generate ENU antenna positions for the square layout.
Returns:
ndarray: Station positions shaped ``(n_antenna_side**2, 3)`` in
metres.
"""
grid_indices = np.arange(0, self.n_antenna_side)
p_east, p_north = np.meshgrid(self.sep_antenna * grid_indices,
self.sep_antenna * grid_indices)
east = p_east.flatten()
north = p_north.flatten()
up = np.zeros_like(east)
return np.vstack((east, north, up)).T
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The single-polarization dipole SEFD in Jy."""
tsys = self.sky_temperature(freq, tsys_sky)
a_eff = self.get_dipole_aeff(freq, self.sep_antenna)
return 2 * const.k_B.value / a_eff * 1e26 * tsys
[docs]
class SkaLow(Telescope):
"""SKA-Low Phase 1 array (AA4)."""
name = 'ska_low'
umin = 30
umax = 250
fov = 3
du = 8
pb_name = name
coord_type = CoordinateType.ENU
location = EarthLocation.from_geodetic(116.7644482, -26.82472208, 365.0, ellipsoid="WGS84")
[docs]
def get_stat_pos_file(self):
"""Return the packaged SKA-Low AA4 ENU station layout."""
return os.path.join(INSTRU_DIR, 'ska1_low_enu_statpos.data')
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The interpolated single-polarization SKA-Low SEFD in Jy."""
# SKA LFAA station specification, arXiv:2003.12744, p. 22.
freqs_spec = np.array([50, 80, 110, 140, 160, 220]) * 1e6
a_eff_over_tsys_spec = 1 * np.array([0.14, 0.46, 1.04, 1.15, 1.2, 1.2])
def t_sky_fct(freqs): return tsys_sky * (3e8 / freqs) ** 2.55
a_eff_fct = scipy.interpolate.interp1d(
freqs_spec, a_eff_over_tsys_spec * t_sky_fct(freqs_spec), kind='slinear',
bounds_error=False, fill_value='extrapolate')
return 2 * const.k_B.value * 1e26 * t_sky_fct(freq) / a_eff_fct(freq)
[docs]
class SkaLowAAstar(SkaLow):
"""SKA-Low AA* layout."""
name = 'ska_low_aastar'
[docs]
def get_stat_pos_file(self):
"""Return the packaged SKA-Low AA* ENU station layout."""
return os.path.join(INSTRU_DIR, 'ska1_low_aastar_enu_statpos.data')
[docs]
class SkaLowAA2(SkaLow):
"""SKA-Low AA2 layout."""
name = 'ska_low_aa2'
[docs]
def get_stat_pos_file(self):
"""Return the packaged SKA-Low AA2 ENU station layout."""
return os.path.join(INSTRU_DIR, 'ska1_low_aa2_enu_statpos.data')
[docs]
class LofarHBA(Telescope):
"""LOFAR High-Band Antenna (HBA) core array."""
name = 'lofar_hba'
umin = 50
umax = 250
fov = 4
du = 8
n_elements_per_stations = 2
pb_name = name
coord_type = CoordinateType.ECEF
location = EarthLocation(lon=6.8670 * units.deg, lat=52.9088 * units.deg, height=15.0 * units.m)
[docs]
def get_stat_pos_file(self):
"""Return the packaged LOFAR ECEF station layout."""
return os.path.join(INSTRU_DIR, 'lofar_statpos.data')
[docs]
def get_sefd(self, freq):
"""The representative single-polarization LOFAR-HBA SEFD in Jy."""
# Typical observed NCP sensitivity between 130 and 160 MHz.
return 4000
[docs]
class A12HBA(Telescope):
"""AARTFAAC-12 LOFAR HBA subarray with 48-element stations."""
name = 'a12_hba'
umin = 10
umax = 200
fov = 24
du = 2
pb_name = name
n_elements_per_stations = 48
coord_type = CoordinateType.ECEF
location = EarthLocation(lon=6.8670 * units.deg, lat=52.9088 * units.deg, height=15.0 * units.m)
[docs]
def get_stat_pos_file(self):
"""Return the packaged AARTFAAC-12 HBA ECEF element layout."""
return os.path.join(INSTRU_DIR, 'aartfaac_a12_hba_statpos.data')
[docs]
def get_sefd(self, freq):
"""The scaled single-polarization AARTFAAC-12 HBA SEFD in Jy."""
# Scale the observed LOFAR-HBA station SEFD to one of its 24 tiles.
return 4000 * 24
[docs]
class A12LBA(Telescope):
"""AARTFAAC-12 LOFAR LBA subarray with wide FoV and drift scan."""
name = 'a12_lba'
umin = 20
umax = 40
fov = 120
du = 1
pb_name = 'ant_1.9_1.1_gaussian'
n_elements_per_stations = 48
only_drift_mode = True
coord_type = CoordinateType.ECEF
location = EarthLocation(lon=6.8670 * units.deg, lat=52.9088 * units.deg, height=15.0 * units.m)
[docs]
def get_stat_pos_file(self):
"""Return the packaged AARTFAAC-12 LBA ECEF element layout."""
return os.path.join(INSTRU_DIR, 'aartfaac_a12_lba_statpos.data')
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The sky-noise dominated single-polarization LBA SEFD in Jy."""
distance_between_dipole = 7
tsys = self.sky_temperature(freq, tsys_sky)
a_eff = self.get_dipole_aeff(freq, distance_between_dipole)
return 2 * const.k_B.value / a_eff * 1e26 * tsys
[docs]
class MWA1(Telescope):
"""First-phase MWA core layout."""
name = 'mwa1'
umin = 18
umax = 80
fov = 30
du = 2
pb_name = 'ant_4_1.05_gaussian'
only_drift_mode = False
coord_type = CoordinateType.ECEF
location = EarthLocation(lon=116.67044463048276 * units.deg, lat=-26.70122102627586 * units.deg,
height=377.8 * units.m)
[docs]
def get_stat_pos_file(self):
"""Return the packaged MWA Phase-I ECEF tile layout."""
return os.path.join(INSTRU_DIR, 'mwa_rev1_statpos.data')
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The sky-noise dominated single-polarization MWA SEFD in Jy."""
distance_between_dipole = 1.1
n_dipole_per_stations = 4 * 4
tsys = self.sky_temperature(freq, tsys_sky)
a_eff = n_dipole_per_stations * self.get_dipole_aeff(freq, distance_between_dipole)
return 2 * const.k_B.value / a_eff * 1e26 * tsys
[docs]
class HERA(Telescope):
"""HERA telescope core array with redundant hexagonal layout."""
name = 'hera'
pb_name = 'ant_14_1.1_gaussian'
only_drift_mode = True
redundant_array = True
redundant_baselines = np.array([14.6, 25.28794179, 29.2, 38.62796914, 43.8,
50.57588358, 52.64104862, 58.4, 63.63992458,
66.90560515, 73., 75.86382537, 77.25593828,
81.2893597, 87.6])
umin = 4
umax = 200
du = 1
coord_type = CoordinateType.ENU
hera_location = EarthLocation(lon=21.43 * units.deg, lat=-30.72 * units.deg,
height=1050.0 * units.m)
def __init__(self, hex_num=11, split_core=True, sep=14.6):
"""Configure a hexagonal HERA layout.
Args:
hex_num: Number of antennas along the central hexagon axis.
split_core: Offset the three core sectors.
sep: Nearest-neighbour separation in metres.
"""
self.hex_num = hex_num
self.split_core = split_core
self.sep = sep
[docs]
def get_stat_pos(self):
"""Generate ENU antenna positions for the hexagonal layout.
Returns:
ndarray: Antenna positions shaped ``(n_antennas, 3)`` in metres.
"""
# Layout construction adapted from HERA-Team/hera_sim.
positions = []
for row in range(self.hex_num - 1, -self.hex_num + self.split_core, -1):
# A split core omits the central row.
for col in range(2 * self.hex_num - abs(row) - 1):
x_pos = self.sep * ((2 - (2 * self.hex_num - abs(row))) / 2 + col)
y_pos = row * self.sep * np.sqrt(3) / 2
positions.append([x_pos, y_pos, 0])
# Hexagonal-grid basis vectors, normalized to ``self.sep``.
up_right = self.sep * np.asarray([0.5, np.sqrt(3) / 2, 0])
up_left = self.sep * np.asarray([-0.5, np.sqrt(3) / 2, 0])
if self.split_core:
new_pos = []
for pos in positions:
# Move each sector away from the core centre.
theta = np.arctan2(pos[1], pos[0])
if pos[0] == 0 and pos[1] == 0:
new_pos.append(pos)
elif -np.pi / 3 < theta < np.pi / 3:
new_pos.append(np.asarray(pos) + (up_right + up_left) / 3)
elif np.pi / 3 <= theta < np.pi:
new_pos.append(np.asarray(pos) + up_left - (up_right + up_left) / 3)
else:
new_pos.append(pos)
positions = new_pos
return np.array(positions)
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The single-polarization dish SEFD in Jy."""
d = 14
eff = 0.78
trxc = 100
lamb = const.c.value / freq
a_eff = self.get_dish_aeff(freq, d, eff)
tsys = tsys_sky * lamb ** 2.55 + trxc
return 2 * const.k_B.value / a_eff * 1e26 * tsys
[docs]
class HERA56(HERA):
"""HERA sub-array with 56 antennas."""
name = 'hera_56'
def __init__(self):
HERA.__init__(self, 5)
[docs]
class HERA120(HERA):
"""HERA sub-array with 120 antennas."""
name = 'hera_120'
def __init__(self):
HERA.__init__(self, 7)
[docs]
class HERA208(HERA):
"""HERA sub-array with 208 antennas."""
name = 'hera_208'
def __init__(self):
HERA.__init__(self, 9)
[docs]
class HERA320(HERA):
"""HERA full array with 320 antennas."""
name = 'hera_320'
def __init__(self):
HERA.__init__(self, 11)
[docs]
class NenuFAR(Telescope):
"""Full NenuFAR array in Nançay (France)."""
name = 'nenufar'
pb_name = 'nenufar'
umin = 6
umax = 60
fov = 16
du = 4
coord_type = CoordinateType.ECEF
location = EarthLocation(lon=2.192400 * units.deg, lat=47.376511 * units.deg,
height=182.096 * units.m)
[docs]
def get_stat_pos_file(self):
"""Return the packaged full-NenuFAR ECEF mini-array layout."""
return os.path.join(INSTRU_DIR, 'nenufar_full_statpos.data')
[docs]
def inst_temperature(self, freq):
"""Return the measured receiver temperature in K.
The interpolation data come from ``nenupy.instru``.
Args:
freq: Frequency in Hz.
Returns:
ndarray or float: Receiver temperature in K.
"""
lna_sky = np.array([
5.0965, 2.3284, 1.0268, 0.4399, 0.2113, 0.1190, 0.0822, 0.0686,
0.0656, 0.0683, 0.0728, 0.0770, 0.0795, 0.0799, 0.0783, 0.0751,
0.0710, 0.0667, 0.0629, 0.0610, 0.0614, 0.0630, 0.0651, 0.0672,
0.0694, 0.0714, 0.0728, 0.0739, 0.0751, 0.0769, 0.0797, 0.0837,
0.0889, 0.0952, 0.1027, 0.1114, 0.1212, 0.1318, 0.1434, 0.1562,
0.1700, 0.1841, 0.1971, 0.2072, 0.2135, 0.2168, 0.2175, 0.2159,
0.2121, 0.2070, 0.2022, 0.1985, 0.1974, 0.2001, 0.2063, 0.2148,
0.2246, 0.2348, 0.2462, 0.2600, 0.2783, 0.3040, 0.3390, 0.3846,
0.4425, 0.5167, 0.6183, 0.7689, 1.0086, 1.4042, 2.0732
])
lna_freqs = (np.arange(71) + 15) * 1e6
return self.sky_temperature(
freq) * scipy.interpolate.interp1d(lna_freqs, lna_sky, bounds_error=False,
fill_value='extrapolate')(freq)
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The modeled single-polarization NenuFAR SEFD in Jy."""
distance_between_dipole = 5.5
n_dipole_per_stations = 19
tsys = self.sky_temperature(freq, tsys_sky) + self.inst_temperature(freq)
a_eff = n_dipole_per_stations * self.get_dipole_aeff(freq, distance_between_dipole)
return 2 * const.k_B.value / a_eff * 1e26 * tsys
[docs]
class NenuFAR80(NenuFAR):
"""Subset of NenuFAR using 80 mini-arrays."""
name = 'nenufar_80'
[docs]
def get_stat_pos_file(self):
"""Return the packaged 80-mini-array NenuFAR ECEF layout."""
return os.path.join(INSTRU_DIR, 'nenufar80_statpos.data')
[docs]
class OVROLWA(Telescope):
"""OVRO-LWA dipole array in California with large FoV."""
name = 'ovro_lwa'
umin = 2
umax = 60
fov = 120
du = 1
pb_name = 'ant_1.9_1.1_gaussian'
only_drift_mode = True
coord_type = CoordinateType.ENU
location = EarthLocation(lon=118.275 * units.deg, lat=37.2337 * units.deg,
height=1222.0 * units.m)
[docs]
def get_stat_pos_file(self):
"""Return the packaged OVRO-LWA ENU station layout."""
return os.path.join(INSTRU_DIR, 'ovro-lwa_enu_statpos.data')
[docs]
def get_sefd(self, freq, tsys_sky=60):
"""The sky-noise dominated single-polarization OVRO-LWA SEFD in Jy."""
distance_between_dipole = 5
tsys = self.sky_temperature(freq, tsys_sky)
a_eff = self.get_dipole_aeff(freq, distance_between_dipole)
return 2 * const.k_B.value / a_eff * 1e26 * tsys
[docs]
class TelescopeSimu:
"""Simulate a telescope's UV coverage over an observation."""
def __init__(self, telescop: Telescope, freqs, dec_deg, hal, har, umin=None, umax=None,
timeres=100, remove_intra_baselines=True):
"""Configure an observation and its baseline selection.
Args:
telescop (Telescope): Array layout and sensitivity model.
freqs (array-like): Observing frequencies in Hz.
dec_deg (float): Pointing declination in degrees.
hal (float): Start hour angle in hours.
har (float): End hour angle in hours.
umin (float): Minimum baseline length in wavelengths.
umax (float): Maximum baseline length in wavelengths.
timeres (float): Time sampling in seconds.
remove_intra_baselines (bool): Exclude baselines within a station.
"""
self.telescop = telescop
self.freqs = freqs
self.dec_deg = dec_deg
self.hal = hal
self.har = har
self.umin = umin
if self.umin is None:
self.umin = telescop.umin
self.umax = umax
if self.umax is None:
self.umax = telescop.umax
self.timeres = timeres
self.remove_intra_baselines = remove_intra_baselines
[docs]
@staticmethod
def from_dict(d, freqs):
"""Build a simulation from its ``PE*`` metadata fields.
Args:
d: Mapping containing the fields returned by :meth:`to_dict`.
freqs: Observing frequencies in Hz; these are not stored in the
metadata mapping.
Returns:
TelescopeSimu: The reconstructed observation.
"""
def get_d_value(name, default=None):
if default is None and name not in d:
raise ValueError(f'{name} missing to initialize TelescopeSimu')
return d.get(name, default)
instru = get_d_value('PEINSTRU')
dec_deg = get_d_value('PEOBSDEC')
hal = get_d_value('PEOBSHAL')
har = get_d_value('PEOBSHAR')
timeres = get_d_value('PEOBSRES')
remove_intra_baselines = get_d_value('PEREMINT')
telescop = Telescope.from_name(instru)
umin = get_d_value('PEOBSUMI', telescop.umin)
umax = get_d_value('PEOBSUMA', telescop.umax)
return TelescopeSimu(telescop, freqs, dec_deg, hal, har, umin=umin, umax=umax,
timeres=timeres, remove_intra_baselines=remove_intra_baselines)
[docs]
def to_dict(self):
"""Serialize the observation as ``PE*`` metadata fields.
Returns:
dict: Telescope name, pointing, hour-angle range, baseline limits,
time resolution, and intra-station selection.
"""
return {
'PEINSTRU': self.telescop.name, 'PEOBSDEC': self.dec_deg, 'PEOBSHAL': self.hal,
'PEOBSHAR': self.har, 'PEOBSUMI': self.umin, 'PEOBSUMA': self.umax,
'PEOBSRES': self.timeres, 'PEREMINT': self.remove_intra_baselines}
[docs]
def get_XYZ_positions(self):
"""Station positions in interferometric XYZ coordinates.
Input positions are converted according to :attr:`Telescope.coord_type`
and the array reference location.
Returns:
ndarray: Positions shaped ``(n_elements, 3)`` in metres.
"""
statpos = self.telescop.get_stat_pos()
coord_type = self.telescop.coord_type
location = self.telescop.location
long_rad = location.lon.to(units.rad).value
if coord_type == CoordinateType.ENU:
ecef = enu_to_ecef(location, statpos)
R = ecef_to_xyz_matrix(long_rad)
return ecef @ R.T
elif coord_type == CoordinateType.ECEF:
R = ecef_to_xyz_matrix(long_rad)
return statpos @ R.T
elif coord_type == CoordinateType.XYZ:
return statpos
else:
raise ValueError(f"Unsupported coordinate type: {coord_type}")
[docs]
def simu_uv(self, include_conj=True):
"""Simulate baseline coordinates over the observation.
Physical baselines are sampled from ``hal`` up to, but excluding,
``har`` at ``timeres`` intervals. A baseline is retained if it enters
the requested wavelength range at any supplied frequency.
Args:
include_conj: Include the conjugate ``(-u, -v, w)`` samples.
Returns:
tuple of ndarray: Flattened ``(u, v, w)`` coordinates in metres,
containing every selected baseline and time sample.
"""
from ps_eor import psutil
def m2a(m): return np.squeeze(np.asarray(m))
lambs = const.c.value / self.freqs
umin_meter = (self.umin * lambs).min()
umax_meter = (self.umax * lambs).max()
timev = np.arange(self.hal * 3600, self.har * 3600, self.timeres)
statpos = self.telescop.get_stat_pos()
nstat = statpos.shape[0]
print('Simulating UV coverage ...')
# Each unordered antenna pair contributes one physical baseline.
stncom = np.array(list(itertools.combinations(np.arange(0, nstat), 2)))
print(f'Number of elements: {nstat}')
print(f'Number of baselines: {stncom.shape[0]}')
if self.remove_intra_baselines and self.telescop.n_elements_per_stations > 1:
n_stations = nstat // self.telescop.n_elements_per_stations
station_id = np.repeat(np.arange(n_stations), self.telescop.n_elements_per_stations)
stncom_stations = np.array(list(itertools.combinations(station_id, 2)))
idx = np.array([a == b for a, b, in stncom_stations]).astype(bool)
stncom = stncom[~idx]
print(f'Discarding {idx.sum()} intra-baselines')
b1, b2 = zip(*stncom, strict=False)
uu = []
vv = []
ww = []
pr = psutil.progress_report(len(timev))
for i, tt in enumerate(timev):
pr(i)
ha_rad = (tt / 3600.) * (15. / 180) * np.pi
dec_rad = self.dec_deg * (np.pi / 180)
XYZ = self.get_XYZ_positions()
R = xyz_to_uvw_matrix(ha_rad, dec_rad)
UVW = XYZ @ R.T
bu = m2a(UVW[b1, 0] - UVW[b2, 0])
bv = m2a(UVW[b1, 1] - UVW[b2, 1])
bw = m2a(UVW[b1, 2] - UVW[b2, 2])
ru = np.sqrt(bu ** 2 + bv ** 2)
idx = (ru > umin_meter) & (ru < umax_meter)
uu.extend(bu[idx])
vv.extend(bv[idx])
ww.extend(bw[idx])
if include_conj:
uu.extend(- bu[idx])
vv.extend(- bv[idx])
ww.extend(bw[idx])
return np.array(uu), np.array(vv), np.array(ww)
[docs]
def redundant_gridding(self, max_distance=1):
"""Cluster physically redundant baselines.
This path keeps coordinates in metres so a redundant group represents
the same physical baseline across the band.
Args:
max_distance: DBSCAN clustering distance in metres.
Returns:
SimuGridded: Metre-coordinate weights, with one column per
redundant baseline group.
"""
import sklearn.cluster
from ps_eor import datacube, psutil
uu_meter, vv_meter, _ = self.simu_uv()
X = np.array([uu_meter, vv_meter]).T
c = sklearn.cluster.DBSCAN(eps=max_distance, min_samples=1)
c.fit(X)
_c_id, idx, counts = np.unique(c.labels_, return_index=True, return_counts=True)
uu_meter_grid = uu_meter[idx]
vv_meter_grid = vv_meter[idx]
meta = datacube.ImageMetaData.from_res(0.01, (100, 100))
meta.wcs.wcs.cdelt[2] = psutil.robust_freq_width(self.freqs)
meta.set('PEINTTIM', self.timeres)
meta.set('PETOTTIM', (self.har - self.hal) * 3600)
w_cube = datacube.CartWeightsCubeMeter(np.repeat(counts[None, :], len(self.freqs), 0),
uu_meter_grid, vv_meter_grid, self.freqs, meta)
return SimuGridded(w_cube, self)
[docs]
def image_gridding(self, fov_deg, oversampling_factor=4, min_weight=10, win_fct=None):
"""Bin simulated baselines onto a regular UV grid.
``win_fct`` is recorded in the output metadata for later noise and
power-spectrum corrections; it is not applied to the weights.
Args:
fov_deg: Image field of view in degrees. It sets the UV-cell size.
oversampling_factor: Number of image pixels per finest requested
angular scale.
min_weight: Keep cells with at least this many samples at every
frequency.
win_fct: Optional :class:`~ps_eor.datacube.WindowFunction` to record
in the output metadata.
Returns:
SimuGridded: Frequency-dependent counts on the selected UV cells.
"""
from ps_eor import datacube, psutil
uu_meter, vv_meter, _ = self.simu_uv()
du = 1 / np.radians(fov_deg)
res = 1 / (oversampling_factor * self.umax)
n_u = int(np.ceil(1 / (res * du)))
shape = (n_u, n_u)
g_uu, g_vv = psutil.get_uv_grid(shape, res)
ranges = [g_uu.min() - du / 2, g_uu.max() + du / 2]
print('Gridding UV coverage ...')
weights = []
pr = psutil.progress_report(len(self.freqs))
for i, lamb in enumerate(const.c.value / self.freqs):
pr(i)
w = histogram2d(uu_meter / lamb, vv_meter / lamb, bins=n_u, range=[ranges] * 2)
weights.append(w)
weights = np.array(weights)
weights = weights.reshape(len(self.freqs), -1)
g_uu = g_uu.flatten()
g_vv = g_vv.flatten()
ru = np.sqrt(g_uu ** 2 + g_vv ** 2)
idx = (weights.min(axis=0) >= min_weight) & (ru >= self.umin) & (ru <= self.umax)
weights = weights[:, idx]
g_uu = g_uu[idx]
g_vv = g_vv[idx]
meta = datacube.ImageMetaData.from_res(res, shape)
meta.wcs.wcs.cdelt[2] = psutil.robust_freq_width(self.freqs)
meta.set('PEINTTIM', self.timeres)
meta.set('PETOTTIM', (self.har - self.hal) * 3600)
if win_fct is not None:
win_fct.to_meta(meta)
w_cube = datacube.CartWeightCube(weights, g_uu, g_vv, self.freqs, meta)
return SimuGridded(w_cube, self)
[docs]
class SimuGridded:
"""Gridded UV weights together with their telescope simulation."""
def __init__(self, weights, telescope_simu):
"""Store gridded weights and the simulation that produced them.
Args:
weights: A :class:`~ps_eor.datacube.CartWeightCube`, or the
metre-coordinate equivalent for a redundant array.
telescope_simu: Observation used to generate ``weights``.
"""
from ps_eor import psutil
self.weights = weights
self.telescope_simu = telescope_simu
self.name = self.telescope_simu.telescop.name
self.z = psutil.freq_to_z(self.weights.freqs.mean())
[docs]
def save(self, filename):
"""Save the weights and simulation metadata.
Args:
filename: Destination HDF5 filename.
"""
self.weights.meta.update(self.telescope_simu.to_dict())
self.weights.save(filename)
[docs]
@staticmethod
def load(filename):
"""Load weights and reconstruct their telescope simulation.
Args:
filename: HDF5 file created by :meth:`save`.
Returns:
SimuGridded: The restored coverage and observation settings.
"""
from ps_eor import datacube
weights = datacube.CartWeightCube.load(filename)
telescope_simu = TelescopeSimu.from_dict(weights.meta.kargs, weights.freqs)
if telescope_simu.telescop.redundant_array:
weights = datacube.CartWeightsCubeMeter(
weights.data, weights.uu, weights.vv, weights.freqs, weights.meta,
weights.uv_scale)
return SimuGridded(weights, telescope_simu)
[docs]
def get_slice(self, freq_start=None, freq_end=None):
"""Select an inclusive frequency interval.
Metre-coordinate redundant weights are converted to wavelength
coordinates at the interval's mean frequency.
Args:
freq_start, freq_end: Interval limits in Hz. Omitted limits use the
first or last available channel.
Returns:
SimuGridded: A new frequency-sliced simulation.
"""
from ps_eor import datacube
if freq_start is None:
freq_start = self.weights.freqs[0]
if freq_end is None:
freq_end = self.weights.freqs[-1]
weights = self.weights.get_slice(freq_start, freq_end)
if isinstance(weights, datacube.CartWeightsCubeMeter):
m_freq = (freq_end + freq_start) / 2.
weights = weights.get_cube(m_freq)
return SimuGridded(weights, self.telescope_simu)
[docs]
def get_sefd(self):
"""The single-polarization SEFD across the simulated frequencies.
Returns:
ndarray: SEFD in Jy, with at least one dimension.
"""
return np.atleast_1d(self.telescope_simu.telescop.get_sefd(self.weights.freqs))
[docs]
def get_i_sefd(self):
"""The Stokes-I SEFD across the simulated frequencies.
Returns:
ndarray: Stokes-I SEFD in Jy.
"""
return 1 / np.sqrt(2) * self.get_sefd()
[docs]
def get_ps_gen(self, filter_kpar_min=None, filter_wedge_theta=0):
"""Create a Cartesian power-spectrum estimator for this coverage.
The estimator inherits the telescope baseline limits and primary beam,
uses the simulated weights by default, and spans the full simulated
frequency interval. Pass the output of :meth:`get_noise_std_cube`
directly to its ``get_ps``, ``get_ps2d``, ``get_ps3d``, or ``get_all``
methods to estimate thermal-noise power. The returned spectrum's
``.err`` attribute is the corresponding one-sigma sensitivity.
Args:
filter_kpar_min: Optional minimum line-of-sight mode retained by the
estimator.
filter_wedge_theta: Optional wedge angle in degrees.
Returns:
ps_eor.pspec.PowerSpectraCart: Configured estimator.
"""
from ps_eor import datacube, pspec
du = 0.75 / self.weights.meta.theta_fov
if self.telescope_simu.telescop.redundant_array:
mfreq = self.weights.freqs.mean()
b = self.telescope_simu.telescop.redundant_baselines
el = 2 * np.pi * b / (const.c.value / mfreq)
else:
el = 2 * np.pi * (np.arange(self.weights.ru.min(), self.weights.ru.max(), du))
ps_conf = pspec.PowerSpectraConfig(el, window_fct='boxcar')
ps_conf.filter_kpar_min = filter_kpar_min
ps_conf.filter_wedge_theta = filter_wedge_theta
ps_conf.du = self.telescope_simu.telescop.du
ps_conf.umin = self.telescope_simu.umin
ps_conf.umax = self.telescope_simu.umax
ps_conf.weights_by_default = True
eor_bin_list = pspec.EorBinList(self.weights.freqs)
eor_bin_list.add_freq(1, self.weights.freqs.min() * 1e-6, self.weights.freqs.max() * 1e-6)
eor = eor_bin_list.get(1, self.weights.freqs)
pb = datacube.PrimaryBeam.from_name(self.telescope_simu.telescop.pb_name)
return pspec.PowerSpectraCart(eor, ps_conf, pb)
[docs]
def get_noise_std_cube(self, total_time_sec, sefd=None, min_weight=1):
"""Estimate the thermal-noise standard deviation in every UV cell.
A window recorded during gridding is included through its equivalent
noise-bandwidth correction. The returned cube can be passed directly
to the estimator from :meth:`get_ps_gen`; drawing an explicit noise
realization is not required for the expected noise spectrum. Its
``.err`` attribute gives the corresponding one-sigma sensitivity.
Args:
total_time_sec: Target total integration time in seconds.
sefd: Optional Stokes-I SEFD in Jy, scalar or one value per
frequency. By default use the telescope model.
min_weight: Remove cells whose weight is below this value.
Returns:
ps_eor.datacube.NoiseStdCube: Real standard deviation per frequency
and UV cell.
"""
if sefd is None:
sefd = self.get_i_sefd()
fake_apply_win_fct = 'PEWINFCT' in self.weights.meta
noise_std = self.weights.get_noise_std_cube(sefd, total_time_sec,
fake_apply_win_fct=fake_apply_win_fct)
noise_std.filter_min_weight(min_weight)
return noise_std