SSB-47: redefinition of l, m before solving for the gains

CAL/CalibrationProcessing/lib/processing/solver.py

@@ -20,7 +20,7 @@ def _rotate_antenna_coordinates(dataset):
     station_position = numpy.array(dataset.target_station_position)
     antenna_field_position = numpy.array(dataset.antenna_field_position)
 
-    antenna_position_offset = numpy.array([antenna_position - station_position for
+    antenna_position_offset = numpy.array([station_position - antenna_position for
                                            antenna_position in antenna_field_position])
 
     station_rotation_matrix = dataset.rotation_matrix
@@ -42,7 +42,8 @@ def _compute_expected_phase_delay(l, m, x, y, frequency):
     return phase
 
 
-def _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, frequency):
+def _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, central_beam,
+                                                          frequency):
     """
     Compute the pointing matrix of a given station, at a given frequency and for a specific beam.
     :param dataset: datatable's dataset
@@ -51,6 +52,8 @@ def _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, fr
     :type datatable: dict(dict(numpy.ndarray))
     :param frequency: frequency at which to compute the pointing matrix
     :type frequency: float
+    :param central_beam: central beam
+    :type central_beam: str
     :return: the pointing matrix
     :rtype: numpy.matrix
     """
@@ -64,8 +67,8 @@ def _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, fr
 
     for i in range(n_beams):
         for j in range(n_antennas):
-            l = datatable[str(i)]['mean']['l'][0]
-            m = datatable[str(i)]['mean']['m'][0]
+            l = datatable[central_beam]['mean']['l'][0] - datatable[str(i)]['mean']['l'][0]
+            m = datatable[central_beam]['mean']['m'][0] - datatable[str(i)]['mean']['m'][0]
             x, y = rotated_coordinates[j, 0:2]
 
             phase = _compute_expected_phase_delay(l, m, x, y, frequency)
@@ -229,7 +232,7 @@ def _solve_gains(visibilities, matrix, **kwargs):
         return __empty
 
 
-def _solve_gains_per_frequency(dataset, datatable, frequency, direct_complex=False, **kwargs):
+def _solve_gains_per_frequency(dataset, datatable, frequency, direct_complex=True, **kwargs):
     """
     SOLVE THE EQUATION M * G = V FOR G
     :param dataset:
@@ -239,13 +242,27 @@ def _solve_gains_per_frequency(dataset, datatable, frequency, direct_complex=Fal
     :type frequency: float
     :return:
     """
-    matrix = _compute_pointing_matrices_per_station_frequency_beam(dataset, datatable, frequency)
+
+    central_beam = dataset.central_beamlets[str(frequency)]
+    matrix = _compute_pointing_matrices_per_station_frequency_beam(dataset,
+                                                                   datatable,
+                                                                   central_beam,
+                                                                   frequency)
     n_beams = len(datatable)
 
     result = dict()
+    #flags = numpy.array(
+    #    [datatable[str(i)]['mean']['flag'] for i in range(n_beams)])
+    #__empty = dict(gains=numpy.array(numpy.nan),
+    #               residual=numpy.array(numpy.nan),
+    #               relative_error=numpy.array(numpy.nan),
+    #               flag=numpy.array(True))
+    #is_frequency_flagged = flags
+
     for polarization in ['XX', 'XY', 'YX', 'YY']:
         visibilities = numpy.matrix(
             [datatable[str(i)]['mean'][polarization] for i in range(n_beams)])
+
         if direct_complex is True:
             result[polarization] = _solve_gains(visibilities, matrix, **kwargs)
         else: