From 258d52c3c2a01fa64a3f674129f78b8baeb87ad4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Thomas=20J=C3=BCrges?= <jurges@astron.nl>
Date: Fri, 5 Oct 2018 15:55:45 +0000
Subject: [PATCH] Task SSB-43: Added code for the central beamlets.
---
.../lib/datacontainers/holography_dataset.py | 83 ++++++++++++++++++-
.../holography_dataset_definitions.py | 1 +
2 files changed, 82 insertions(+), 2 deletions(-)
diff --git a/CAL/CalibrationCommon/lib/datacontainers/holography_dataset.py b/CAL/CalibrationCommon/lib/datacontainers/holography_dataset.py
index d1a53b791b3..ac3133d845a 100644
--- a/CAL/CalibrationCommon/lib/datacontainers/holography_dataset.py
+++ b/CAL/CalibrationCommon/lib/datacontainers/holography_dataset.py
@@ -1,5 +1,5 @@
import logging
-
+import math
import h5py
from lofar.calibration.common.datacontainers.holography_observation import HolographyObservation
@@ -47,6 +47,10 @@ class HolographyDataset():
# list of beamlet numbers
self.beamlets = list()
+
+ # The central beam or beamlet for each frequency
+ self.central_beamlets = dict()
+
# coordinates of the antenna position in the target
self.antenna_field_position = []
# list of reference station names
@@ -149,6 +153,75 @@ class HolographyDataset():
return False
return equality
+ def find_central_beamlets(self, source, ra_dec, frequencies, beamlets):
+ '''
+ This function finds the central beamlet of a target station for every
+ frequency.
+ @param source: A dict containing right ascension (source["RA"]) and
+ declination (source["DEC"]) of the source.
+ @param ra_dec: A dict[frequency][beamlet] that contains the right
+ ascension (ra_dec[frequency][beamlet][0]) and declination
+ (ra_dec[frequency][beamlet][1]) for every beamlet.
+ @param frequencies: A dict that contains the frequencies at which the
+ holography observation was performed.
+ @param beamlets: A list of beamlet numbers.
+ '''
+ logger.debug("Finding the central beamlets...")
+ source_ra = source["RA"]
+ source_dec = source["DEC"]
+ # Store each pseudo distance in a dict of dicts: pd[frequency][beamlet]
+ pseudo_distance = dict()
+ # Iterate over all frequencies...
+ for frequency in frequencies:
+ frequency_string = str(frequency)
+ pseudo_distance[frequency_string] = dict()
+
+ # and iterate over the beamlets.
+ for beamlet in beamlets:
+ beamlet_string = str(beamlet)
+ (ra, dec, epoch) = ra_dec[frequency_string][beamlet_string]
+ # calculate a pseudo distance that is an indicator if this
+ # beamlet is closer to the source than the ones before.
+ # Note that I am too lazy to calculate the spherical distance.
+ # But keep in mind that this pseudo distance like its cousin
+ # the geometrical distance are not really cutting it since RA
+ # and DEC are coordinates of a spherical coordinate system.
+ # But here the pseudo distance is good enough.
+ pseudo_distance[frequency][beamlet] = math.abs(ra - source_ra) + math.abs(dec - source_dec)
+
+ # OK. Done with all the iteration business. Now check if across all
+ # frequencies the same beamlet is the central one. It is allowed that
+ # the central beamlets are not always the same but it would be better
+ # and let the astronomers sleep better if it would always be the same.
+ # And also check if there is actually only one central beamlet per
+ # frequency.
+ central_beamlet = dict()
+ for (frequency_string, beamlets) in pseudo_distance.items():
+ values = beamlets.values()
+ keys = beamlets.keys()
+ minimal_distance = min(values)
+ # If I find a second beamlet that is 'closest' to the source for the
+ # same frequency then that is a problem. The outset was that there
+ # is only a single central beamlet and not a bunch of 'em.
+ if values.count(minimal_distance) == 1:
+ # All is good. Only one central beamlet.
+ central_beamlet[frequency_string] = keys[values.index(minimal_distance)]
+ else:
+ text = "Found %d beamlets that have the same distance from the source position. Therefore a unique central beamlet does not exist." % (values.count(minimal_distance))
+ logger.error(text)
+ raise ValueError(text)
+
+ # Now check if the central beamlet is the same across all frequencies.
+ # I do that by creating a set from the central beamlets of all stations.
+ # If there is only one central beamlet for all frequencies, the size of
+ # the set will be 1.
+ central_beamlet_set = set(central_beamlet.values())
+ if len(central_beamlet_set) == 1:
+ logger.debug("All is good, unicorns everywhere, there is only one central beamlet \"%s\" for all frequencies.", self.central_beamlet)
+ else:
+ logger.warn("Multiple central beamlets have been identified: ", central_beamlet)
+ return central_beamlet
+
def load_from_beam_specification_and_ms(self, station_name, list_of_hbs_ms_tuples):
"""
Loads the dataset from the specification files and the measurements for the given station
@@ -160,6 +233,7 @@ class HolographyDataset():
logger.info("Creating a holography data set for station \"%s\"...", station_name)
self.__collect_preliminary_information(station_name, list_of_hbs_ms_tuples)
self.__read_data(station_name, list_of_hbs_ms_tuples)
+ self.central_beamlets = self.find_central_beamlets(self.source_position, self.ra_dec, self.frequencies, self.beamlets)
logger.info("Creation of a holography data set for station \"%s\" done.", station_name)
def __read_data(self, station_name, list_of_hbs_ms_tuples):
@@ -216,6 +290,7 @@ class HolographyDataset():
target_stations = set()
reference_stations = set()
beamlets = set()
+ central_beamlet = None
virtual_pointing = dict()
frequencies = set()
sas_ids = set()
@@ -339,6 +414,7 @@ class HolographyDataset():
logger.info("Target station position = %s", hds.target_station_position)
logger.info("Source name = %s", hds.source_name)
logger.info("Source position = %s", hds.source_position)
+ logger.info("Central beamlet = %s", hds.central_beamlet)
logger.info("Start time = %s", hds.start_time)
logger.info("End time = %s", hds.end_time)
logger.info("Rotation matrix = %s", hds.rotation_matrix)
@@ -444,12 +520,13 @@ class HolographyDataset():
result.target_station_name = f.attrs[HDS_TARGET_STATION_NAME]
result.target_station_position = list(f.attrs[HDS_TARGET_STATION_POSITION])
result.source_name = f.attrs[HDS_SOURCE_NAME]
- result.source_position = tuple(f.attrs[HDS_SOURCE_POSITION].tolist())
+ result.source_position = numpy.array(f.attrs[HDS_SOURCE_POSITION])
(result.start_time, result.end_time) = f.attrs[HDS_OBSERVATION_TIME]
result.rotation_matrix = f.attrs[HDS_ROTATION_MATRIX]
result.antenna_field_position = f.attrs[HDS_ANTENNA_FIELD_POSITION].tolist()
result.reference_stations = list(f[HDS_REFERENCE_STATION])
result.frequencies = list(f[HDS_FREQUENCY])
+ result.central_beamlets = f.attrs[HDS_CENTRAL_BEAMLETS]
result.ra_dec = dict()
for frequency in f["RA_DEC"].keys():
@@ -521,6 +598,8 @@ class HolographyDataset():
f.attrs[HDS_ROTATION_MATRIX] = self.rotation_matrix
f.attrs[HDS_ANTENNA_FIELD_POSITION] = self.antenna_field_position
+ f.attrs[HDS_CENTRAL_BEAMLETS] = self.central_beamlets
+
# Store the list of reference stations and frequencies. We just
# want to keep 'em around for quick reference.
f[HDS_REFERENCE_STATION] = self.reference_stations
diff --git a/CAL/CalibrationCommon/lib/datacontainers/holography_dataset_definitions.py b/CAL/CalibrationCommon/lib/datacontainers/holography_dataset_definitions.py
index 837300bdd62..56123974a28 100644
--- a/CAL/CalibrationCommon/lib/datacontainers/holography_dataset_definitions.py
+++ b/CAL/CalibrationCommon/lib/datacontainers/holography_dataset_definitions.py
@@ -18,6 +18,7 @@ HDS_ROTATION_MATRIX = "Rotation_matrix"
HDS_ANTENNA_FIELD_POSITION = "Antenna_field_position"
HDS_SPECIFIED_REFERENCE_STATION = "Specified_Reference_station"
HDS_SPECIFIED_FREQUENCY = "Specified_frequency"
+HDS_CENTRAL_BEAMLETS = "Central_beamlets"
# GROUP "RA DEC"
HDS_SPECIFIED_RA_DEC = "Specified_RA_DEC"
--
GitLab