diff --git a/CDB/LOFAR_ConfigDb.json b/CDB/LOFAR_ConfigDb.json
index 9c243aab268d9ca13f37bfe3d7af9cfbba8222de..de32af15daedc471ada0abd83ae9c1ec4e993346 100644
--- a/CDB/LOFAR_ConfigDb.json
+++ b/CDB/LOFAR_ConfigDb.json
@@ -90,6 +90,12 @@
"DigitalBeam": {
"STAT/DigitalBeam/1": {
"properties": {
+ "Beam_tracking_interval": [
+ "1.0"
+ ],
+ "Beam_tracking_preparation_period": [
+ "0.5"
+ ]
}
}
}
diff --git a/tangostationcontrol/tangostationcontrol/beam/delays.py b/tangostationcontrol/tangostationcontrol/beam/delays.py
index 850213c0cc64cdb72cb567a9b8c5334a464f4965..b1d7f9976ec5de3e75784dac484c9c33584f2806 100644
--- a/tangostationcontrol/tangostationcontrol/beam/delays.py
+++ b/tangostationcontrol/tangostationcontrol/beam/delays.py
@@ -1,5 +1,4 @@
import casacore.measures
-from math import sin, cos
import numpy
import datetime
@@ -29,26 +28,30 @@ class delay_calculator:
result = self.measure.do_frame(frame_time)
assert result == True
- def get_direction_vector(self, pointing: dict):
- """ Compute a direction vector for the given pointing, relative to the measure. """
- angles = self.measure.measure(pointing, "ITRF")
- angle0, angle1 = angles["m0"]["value"], angles["m1"]["value"]
+ def get_direction_vector(self, pointing: numpy.ndarray) -> numpy.ndarray:
+ """ Compute direction vector for a given pointing, relative to the measure. """
- # Convert polar to carthesian coordinates
- # see also https://github.com/casacore/casacore/blob/e793b3d5339d828a60339d16476bf688a19df3ec/casa/Quanta/MVDirection.cc#L67
- direction_vector = numpy.array(
- (cos(angle0) * cos(angle1),
- sin(angle0) * cos(angle1),
- sin(angle1)))
+ return self.get_direction_vector_bulk([pointing]).flatten()
+
+ def get_direction_vector_bulk(self, pointings: numpy.ndarray) -> numpy.ndarray:
+ """ Compute direction vectors for the given pointings, relative to the measure. """
- return direction_vector
+ angles0 = numpy.empty(len(pointings))
+ angles1 = numpy.empty(len(pointings))
+ for idx, pointing in enumerate(pointings):
+ angles = self.measure.measure(pointing, "ITRF")
+ angles0[idx] = angles["m0"]["value"]
+ angles1[idx] = angles["m1"]["value"]
- def _get_delay(self, reference_direction_vector: numpy.array, relative_itrf: numpy.array) -> float:
- """ Compute the delay to apply, in seconds, to align signals for a distance of `relative_itrf` in the
- direction of `reference_direction_vector`. """
- speed_of_light = 299792458.0
+ # Convert polar to carthesian coordinates
+ # see also https://github.com/casacore/casacore/blob/e793b3d5339d828a60339d16476bf688a19df3ec/casa/Quanta/MVDirection.cc#L67
+ direction_vectors = numpy.array(
+ [numpy.cos(angles0) * numpy.cos(angles1),
+ numpy.sin(angles0) * numpy.cos(angles1),
+ numpy.sin(angles1)])
- return numpy.dot(reference_direction_vector, relative_itrf) / speed_of_light
+ # Return array [directions][angles]
+ return direction_vectors.T
def is_valid_direction(self, direction):
try:
@@ -58,18 +61,34 @@ class delay_calculator:
return True
- def convert(self, direction, antenna_itrf: list([float])):
- # explicitly check validity to get a proper error message
- if not self.is_valid_direction(direction):
- raise ValueError(f"Invalid direction: {direction}")
+ def convert(self, direction, antenna_absolute_itrf: list([float])):
+ """ Get the delays for a direction and *absolute* antenna positions.
+
+ Returns delays[antenna]. """
- # obtain the direction vector for a specific pointing
- pointing = self.measure.direction(*direction)
+ return self.convert_bulk([direction], numpy.array(antenna_absolute_itrf) - self.reference_itrf).flatten()
- reference_dir_vector = self.get_direction_vector(pointing)
+ def convert_bulk(self, directions: numpy.ndarray, antenna_relative_itrfs: numpy.ndarray) -> numpy.ndarray:
+ """ Get the delays for each direction and each *relative* antenna position.
+
+ Returns delays[antenna][direction]. """
- # # compute the delays for an antennas w.r.t. the reference position
- antenna_relative_itrf = [subtract(pos, self.reference_itrf) for pos in antenna_itrf]
- delays = [self._get_delay(reference_dir_vector, relative_itrf) for relative_itrf in antenna_relative_itrf]
+ # obtain the direction vector for each pointing
+ try:
+ pointings = [self.measure.direction(*direction) for direction in directions]
+ except (RuntimeError, TypeError) as e:
+ raise ValueError("Invalid direction") from e
+
+ direction_vectors = self.get_direction_vector_bulk(pointings)
+
+ # compute the corresponding delays for all directions
+ def get_delay_all_directions(relative_itrf):
+ # Dot product between relative position and angle vector determines
+ # distance. Divide by speed of light to obtain delay.
+ return numpy.inner(relative_itrf, direction_vectors) / 299792458.0
+
+ # apply for each position
+ delays = numpy.apply_along_axis(get_delay_all_directions, 1, antenna_relative_itrfs)
return delays
+
diff --git a/tangostationcontrol/tangostationcontrol/clients/attribute_wrapper.py b/tangostationcontrol/tangostationcontrol/clients/attribute_wrapper.py
index 7f9c8d07fb3e5ceca540c2a03294dc9a92450d66..d3e6ee91b1177cb6e7306e07c448528cc476c8df 100644
--- a/tangostationcontrol/tangostationcontrol/clients/attribute_wrapper.py
+++ b/tangostationcontrol/tangostationcontrol/clients/attribute_wrapper.py
@@ -90,8 +90,11 @@ class attribute_wrapper(attribute):
"""
read_func_wrapper reads the attribute value, stores it and returns it"
"""
- device.value_dict[self] = self.read_function()
- return device.value_dict[self]
+ try:
+ device.value_dict[self] = self.read_function()
+ return device.value_dict[self]
+ except Exception as e:
+ raise e.__class__(f"Could not read attribute {comms_annotation}") from e
self.fget = read_func_wrapper
diff --git a/tangostationcontrol/tangostationcontrol/devices/sdp/digitalbeam.py b/tangostationcontrol/tangostationcontrol/devices/sdp/digitalbeam.py
index 9eb7533375f669acbcf652a5352deb47fdc5653e..2e985b8f10d66341a76bbe6b192dc3bf69677c8d 100644
--- a/tangostationcontrol/tangostationcontrol/devices/sdp/digitalbeam.py
+++ b/tangostationcontrol/tangostationcontrol/devices/sdp/digitalbeam.py
@@ -89,8 +89,8 @@ class DigitalBeam(beam_device):
# a delay calculator
self.delay_calculator = delay_calculator(reference_itrf)
- # absolute positions of each antenna element
- self.antenna_positions = antenna_itrf
+ # relatvie positions of each antenna
+ self.relative_antenna_positions = antenna_itrf - reference_itrf
# use all antennas unless specified otherwise
self._antenna_select = numpy.ones((self.NUM_ANTENNAS, self.NUM_BEAMLETS), dtype=numpy.bool_)
@@ -112,9 +112,8 @@ class DigitalBeam(beam_device):
d = self.delay_calculator
d.set_measure_time(timestamp)
- for beamlet in range(self.NUM_BEAMLETS):
- # calculate the delays based on the set reference position, the set time and now the set direction and antenna positions
- delays[:, beamlet] = d.convert(pointing_direction[beamlet], self.antenna_positions)
+ # calculate the delays based on the set reference position, the set time and now the set direction and antenna positions
+ delays = d.convert_bulk(pointing_direction, self.relative_antenna_positions)
return delays
diff --git a/tangostationcontrol/tangostationcontrol/test/beam/test_delays.py b/tangostationcontrol/tangostationcontrol/test/beam/test_delays.py
index b984660d677c81349c134d9be7c7441924253fc6..f2c537aafea96a8c28388bff40637bd8ca5eea73 100644
--- a/tangostationcontrol/tangostationcontrol/test/beam/test_delays.py
+++ b/tangostationcontrol/tangostationcontrol/test/beam/test_delays.py
@@ -1,12 +1,13 @@
import datetime
+import time
+import logging
+import numpy
+import numpy.testing
from tangostationcontrol.beam.delays import delay_calculator
from tangostationcontrol.test import base
-
-
-
class TestDelays(base.TestCase):
def test_init(self):
"""
@@ -94,7 +95,27 @@ class TestDelays(base.TestCase):
# calculate the delays based on the set reference position, the set time and now the set direction and antenna positions.
delays = d.convert(direction, antenna_itrf)
- self.assertListEqual(delays, [0.0], msg=f"delays = {delays}")
+ self.assertListEqual(delays.tolist(), [0.0], msg=f"delays = {delays}")
+
+ def test_regression(self):
+ reference_itrf = [3826577.066, 461022.948, 5064892.786] # CS002LBA, in ITRF2005 epoch 2012.5
+ d = delay_calculator(reference_itrf)
+
+ # set the antenna position identical to the reference position
+ antenna_itrf = [[3826923.503, 460915.488, 5064643.517]] # CS001LBA, in ITRF2005 epoch 2012.5
+
+ # # set the timestamp to solve for
+ timestamp = datetime.datetime(2000, 1, 1, 0, 0, 0)
+ d.set_measure_time(timestamp)
+
+ # # obtain the direction vector for a specific pointing
+ direction = "J2000", "0deg", "90deg"
+
+ # calculate the delays based on the set reference position, the set time and now the set direction and antenna positions.
+ delays = d.convert(direction, antenna_itrf)
+
+ # check for regression
+ self.assertAlmostEqual(-8.31467564781444e-07, delays[0], delta=1.0e-10)
def test_light_second_delay(self):
"""
@@ -123,3 +144,39 @@ class TestDelays(base.TestCase):
self.assertAlmostEqual(0.1, delays[0], 6, f"delays[0] = {delays[0]}")
+
+ def test_convert_bulk(self):
+ d = delay_calculator([0, 0, 0])
+ timestamp = datetime.datetime(2022, 3, 1, 0, 0, 0) # timestamp does not actually matter, but casacore doesn't know that.
+ d.set_measure_time(timestamp)
+
+ # generate different positions and directions
+ positions = numpy.array([[i,2,3] for i in range(5)])
+ directions = numpy.array([["J2000", f"{i}deg", f"{i}deg"] for i in range(90)])
+
+ bulk_result = d.convert_bulk(directions, positions)
+
+ # verify parallellisation along direction axis
+ for i, single_dir in enumerate(directions):
+ single_dir_result = d.convert_bulk([single_dir], positions)
+ numpy.testing.assert_equal(single_dir_result[:, 0], bulk_result[:, i])
+
+ # verify parallellisation along position axis
+ for i, single_pos in enumerate(positions):
+ single_pos_result = d.convert_bulk(directions, [single_pos])
+ numpy.testing.assert_equal(single_pos_result[0, :], bulk_result[i, :])
+
+ def test_convert_bulk_speed(self):
+ d = delay_calculator([0, 0, 0])
+ timestamp = datetime.datetime(2022, 3, 1, 0, 0, 0) # timestamp does not actually matter, but casacore doesn't know that.
+ d.set_measure_time(timestamp)
+
+ positions = numpy.array([[1,2,3]] * 96)
+ directions = numpy.array([["J2000", "0deg", "0deg"]] * 488)
+
+ count = 10
+ before = time.monotonic_ns()
+ for _ in range(count):
+ _ = d.convert_bulk(directions, positions)
+ after = time.monotonic_ns()
+ logging.error(f"convert_bulk averaged {(after - before)/count/1e6} ms to convert 488 directions for 96 antennas.")