Merge branch 'L2SS-479_09_dec_2021-convert-pointing-to-delays' into 'master'

Resolve L2SS-479 "09 dec 2021 convert pointing to delays"

Closes L2SS-479

See merge request !194

tangostationcontrol/tangostationcontrol/beam/delays.py

+import casacore.measures
+from math import sin, cos
+import numpy
+import datetime
+
+def subtract(a, b):
+    return numpy.array([x - y for x, y in zip(a, b)])
+
+
+class delay_calculator:
+    def __init__(self, itrf: list([float])):
+
+        """ Create a measure object, configured for the specified terrestrial location. """
+
+        measure = casacore.measures.measures()
+        frame_location = measure.position("ITRF", *[f"{x}m" for x in itrf])
+
+        result = measure.do_frame(frame_location)
+        assert result == True
+
+        self.reference_itrf = itrf
+        self.measure = measure
+        self.measure_time = None
+
+    def set_measure_time(self, utc_time: datetime.datetime):
+        """ Configure the measure object for the specified time. """
+        frame_time = self.measure.epoch("UTC", utc_time.isoformat(' '))
+
+        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"]
+
+        # 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 direction_vector
+
+    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
+
+        return numpy.dot(reference_direction_vector, relative_itrf) / speed_of_light
+
+    def convert(self, direction, antenna_itrf: list([float])):
+        # obtain the direction vector for a specific pointing
+        pointing = self.measure.direction(*direction)
+        reference_dir_vector = self.get_direction_vector(pointing)
+
+        # # 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]
+
+        return delays

tangostationcontrol/tangostationcontrol/beam/test_delays.py

+from delays import *
+import pprint
+
+if __name__ == '__main__':
+    # # create a frame tied to the reference position
+    reference_itrf = [3826577.066, 461022.948, 5064892.786] # CS002LBA, in ITRF2005 epoch 2012.5
+    d = delay_calculator(reference_itrf)
+
+    # # set the timestamp to solve for
+    timestamp = datetime.datetime(2021,1,1,0,0,5)
+    d.set_measure_time(timestamp)
+
+    # compute the delays for an antennas w.r.t. the reference position
+    antenna_itrf = [[3826923.546, 460915.441, 5064643.489]] # CS001LBA, in ITRF2005 epoch 2012.5
+
+    # # obtain the direction vector for a specific pointing
+    direction = "J2000","0deg","0deg"
+
+    # 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)
+
+    # print the delays
+    # pprint.pprint(delays)
+
+
+    #test changing the time
+
+    print(f"Changing timestamp test\nBase parametres: Direction: {direction}, position: {antenna_itrf}")
+    for i in range(10):
+        # # set the timestamp to solve for
+        timestamp = datetime.datetime(2021,1,1,0,i,5)
+        d.set_measure_time(timestamp)
+
+        delays = d.convert(direction, antenna_itrf)
+
+        # print the delays
+        print(f"Timestamp: {timestamp}:   {delays}")
+
+
+    # reset time
+    timestamp = datetime.datetime(2021, 1, 1, 0, 0, 5)
+    d.set_measure_time(timestamp)
+
+
+    #test changing the antenna position
+    print(f"Changing Antenna position test.\nBase parametres: Time: {timestamp} Direction: {direction}")
+    for i in range(10):
+        antenna_itrf = [[3826577.066 + i, 461022.948, 5064892.786]]  # CS002LBA, in ITRF2005 epoch 2012.5
+
+        delays = d.convert(direction, antenna_itrf)
+
+        # print the delays
+        print(f"Antenna position: {antenna_itrf}:   {delays}")
+
+    # test changing the direction
+
+    antenna_itrf = [[3826923.546, 460915.441, 5064643.489]]  # CS001LBA, in ITRF2005 epoch 2012.5
+    print(f"Changing direction test.\nBase parametres: Time: {timestamp} , position: {antenna_itrf}")
+
+    for i in range(10):
+        direction = "J2000", f"{i}deg", "0deg"
+
+        delays = d.convert(direction, antenna_itrf)
+
+        # print the delays
+        print(f"Direction: {direction}:  {delays}")