LOFAR telescope: correcting for the array factor in a given direction

The EveryBeam Python bindings can be used to evaluate the station response in an arbitrary equatorial direction. This demo illustrates in a step-by-step manner how to do this. The complete code is shown at the bottom of this page.

It all starts with importing the everybeam python module, along with any other libraries that are needed for your project, and by specifying a path to the Measurement Set that will be used:

from astropy.coordinates import AltAz, EarthLocation, ITRS, SkyCoord
from astropy.time import Time

import astropy.units as u
import everybeam as eb
import numpy as np

# Set path to LOFAR LBA MS and load telescope
ms_path = "/path/to/my.ms"

The telescope can now be loaded with eb.load_telescope. This should return an instance of a LOFAR telescope.

telescope = eb.load_telescope(ms_path)
assert type(telescope) == eb.LOFAR

Please note that since no optional arguments were passed to load_telescope, no differential beam will applied and the default element response model ("Hamaker") will be used.

In order to evaluate the array factor, we need to pass some user input. This includes a time, station index, frequency, the direction for which we want to know the array factor, and the reference direction (where the array factor is unity). For the array_factor method these coordinates need to be in ITRF coordinates. We therefore also need to convert from equatorial IRCS (J2000) coordinates to ITRF coordinates.

First, let’s read observation time, frequency and pointing direction information from the Measurement Set.

import casacore.tables as pt
# Time slots at which to evaluate the beam response.
ms_times = pt.taql('SELECT DISTINCT TIME FROM {ms:s}'.format(ms=ms_path))
times = ms_times.getcol('TIME')

# Frequencies at which to evaluate the beam response.
ms_freqs = pt.taql('SELECT CHAN_FREQ FROM {ms:s}::SPECTRAL_WINDOW'.format(ms=ms_path))
freqs = ms_freqs.getcol('CHAN_FREQ').squeeze()

# Obtain the reference direction from the Measurement Set.
ms_dirs = pt.taql('SELECT REFERENCE_DIR FROM {ms:s}::FIELD'.format(ms=ms_path))
ra_ref, dec_ref = ms_dirs.getcol('REFERENCE_DIR').squeeze()

Next, define a few (phase) directions for which we would like to calculate the array factor. We need two separate arrays for RA and Dec. The directions should also be in ICRS (J2000) in units of radians.

# Choose three random directions (units in radians), and create separate RA and DEC arrays.
ra, dec = list(zip((1.34, 1.56), (0.78, -0.14), (-0.57, 0.38)))

Now that we have our initial information we can calculate the ITRF coordinates. We use AstroPy’s coordinate transformations to obtain the coordinates as cartesian coordinates in ITRS. Let’s define a function for this, because we will need it more often.

def radec_to_xyz(ra, dec, time):
    """
    Convert RA and Dec ICRS coordinates to ITRS cartesian coordinates.

    Args:
        ra (astropy.coordinates.Angle): Right ascension
        dec (astropy.coordinates.Angle): Declination
        time (float): MJD time in seconds

    Returns:
        pointing_xyz (ndarray): NumPy array containing the ITRS X, Y and Z coordinates
    """
    obstime = Time(time/3600/24, scale='utc', format='mjd')
    dir_pointing = SkyCoord(ra, dec)
    dir_pointing_itrs = dir_pointing.transform_to(ITRS(obstime=obstime))
    return np.asarray(dir_pointing_itrs.cartesian.xyz.transpose())

Note

Our radec_to_xyz function can take arrays of RA and Dec, and an array of times, but not simultaneously. So, either calculate the ITRF coordinates for a single direction for an array of MJD times, or calculate the ITRF coordinates for an array of directions for a single MJD time.

Next, convert the reference direction and our (phase) directions of interest to ITRF. Note that the coordinates in ITRF are time-dependent, due to the rotation of the earth. So let’s use the first time slot, and select five stations.

# ITRF coordinates of the reference direction.
reference_xyz = radec_to_xyz(ra_ref * u.rad, dec_ref * u.rad, times[0])
# ITRF coordinates of the phase centre to correct the array factor for.
phase_xyz = radec_to_xyz(ra * u.rad, dec * u.rad, times[0])
# Station IDs for which we want to calculate the array factor.
station_ids = [0, 3, 6, 14, 23]

We’re all set now to compute the array factor response using the array_factor method. Let’s calculate it for the given station ID and the first time slot:

array_factor = telescope.array_factor(times[0], station_ids, freqs, phase_xyz, reference_xyz)

This returns a 5-dimensional array. The outermost (first) dimension is the station ID (len(station_ids)), the second is the frequency (len(freqs)), the third is the array of phase directions (len(ra) or len(dec)), and the two innermost (fourth and fifth) are a 2x2 array with the response of the XX, XY, YX and YY.

Note

array_factor will always return a squeezed array, if you use scalar values for station_ids, freqs, or phase_xyz. Hence, it will return a 2-dimensional array if you calculate the array factor for a single station, for a single frequency, for a single phase direction.

Since EveryBeam’s methods do not take arrays of time, you will need to calculate the array factor for the full time range using a loop:

# Create an empty 6-dim complex numpy array to hold the array factor for each time slot
timeslices = np.empty((len(times), len(station_ids), len(freqs), len(ra), 2, 2), dtype=np.complex128)
for idx, time in enumerate(times):
    reference_xyz = radec_to_xyz(ra_ref * u.rad, dec_ref * u.rad, time)
    phase_xyz = radec_to_xyz(ra * u.rad, dec * u.rad, time)
    beam = telescope.array_factor(time, station_ids, freqs, phase_xyz, reference_xyz)
    timeslices[idx] = beam

Caution

The code above merely serves an educational purpose, demonstrating how one can collect all the desired information in one big NumPy array. In practice, you will neither want to calculate the array factor for every time slot in the Measurement Set, nor for every frequency channel. The array factor will change slowly over time and frequency. Hence, calculating it once every 15 minutes or so, for just a couple of frequencies, usually suffices in practice.

A complete overview of the code is shown below:

#!/usr/bin/env python3
# Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
# SPDX-License-Identifier: GPL-3.0-or-later

import os

import astropy.units as u
import casacore.tables as pt
import everybeam as eb
import numpy as np
from astropy.coordinates import ITRS, SkyCoord
from astropy.time import Time


def radec_to_xyz(ra, dec, time):
    """
    Convert RA and Dec ICRS coordinates to ITRS cartesian coordinates.

    Args:
        ra (astropy.coordinates.Angle): Right ascension
        dec (astropy.coordinates.Angle): Declination
        time (float): MJD time in seconds

    Returns:
        pointing_xyz (ndarray): NumPy array containing the ITRS X, Y and Z coordinates
    """
    obstime = Time(time / 3600 / 24, scale="utc", format="mjd")
    dir_pointing = SkyCoord(ra, dec)
    dir_pointing_itrs = dir_pointing.transform_to(ITRS(obstime=obstime))
    return np.asarray(dir_pointing_itrs.cartesian.xyz.transpose())


# Set path to LOFAR LBA MS and load telescope
ms_path = os.path.join(os.environ["DATA_DIR"], "LOFAR_LBA_MOCK.ms")

telescope = eb.load_telescope(ms_path)
assert type(telescope) == eb.LOFAR

# Time slots at which to evaluate the beam response.
ms_times = pt.taql("SELECT DISTINCT TIME FROM {ms:s}".format(ms=ms_path))
times = ms_times.getcol("TIME")

# Frequencies at which to evaluate the beam response.
ms_freqs = pt.taql(
    "SELECT CHAN_FREQ FROM {ms:s}::SPECTRAL_WINDOW".format(ms=ms_path)
)
freqs = ms_freqs.getcol("CHAN_FREQ").squeeze()

# Obtain the reference direction from the Measurement Set.
ms_dirs = pt.taql("SELECT REFERENCE_DIR FROM {ms:s}::FIELD".format(ms=ms_path))
ra_ref, dec_ref = ms_dirs.getcol("REFERENCE_DIR").squeeze()

# Choose three random directions (units in radians), and create separate RA and DEC arrays.
ra, dec = list(zip((1.34, 1.56), (0.78, -0.14), (-0.57, 0.38)))

# ITRF coordinates of the reference direction.
reference_xyz = radec_to_xyz(ra_ref * u.rad, dec_ref * u.rad, times[0])

# ITRF coordinates of the phase centre to correct the array factor for.
phase_xyz = radec_to_xyz(ra * u.rad, dec * u.rad, times[0])

# Station IDs for which we want to calculate the array factor.
station_ids = [0, 3, 6, 14, 23]

# Compute the array factor response
array_factor = telescope.array_factor(
    times[0], station_ids, freqs, phase_xyz, reference_xyz
)
print(
    f"\n******** array_factor {array_factor.shape} ********\n\n{array_factor}\n"
)

# Create an empty 6-dim complex numpy array to hold the array factor for each time slot
timeslices = np.empty(
    (len(times), len(station_ids), len(freqs), len(ra), 2, 2),
    dtype=np.complex128,
)
for idx, time in enumerate(times):
    reference_xyz = radec_to_xyz(ra_ref * u.rad, dec_ref * u.rad, time)
    phase_xyz = radec_to_xyz(ra * u.rad, dec * u.rad, time)
    beam = telescope.array_factor(
        time, station_ids, freqs, phase_xyz, reference_xyz
    )
    timeslices[idx] = beam
print(f"\n******** timeslices {timeslices.shape} ********\n\n{timeslices}\n")

print(
    f"timeslices[0].all() == array_factor.all(): {timeslices[0].all() == array_factor.all()}"
)