diff --git a/tangostationcontrol/tangostationcontrol/beam/delays.py b/tangostationcontrol/tangostationcontrol/beam/delays.py new file mode 100644 index 0000000000000000000000000000000000000000..771fad6c90a59ebe002867fedfc55446ead51000 --- /dev/null +++ b/tangostationcontrol/tangostationcontrol/beam/delays.py @@ -0,0 +1,62 @@ +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 diff --git a/tangostationcontrol/tangostationcontrol/beam/test_delays.py b/tangostationcontrol/tangostationcontrol/beam/test_delays.py new file mode 100644 index 0000000000000000000000000000000000000000..32eee0838224abfe16ad47b11f3bd0c9af0f9e93 --- /dev/null +++ b/tangostationcontrol/tangostationcontrol/beam/test_delays.py @@ -0,0 +1,66 @@ +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}")