diff --git a/tangostationcontrol/tangostationcontrol/beam/test.py b/tangostationcontrol/tangostationcontrol/beam/test.py new file mode 100644 index 0000000000000000000000000000000000000000..59b6a87e5f9608e1abaf47a1e7dd23d0617d3296 --- /dev/null +++ b/tangostationcontrol/tangostationcontrol/beam/test.py @@ -0,0 +1,83 @@ +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 tempname: + def __init__(self, itrf: list([float])): + + """ Return a measure object, configured for the specified 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, 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 + +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 = tempname(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 + import pprint + pprint.pprint(delays)