diff --git a/CHANGELOG.md b/CHANGELOG.md index d1b2fcd53..445454b71 100644 --- a/CHANGELOG.md +++ b/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" diff --git a/docs/Images/airy_beam.png b/docs/Images/airy_beam.png index c4a3d7e31..258eec516 100644 Binary files a/docs/Images/airy_beam.png and b/docs/Images/airy_beam.png differ diff --git a/docs/Images/amplitude_waterfall.png b/docs/Images/amplitude_waterfall.png index 885345e29..495ecc378 100644 Binary files a/docs/Images/amplitude_waterfall.png and b/docs/Images/amplitude_waterfall.png differ diff --git a/docs/Images/dipole_mwa_power.png b/docs/Images/dipole_mwa_power.png index fdf3d7271..a46fcd902 100644 Binary files a/docs/Images/dipole_mwa_power.png and b/docs/Images/dipole_mwa_power.png differ diff --git a/docs/Images/short_dipole_beam.png b/docs/Images/short_dipole_beam.png index cd7e9fc06..27e291f4e 100644 Binary files a/docs/Images/short_dipole_beam.png and b/docs/Images/short_dipole_beam.png differ diff --git a/docs/analytic_beam_tutorial.rst b/docs/analytic_beam_tutorial.rst index 1e7cf0c9f..a4d4b3cce 100644 --- a/docs/analytic_beam_tutorial.rst +++ b/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 diff --git a/docs/beam_interface_tutorial.rst b/docs/beam_interface_tutorial.rst index 7788a5943..86cc2bc39 100644 --- a/docs/beam_interface_tutorial.rst +++ b/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') diff --git a/docs/uvdata_tutorial.rst b/docs/uvdata_tutorial.rst index c7a1ca740..bac35351d 100644 --- a/docs/uvdata_tutorial.rst +++ b/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])) diff --git a/src/pyuvdata/beam_interface.py b/src/pyuvdata/beam_interface.py index 0bd0bde9a..23a4e0f2b 100644 --- a/src/pyuvdata/beam_interface.py +++ b/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 ------- diff --git a/src/pyuvdata/uvbeam/mwa_beam.py b/src/pyuvdata/uvbeam/mwa_beam.py index aac26638f..ae0246b05 100644 --- a/src/pyuvdata/uvbeam/mwa_beam.py +++ b/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 diff --git a/src/pyuvdata/uvbeam/uvbeam.py b/src/pyuvdata/uvbeam/uvbeam.py index b07a19e26..d8e52ccaf 100644 --- a/src/pyuvdata/uvbeam/uvbeam.py +++ b/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 diff --git a/tests/uvbeam/test_mwa_beam.py b/tests/uvbeam/test_mwa_beam.py index dea7eb500..0dfdaa33c 100644 --- a/tests/uvbeam/test_mwa_beam.py +++ b/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