Fix MWA Efield beam sign

CHANGELOG.md

@@ -9,6 +9,10 @@ the data array are single precision. Default is set to write them as doubles.
 - ATA has been added to the list of known telescopes.
 
 ### Fixed
+- A sign flip of the MWA beam response to azimuthally aligned polarization. This
+was caused by the difference in coordinates used in UVBeam and the mwa_pb repo,
+the coordinates themselves were properly accounted for but the flip in
+orientation of the basis vectors was not properly handled.
 - Bug in the `uvh5._add` method that expected keys which could be
 optional for some phase-center catalog-entry types. This was noticed as an issue when
 adding two uvh5 objects that were of type "sidereal" which did not have "info_source"

docs/analytic_beam_tutorial.rst

@@ -69,21 +69,28 @@ are included, the array returned from the ``power_eval`` method will be complex.
   ...   beam_vals[0,0,0],
   ...   norm=LogNorm(vmin = 1e-8, vmax =1),
   ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
   ... )
   >>> _ = ax[0].set_title(f"Airy beam {freqs[0]*1e-6} MHz")
-  >>> _ = ax[0].set_xlabel("direction cosine l")
-  >>> _ = ax[0].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(bp_low, ax=ax[0], fraction=0.046, pad=0.04)
+  >>> _ = fig.colorbar(bp_low, ax=ax[0], fraction=0.046, pad=0.04, location="left")
 
   >>> bp_high = ax[1].imshow(
   ...   beam_vals[0,0,-1],
   ...   norm=LogNorm(vmin = 1e-8, vmax =1),
   ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
   ... )
   >>> _ = ax[1].set_title(f"Airy beam {freqs[-1]*1e-6} MHz")
-  >>> _ = ax[1].set_xlabel("direction cosine l")
-  >>> _ = ax[1].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(bp_high, ax=ax[1], fraction=0.046, pad=0.04)
+  >>> _ = fig.colorbar(bp_high, ax=ax[1], fraction=0.046, pad=0.04, location="left")
+
+  >>> for ind in range(2):
+  ...   _ = ax[ind].set_xticks([0], labels=["North"])
+  ...   _ = ax[ind].set_yticks([0], labels=["East"])
+  ...   _ = ax[ind].yaxis.set_label_position("right")
+  ...   _ = ax[ind].yaxis.tick_right()
+  ...   _ = ax[ind].xaxis.set_label_position("top")
+  ...   _ = ax[ind].xaxis.tick_top()
+
   >>> fig.tight_layout()
   >>> plt.show()  # doctest: +SKIP
   >>> plt.savefig("Images/airy_beam.png", bbox_inches='tight')
@@ -144,29 +151,46 @@ part, but in general E-Field beams can be complex, so a complex array is returne
 
   >>> fig, ax = plt.subplots(2, 2)
 
-  >>> be00 = ax[0,0].imshow(beam_vals[0,0,0].real, extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)])
+  >>> be00 = ax[0,0].imshow(
+  ...   beam_vals[0,0,0].real,
+  ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
+  ... )
   >>> _ = ax[0,0].set_title("E/W dipole azimuth response")
-  >>> _ = ax[0,0].set_xlabel("direction cosine l")
-  >>> _ = ax[0,0].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(be00, ax=ax[0,0])
+  >>> _ = fig.colorbar(be00, ax=ax[0,0], location="left")
 
-  >>> be10 = ax[1,0].imshow(beam_vals[1,0,0].real, extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)])
+  >>> be10 = ax[1,0].imshow(
+  ...   beam_vals[1,0,0].real,
+  ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
+  ... )
   >>> _ = ax[1,0].set_title("E/W dipole zenith angle response")
-  >>> _ = ax[1,0].set_xlabel("direction cosine l")
-  >>> _ = ax[1,0].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(be00, ax=ax[1,0])
+  >>> _ = fig.colorbar(be10, ax=ax[1,0], location="left")
 
-  >>> be01 = ax[0,1].imshow(beam_vals[0,1,0].real, extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)])
+  >>> be01 = ax[0,1].imshow(
+  ...   beam_vals[0,1,0].real,
+  ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
+  ... )
   >>> _ = ax[0,1].set_title("N/S dipole azimuth response")
-  >>> _ = ax[0,1].set_xlabel("direction cosine l")
-  >>> _ = ax[0,1].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(be00, ax=ax[0,1])
+  >>> _ = fig.colorbar(be01, ax=ax[0,1], location="left")
 
-  >>> be11 = ax[1,1].imshow(beam_vals[1,1,0].real, extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)])
+  >>> be11 = ax[1,1].imshow(
+  ...   beam_vals[1,1,0].real,
+  ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
+  ... )
   >>> _ = ax[1,1].set_title("N/S dipole zenith angle response")
-  >>> _ = ax[1,1].set_xlabel("direction cosine l")
-  >>> _ = ax[1,1].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(be00, ax=ax[1,1])
+  >>> _ = fig.colorbar(be11, ax=ax[1,1], location="left")
+
+  >>> for row_ind in range(2):
+  ...   for col_ind in range(2):
+  ...      _ = ax[row_ind,col_ind].set_xticks([0], labels=["North"])
+  ...      _ = ax[row_ind,col_ind].set_yticks([0], labels=["East"])
+  ...      _ = ax[row_ind,col_ind].yaxis.set_label_position("right")
+  ...      _ = ax[row_ind,col_ind].yaxis.tick_right()
+  ...      _ = ax[row_ind,col_ind].xaxis.set_label_position("top")
+  ...      _ = ax[row_ind,col_ind].xaxis.tick_top()
 
   >>> fig.tight_layout()
   >>> plt.show()  # doctest: +SKIP

docs/beam_interface_tutorial.rst

@@ -79,21 +79,28 @@ that is attached to it is an analytic beam or a UVBeam.
   ...   dipole_beam_vals[0].real,
   ...   norm=LogNorm(vmin = 1e-4, vmax =1),
   ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
   ... )
   >>> _ = ax[0].set_title(f"E/W Dipole power beam")
-  >>> _ = ax[0].set_xlabel("direction cosine l")
-  >>> _ = ax[0].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(bp_dip, ax=ax[0], fraction=0.046, pad=0.04)
+  >>> _ = fig.colorbar(bp_dip, ax=ax[0], fraction=0.046, pad=0.04, location="left")
 
   >>> bp_mwa = ax[1].imshow(
   ...   mwa_beam_vals[0].real,
   ...   norm=LogNorm(vmin = 1e-4, vmax =1),
   ...   extent=[np.min(l_arr), np.max(l_arr), np.min(m_arr), np.max(m_arr)],
+  ...   origin="lower",
   ... )
   >>> _ = ax[1].set_title(f"MWA E/W power beam")
-  >>> _ = ax[1].set_xlabel("direction cosine l")
-  >>> _ = ax[1].set_ylabel("direction cosine m")
-  >>> _ = fig.colorbar(bp_mwa, ax=ax[1], fraction=0.046, pad=0.04)
+  >>> _ = fig.colorbar(bp_mwa, ax=ax[1], fraction=0.046, pad=0.04, location="left")
+
+  >>> for ind in range(2):
+  ...   _ = ax[ind].set_xticks([0], labels=["North"])
+  ...   _ = ax[ind].set_yticks([0], labels=["East"])
+  ...   _ = ax[ind].yaxis.set_label_position("right")
+  ...   _ = ax[ind].yaxis.tick_right()
+  ...   _ = ax[ind].xaxis.set_label_position("top")
+  ...   _ = ax[ind].xaxis.tick_top()
+
   >>> fig.tight_layout()
   >>> plt.show()  # doctest: +SKIP
   >>> plt.savefig("Images/dipole_mwa_power.png", bbox_inches='tight')

docs/uvdata_tutorial.rst

@@ -919,7 +919,7 @@ Note: there is now support for reading in only part of a uvfits, uvh5 or miriad
 
   >>> # Amplitude waterfall for all spectral channels and 0th polarization
   >>> fig, ax = plt.subplots(1, 1)
-  >>> _ = ax.imshow(np.abs(waterfall_data), interpolation='none')
+  >>> _ = ax.imshow(np.abs(waterfall_data), interpolation='none', origin="lower")
   >>> _ = ax.set_yticks([0, waterfall_times.size - 1])
   >>> _ = ax.set_yticklabels([waterfall_times[0], waterfall_times[1]])
   >>> freq_tick_inds = np.concatenate((np.arange(0, uvd.Nfreqs, 16), [uvd.Nfreqs-1]))

src/pyuvdata/beam_interface.py

@@ -266,13 +266,15 @@ def compute_response(
         spline_opts : dict
             Provide options to numpy.RectBivariateSpline. This includes spline
             order parameters `kx` and `ky`, and smoothing parameter `s`. Only
-            applies if beam is a UVBeam and interpolation_function is "az_za_simple".
+            applies if beam is a UVBeam and interpolation_function is "az_za_simple"
+            or "az_za_map_coordinates".
         check_azza_domain : bool
             Whether to check the domain of az/za to ensure that they are covered by the
             intrinsic data array. Checking them can be quite computationally expensive.
             Conversely, if the passed az/za are outside of the domain, they will be
             silently extrapolated and the behavior is not well-defined. Only
-            applies if beam is a UVBeam and interpolation_function is "az_za_simple".
+            applies if beam is a UVBeam and interpolation_function is "az_za_simple"
+            or "az_za_map_coordinates".
 
         Returns
         -------

src/pyuvdata/uvbeam/mwa_beam.py

@@ -199,14 +199,17 @@ class MWABeam(UVBeam):
     This class should not be interacted with directly, instead use the
     read_mwa_beam method on the UVBeam class.
 
-    This is based on https://github.com/MWATelescope/mwa_pb/ but we don’t import
-    that module because it’s not python 3 compatible.
-
-    Note that the azimuth convention in for the UVBeam object is different than the
-    azimuth convention in the mwa_pb repo. In that repo, the azimuth convention is
-    changed from the native FEKO convention (the FEKO convention is the same as the
-    UVBeam convention). The convention in the mwa_pb repo has a different zero point
-    and a different direction (so it is in a left handed coordinate system).
+    This is based on https://github.com/MWATelescope/mwa_pb/ but we don't import
+    that module because it's not python 3 compatible.
+
+    Note that the azimuth convention for the UVBeam object is different than the
+    azimuth convention in the mwa_pb repo. In that repo, the azimuth convention
+    is changed from the native FEKO convention that the underlying data file is
+    in. The FEKO convention that the data file is in is the same as the UVBeam
+    convention, so we do not need to do a conversion here. The convention in the
+    mwa_pb repo is North through East, so it has a different zero point and a
+    different direction (so it is in a left handed coordinate system looking
+    down at the beam, a right handed coordinate system looking up at the sky).
 
     """
 
@@ -475,7 +478,10 @@ def _get_response(self, *, freqs_hz, pol_names, beam_modes, phi_arr, theta_arr):
                 Sigma_P = np.inner(phi_comp, emn_P_sum)
                 Sigma_T = np.inner(phi_comp, emn_T_sum)
 
-                jones[pol_i, 0, freq_i] = -Sigma_P
+                # we do not want a minus sign on Sigma_P unlike in mwa_pb because
+                # that minus sign is associated with the coordinate conversion
+                # they do that we do not want.
+                jones[pol_i, 0, freq_i] = Sigma_P
                 jones[pol_i, 1, freq_i] = Sigma_T
 
         return jones

src/pyuvdata/uvbeam/uvbeam.py

@@ -2064,7 +2064,7 @@ def interp(
             or for frequency only interpolation.
         reuse_spline : bool
             Save the interpolation functions for reuse. Only applies for
-            `az_za_simple` interpolation.
+            `az_za_simple` and `az_za_map_coordinates` interpolation.
         spline_opts : dict
             Provide options to numpy.RectBivariateSpline. This includes spline
             order parameters `kx` and `ky`, and smoothing parameter `s`.
@@ -2083,7 +2083,7 @@ def interp(
             intrinsic data array. Checking them can be quite computationally expensive.
             Conversely, if the passed az/za are outside of the domain, they will be
             silently extrapolated and the behavior is not well-defined. Only
-            applies for `az_za_simple` interpolation.
+            applies for `az_za_simple` and `az_za_map_coordinates` interpolation.
         return_basis_vector : bool
             Whether to return the interpolated basis vectors. Prior to v3.1.1 these
             were always returned. In v3.3+ they will _not_ be returned by default

tests/uvbeam/test_mwa_beam.py

@@ -113,6 +113,18 @@ def test_mwa_orientation(mwa_beam_1ppd):
     max_za_response = np.max(np.abs(efield_beam.data_array[1, north_ind, 0, za_val, :]))
     assert max_az_response > max_za_response
 
+    # check the sign of the responses are as expected close to zenith
+    efield_beam = mwa_beam_1ppd
+    za_val = np.nonzero(np.isclose(power_beam.axis2_array, 2.0 * np.pi / 180))
+
+    # first check zenith angle aligned response
+    assert efield_beam.data_array[1, east_ind, 0, za_val, east_az] > 0
+    assert efield_beam.data_array[1, north_ind, 0, za_val, north_az] > 0
+
+    # then check azimuthal aligned response
+    assert efield_beam.data_array[0, north_ind, 0, za_val, east_az] > 0
+    assert efield_beam.data_array[0, east_ind, 0, za_val, north_az] < 0
+
 
 def test_freq_range(mwa_beam_1ppd):
     beam1 = mwa_beam_1ppd