Task SSB-43: Added code for the central beamlets.

CAL/CalibrationCommon/lib/datacontainers/holography_dataset.py

 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

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"