Add support for dish telescopes

base/PredictBuffer.h

@@ -19,7 +19,7 @@ namespace base {
 class PredictBuffer {
  public:
   void resize(size_t n_threads, size_t n_correlations, size_t n_channels,
-              size_t n_baselines, size_t n_stations, bool include_beam,
+              size_t n_baselines, size_t n_stations_beam, bool include_beam,
               bool full_beam) {
     if (include_beam) {
       // The full buffer is not used when Stokes I is used -- conditionally
@@ -32,9 +32,9 @@ class PredictBuffer {
 
       for (size_t i = 0; i != n_threads; ++i) {
         if (full_beam) {
-          full_beam_values_[i].resize(n_stations * n_channels);
+          full_beam_values_[i].resize(n_stations_beam * n_channels);
         } else {
-          scalar_beam_values_[i].resize(n_stations * n_channels);
+          scalar_beam_values_[i].resize(n_stations_beam * n_channels);
         }
       }
     }

base/Telescope.cc

@@ -3,8 +3,6 @@
 
 #include "Telescope.h"
 
-#include <EveryBeam/telescope/phasedarray.h>
-
 // Support both old and new return value types of PhasedArray::GetStation().
 // TODO(AST-1001): Remove support for the old type in 2023.
 namespace {
@@ -25,14 +23,14 @@ namespace base {
 std::vector<size_t> SelectStationIndices(
     const everybeam::telescope::Telescope& telescope,
     const std::vector<std::string>& station_names) {
+  if (IsDish(telescope)) {
+    // For a dish telescope, the beam of all stations is assumed to be identical
+    return {0};
+  }
+
+  // If the telescope is not a dish, it is a phased array
   auto phased_array =
       dynamic_cast<const everybeam::telescope::PhasedArray*>(&telescope);
-  if (phased_array == nullptr) {
-    throw std::runtime_error(
-        "Currently, only PhasedArray telescopes accepted as input, i.e. "
-        "OSKAR or LOFAR. Support for other telescopes may become available "
-        "soon.");
-  }
 
   // Copy only those stations for which the name matches.
   // Note: the order of the station names in both vectors match,

base/Telescope.h

@@ -5,6 +5,8 @@
 #define DP3_BASE_TELESCOPE_H_
 
 #include <EveryBeam/load.h>
+#include <EveryBeam/telescope/phasedarray.h>
+#include <EveryBeam/telescope/dish.h>
 
 namespace dp3 {
 namespace base {
@@ -24,10 +26,20 @@ inline std::unique_ptr<everybeam::telescope::Telescope> GetTelescope(
   return telescope;
 }
 
+inline bool IsPhasedArray(const everybeam::telescope::Telescope& telescope) {
+  auto phased_array =
+      dynamic_cast<const everybeam::telescope::PhasedArray*>(&telescope);
+  return phased_array != nullptr;
+}
+
+inline bool IsDish(const everybeam::telescope::Telescope& telescope) {
+  auto dish = dynamic_cast<const everybeam::telescope::Dish*>(&telescope);
+  return dish != nullptr;
+}
+
 /**
  * Find stations in a telescope by name and return their indices.
  * @param telescope The telescope, which contains antennae / stations.
- *        Only PhasedArray telescopes are currently supported.
  * @param station_names A list of station names. The order of the names in this
  *        list should match the order in which they occur in the telescope.
  * @return The indices corresponding to the station names. Because of the

docker/py_wheel.docker

@@ -8,7 +8,7 @@ ARG PYMINOR
 FROM quay.io/casacore/casacore:py${PYMAJOR}${PYMINOR}_master
 
 ARG AOFLAGGER_VERSION=3.2.0
-ARG EVERYBEAM_VERSION=0.5.6
+ARG EVERYBEAM_VERSION=0.5.7
 ARG HDF5_VERSION=1.12.2
 ARG FFTW_VERSION=3.3.8
 ARG LUA_VERSION=5.3.6

docker/ubuntu_20_04_base

@@ -2,7 +2,7 @@ FROM ubuntu:20.04
 
 # TODO: needs to be bumped before next DP3 release
 # ENV IDG_VERSION=0.8
-ENV EVERYBEAM_VERSION=v0.5.6
+ENV EVERYBEAM_VERSION=v0.5.7
 ENV IDG_VERSION=6b61c038883ad3f807d20047c4f9e1a1f0b8d98a
 ENV AOFLAGGER_VERSION=65d5fba4f4c12797386d3fd9cd76734956a8b233
 

docker/ubuntu_22_04_base

@@ -2,7 +2,7 @@ FROM ubuntu:22.04
 
 # TODO: needs to be bumped before next DP3 release
 # ENV IDG_VERSION=0.8
-ENV EVERYBEAM_VERSION=v0.5.6
+ENV EVERYBEAM_VERSION=v0.5.7
 ENV IDG_VERSION=6b61c038883ad3f807d20047c4f9e1a1f0b8d98a
 ENV AOFLAGGER_VERSION=65d5fba4f4c12797386d3fd9cd76734956a8b233
 

steps/ApplyBeam.cc

@@ -148,8 +148,8 @@ void ApplyBeam::updateInfo(const DPInfo& infoIn) {
         static_cast<int>(everybeam::CorrectionMode::kNone));
   }
 
-  const size_t nSt = info().nantenna();
-  const size_t nCh = info().nchan();
+  const size_t n_stations = info().nantenna();
+  const size_t n_channels = info().nchan();
 
   const size_t nThreads = aocommon::ThreadPool::GetInstance().NThreads();
   itsBeamValues.resize(nThreads);
@@ -161,7 +161,7 @@ void ApplyBeam::updateInfo(const DPInfo& infoIn) {
   telescopes_.resize(nThreads);
 
   for (size_t thread = 0; thread < nThreads; ++thread) {
-    itsBeamValues[thread].resize(nSt * nCh);
+    itsBeamValues[thread].resize(n_stations * n_channels);
     itsMeasFrames[thread].set(info().arrayPosCopy());
     itsMeasFrames[thread].set(
         MEpoch(MVEpoch(info().startTime() / 86400), MEpoch::UTC));
@@ -267,13 +267,12 @@ template <typename T>
 void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0,
                           float* weight0, const everybeam::vector3r_t& srcdir,
                           const everybeam::telescope::Telescope* telescope,
-                          std::vector<aocommon::MC2x2>& beamValues, bool invert,
-                          everybeam::CorrectionMode mode, bool doUpdateWeights,
-                          std::mutex* mutex) {
+                          std::vector<aocommon::MC2x2>& beam_values,
+                          bool invert, everybeam::CorrectionMode mode,
+                          bool doUpdateWeights, std::mutex* mutex) {
   // Get the beam values for each station.
-  const size_t nCh = info.chanFreqs().size();
-  const size_t nSt = beamValues.size() / nCh;
-  const size_t nBl = info.nbaselines();
+  const size_t n_channels = info.chanFreqs().size();
+  const size_t n_baselines = info.nbaselines();
 
   std::unique_ptr<everybeam::pointresponse::PointResponse> point_response =
       telescope->GetPointResponse(time);
@@ -281,26 +280,29 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0,
   const std::vector<size_t> station_indices =
       base::SelectStationIndices(*telescope, info.antennaNames());
 
+  const size_t n_stations = station_indices.size();
+
   // Apply the beam values of both stations to the ApplyBeamed data.
-  for (size_t ch = 0; ch < nCh; ++ch) {
+  for (size_t ch = 0; ch < n_channels; ++ch) {
     switch (mode) {
       case everybeam::CorrectionMode::kFull:
       case everybeam::CorrectionMode::kElement:
-        // Fill beamValues for channel ch
-        for (size_t st = 0; st < nSt; ++st) {
-          beamValues[nCh * st + ch] = point_response->Response(
+        // Fill beam_values for channel ch
+        for (size_t st = 0; st < n_stations; ++st) {
+          beam_values[n_channels * st + ch] = point_response->Response(
               mode, station_indices[st], info.chanFreqs()[ch], srcdir, mutex);
           if (invert) {
             // Terminate if the matrix is not invertible.
-            [[maybe_unused]] bool status = beamValues[nCh * st + ch].Invert();
+            [[maybe_unused]] bool status =
+                beam_values[n_channels * st + ch].Invert();
             assert(status);
           }
         }
         break;
       case everybeam::CorrectionMode::kArrayFactor: {
         aocommon::MC2x2 af_tmp;
-        // Fill beamValues for channel ch
-        for (size_t st = 0; st < nSt; ++st) {
+        // Fill beam_values for channel ch
+        for (size_t st = 0; st < n_stations; ++st) {
           af_tmp = point_response->Response(
               mode, station_indices[st], info.chanFreqs()[ch], srcdir, mutex);
 
@@ -308,13 +310,13 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0,
             af_tmp[0] = 1. / af_tmp[0];
             af_tmp[3] = 1. / af_tmp[3];
           }
-          beamValues[nCh * st + ch] = af_tmp;
+          beam_values[n_channels * st + ch] = af_tmp;
         }
         break;
       }
       case everybeam::CorrectionMode::kNone:  // this should not happen
-        for (size_t st = 0; st < nSt; ++st) {
-          beamValues[nCh * st + ch] = aocommon::MC2x2::Unity();
+        for (size_t st = 0; st < n_stations; ++st) {
+          beam_values[n_channels * st + ch] = aocommon::MC2x2::Unity();
         }
         break;
     }
@@ -322,15 +324,23 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time, T* data0,
     // Apply beam for channel ch on all baselines
     // For mode=ARRAY_FACTOR, too much work is done here because we know
     // that r and l are diagonal
-    for (size_t bl = 0; bl < nBl; ++bl) {
-      T* data = data0 + bl * 4 * nCh + ch * 4;
+    for (size_t bl = 0; bl < n_baselines; ++bl) {
+      T* data = data0 + bl * 4 * n_channels + ch * 4;
       const aocommon::MC2x2F mat(data);
-      const aocommon::MC2x2F left(beamValues[nCh * info.getAnt1()[bl] + ch]);
-      const aocommon::MC2x2F right(beamValues[nCh * info.getAnt2()[bl] + ch]);
+
+      // If the beam is the same for all stations (i.e. when n_stations = 1),
+      // all baselines will have the same beam values
+      size_t index_left =
+          (n_stations == 1 ? ch : n_channels * info.getAnt1()[bl] + ch);
+      size_t index_right =
+          (n_stations == 1 ? ch : n_channels * info.getAnt2()[bl] + ch);
+      const aocommon::MC2x2F left(beam_values[index_left]);
+      const aocommon::MC2x2F right(beam_values[index_right]);
       const aocommon::MC2x2F result = left * mat.MultiplyHerm(right);
       result.AssignTo(data);
       if (doUpdateWeights) {
-        ApplyCal::ApplyWeights(left, right, weight0 + bl * 4 * nCh + ch * 4);
+        ApplyCal::ApplyWeights(left, right,
+                               weight0 + bl * 4 * n_channels + ch * 4);
       }
     }
   }
@@ -340,7 +350,7 @@ template void ApplyBeam::applyBeam(
     const DPInfo& info, double time, std::complex<double>* data0,
     float* weight0, const everybeam::vector3r_t& srcdir,
     const everybeam::telescope::Telescope* telescope,
-    std::vector<aocommon::MC2x2>& beamValues, bool invert,
+    std::vector<aocommon::MC2x2>& beam_values, bool invert,
     everybeam::CorrectionMode mode, bool doUpdateWeights, std::mutex* mutex);
 
 void ApplyBeam::applyBeam(const DPInfo& info, double time,
@@ -355,13 +365,13 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time,
                           bool do_update_weights, std::mutex* mutex) {
   // Get the beam values for each station.
   const size_t n_channels = info.chanFreqs().size();
-  const size_t n_stations = beam_values.size() / n_channels;
 
   std::unique_ptr<everybeam::pointresponse::PointResponse> point_response =
       telescope->GetPointResponse(time);
 
   const std::vector<size_t> station_indices =
       base::SelectStationIndices(*telescope, info.antennaNames());
+  const size_t n_stations = station_indices.size();
 
   // Apply the beam values of both stations to the ApplyBeamed data.
 
@@ -419,10 +429,13 @@ void ApplyBeam::applyBeam(const DPInfo& info, double time,
     for (size_t bl = baseline_range.first; bl < baseline_range.second; ++bl) {
       std::complex<double>* data = data0 + bl * 4 * n_channels + ch * 4;
       const aocommon::MC2x2F mat(data);
-      const aocommon::MC2x2F left(
-          beam_values[n_channels * info.getAnt1()[bl] + ch]);
-      const aocommon::MC2x2F right(
-          beam_values[n_channels * info.getAnt2()[bl] + ch]);
+      size_t index_left =
+          (n_stations == 1 ? ch : n_channels * info.getAnt1()[bl] + ch);
+      size_t index_right =
+          (n_stations == 1 ? ch : n_channels * info.getAnt2()[bl] + ch);
+      const aocommon::MC2x2F left(beam_values[index_left]);
+      const aocommon::MC2x2F right(beam_values[index_right]);
+
       const aocommon::MC2x2F result = left * mat.MultiplyHerm(right);
       result.AssignTo(data);
       if (do_update_weights) {
@@ -440,37 +453,42 @@ void ApplyBeam::applyBeamStokesIArrayFactor(
     const DPInfo& info, double time, T* data0,
     const everybeam::vector3r_t& srcdir,
     const everybeam::telescope::Telescope* telescope,
-    std::vector<everybeam::complex_t>& beamValues, bool invert,
+    std::vector<everybeam::complex_t>& beam_values, bool invert,
     everybeam::CorrectionMode mode, std::mutex* mutex) {
   assert(mode == everybeam::CorrectionMode::kArrayFactor);
   // Get the beam values for each station.
-  const size_t nCh = info.chanFreqs().size();
-  const size_t nSt = beamValues.size() / nCh;
-  const size_t nBl = info.nbaselines();
+  const size_t n_channels = info.chanFreqs().size();
+  const size_t n_baselines = info.nbaselines();
 
   std::unique_ptr<everybeam::pointresponse::PointResponse> point_response =
       telescope->GetPointResponse(time);
 
   const std::vector<size_t> station_indices =
       base::SelectStationIndices(*telescope, info.antennaNames());
+  const size_t n_stations = station_indices.size();
 
   // Apply the beam values of both stations to the ApplyBeamed data.
-  for (size_t ch = 0; ch < nCh; ++ch) {
-    // Fill beamValues for channel ch
-    for (size_t st = 0; st < nSt; ++st) {
-      beamValues[nCh * st + ch] = point_response->Response(
+  for (size_t ch = 0; ch < n_channels; ++ch) {
+    // Fill beam_values for channel ch
+    for (size_t st = 0; st < n_stations; ++st) {
+      beam_values[n_channels * st + ch] = point_response->Response(
           everybeam::BeamMode::kArrayFactor, station_indices[st],
           info.chanFreqs()[ch], srcdir, mutex)[0];
       if (invert) {
-        beamValues[nCh * st + ch] = 1. / beamValues[nCh * st + ch];
+        beam_values[n_channels * st + ch] =
+            1. / beam_values[n_channels * st + ch];
       }
     }
 
     // Apply beam for channel ch on all baselines
-    for (size_t bl = 0; bl < nBl; ++bl) {
-      T* data = data0 + bl * nCh + ch;
-      everybeam::complex_t* left = &(beamValues[nCh * info.getAnt1()[bl]]);
-      everybeam::complex_t* right = &(beamValues[nCh * info.getAnt2()[bl]]);
+    for (size_t bl = 0; bl < n_baselines; ++bl) {
+      T* data = data0 + bl * n_channels + ch;
+      size_t index_left =
+          (n_stations == 1 ? 0 : n_channels * info.getAnt1()[bl]);
+      size_t index_right =
+          (n_stations == 1 ? 0 : n_channels * info.getAnt2()[bl]);
+      everybeam::complex_t* left = &(beam_values[index_left]);
+      everybeam::complex_t* right = &(beam_values[index_right]);
       data[0] = left[ch] * std::complex<double>(data[0]) * conj(right[ch]);
 
       // TODO: update weights?
@@ -482,7 +500,7 @@ template void ApplyBeam::applyBeamStokesIArrayFactor(
     const DPInfo& info, double time, std::complex<double>* data0,
     const everybeam::vector3r_t& srcdir,
     const everybeam::telescope::Telescope* telescope,
-    std::vector<everybeam::complex_t>& beamValues, bool invert,
+    std::vector<everybeam::complex_t>& beam_values, bool invert,
     everybeam::CorrectionMode mode, std::mutex* mutex);
 
 void ApplyBeam::applyBeamStokesIArrayFactor(
@@ -496,13 +514,13 @@ void ApplyBeam::applyBeamStokesIArrayFactor(
   assert(mode == everybeam::CorrectionMode::kArrayFactor);
   // Get the beam values for each station.
   const size_t n_channels = info.chanFreqs().size();
-  const size_t n_stations = beam_values.size() / n_channels;
 
   std::unique_ptr<everybeam::pointresponse::PointResponse> point_response =
       telescope->GetPointResponse(time);
 
   const std::vector<size_t> station_indices =
       base::SelectStationIndices(*telescope, info.antennaNames());
+  const size_t n_stations = station_indices.size();
 
   // Apply the beam values of both stations to the ApplyBeamed data.
   if (station_range.first < n_stations) {
@@ -525,10 +543,12 @@ void ApplyBeam::applyBeamStokesIArrayFactor(
     // Apply beam for channel ch on the baselines handeled by this thread
     for (size_t bl = baseline_range.first; bl < baseline_range.second; ++bl) {
       std::complex<double>* data = data0 + bl * n_channels + ch;
-      everybeam::complex_t* left =
-          &(beam_values[n_channels * info.getAnt1()[bl]]);
-      everybeam::complex_t* right =
-          &(beam_values[n_channels * info.getAnt2()[bl]]);
+      size_t index_left =
+          (n_stations == 1 ? 0 : n_channels * info.getAnt1()[bl]);
+      size_t index_right =
+          (n_stations == 1 ? 0 : n_channels * info.getAnt2()[bl]);
+      everybeam::complex_t* left = &(beam_values[index_left]);
+      everybeam::complex_t* right = &(beam_values[index_right]);
       data[0] = left[ch] * std::complex<double>(data[0]) * conj(right[ch]);
 
       // TODO: update weights?

steps/OnePredict.cc

@@ -222,11 +222,14 @@ void OnePredict::initializeThreadData() {
   if (!predict_buffer_) {
     predict_buffer_ = std::make_shared<base::PredictBuffer>();
   }
+  bool is_dish_telescope = false;
   if (apply_beam_) {
     telescope_ = base::GetTelescope(info().msName(), element_response_model_,
                                     use_channel_freq_);
+    is_dish_telescope = base::IsDish(*telescope_);
   }
-  predict_buffer_->resize(nThreads, nCr, nCh, nBl, nSt, apply_beam_,
+  predict_buffer_->resize(nThreads, nCr, nCh, nBl,
+                          (is_dish_telescope ? 1 : nSt), apply_beam_,
                           !stokes_i_only_);
   // Create the Measure ITRF conversion info given the array position.
   // The time and direction are filled in later.

steps/test/integration/tApplyBeam.py

@@ -13,7 +13,7 @@ import sys
 sys.path.append(".")
 
 import testconfig as tcf
-from utils import assert_taql, untar_ms
+from utils import assert_taql, untar_ms, get_taql_result
 
 """
 Tests for applying the beam model.
@@ -25,9 +25,24 @@ Script can be invoked in two ways:
 """
 
 MSIN = "tNDPPP-generic.MS"
+DISH_MSIN = "tDish.MS"
 MSAPPLYBEAM = "tApplyBeam.tab"
 CWD = os.getcwd()
 
+"""
+The tDish.MS is a reduced version of a MEERKAT dataset, generated with the following command:
+
+DP3 \
+msin=MKT_CDFS_2_5_856_963MHz.ms \
+msout=tDish.MS \
+msout.overwrite=true \
+steps=[filter] \
+msin.starttime='29-Jun-2019/05:08:00.581' \
+msin.ntimes=5 \
+msin.nchan=5 \
+filter.blrange="[0,100,0,100]"
+"""
+
 
 @pytest.fixture(autouse=True)
 def source_env():
@@ -37,6 +52,7 @@ def source_env():
     os.chdir(tmpdir)
 
     untar_ms(f"{tcf.RESOURCEDIR}/{MSIN}.tgz")
+    untar_ms(f"{tcf.RESOURCEDIR}/{DISH_MSIN}.tgz")
     untar_ms(f"{tcf.SRCDIR}/{MSAPPLYBEAM}.tgz")
 
     # Tests are executed here
@@ -111,3 +127,46 @@ def test_with_updateweights():
     )
     taql_command = f"select from {MSIN} where all(near(WEIGHT_SPECTRUM, NEW_WEIGHT_SPECTRUM))"
     assert_taql(taql_command)
+
+
+from testconfig import TAQLEXE
+
+
+def test_dish_beam():
+    # Apply the beam with an offset of [0.01,-0.02] to the phase center
+    check_call(
+        [
+            tcf.DP3EXE,
+            f"msin={DISH_MSIN}",
+            f"msout=beam_applied.ms",
+            "msout.overwrite=true",
+            "steps=[applybeam]",
+            "applybeam.usechannelfreq=true",
+            "applybeam.direction=[0.91848969rad,-0.50149271rad]",
+        ]
+    )
+
+    # Assert that there are 25 elements which are zeroes before and after applying the beam
+    zero_elements = f"select t1.DATA[0,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[0,0])==0.0 AND abs(t2.DATA[0,0])==0.0"
+    assert_taql(zero_elements, 25)
+
+    # Assert that all other elements have changed after applying the beam
+    different_elements = f"select t1.DATA[0,0] from {DISH_MSIN} t1, beam_applied.ms t2 where (abs(t1.DATA[0,0]-t2.DATA[0,0])>1e-1) AND abs(t2.DATA[0,0])!=0.0"
+    assert_taql(different_elements, 425)
+
+    # Assert that the beam value is correct per each channel
+    assert_taql(
+        f"select t1.DATA[0,0]/t2.DATA[0,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[0,0] / t2.DATA[0,0] - 0.146382) > 1e-6 AND abs(t2.DATA[0,0])!=0.0"
+    )
+    assert_taql(
+        f"select t1.DATA[1,0]/t2.DATA[1,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[1,0] / t2.DATA[1,0] - 0.146003) > 1e-6 AND abs(t2.DATA[1,0])!=0.0"
+    )
+    assert_taql(
+        f"select t1.DATA[2,0]/t2.DATA[2,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[2,0] / t2.DATA[2,0] - 0.145628) > 1e-6 AND abs(t2.DATA[2,0])!=0.0"
+    )
+    assert_taql(
+        f"select t1.DATA[3,0]/t2.DATA[3,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[3,0] / t2.DATA[3,0] - 0.145258) > 1e-6 AND abs(t2.DATA[3,0])!=0.0"
+    )
+    assert_taql(
+        f"select t1.DATA[4,0]/t2.DATA[4,0] from {DISH_MSIN} t1, beam_applied.ms t2 where abs(t1.DATA[4,0] / t2.DATA[4,0] - 0.144933) > 1e-6 AND abs(t2.DATA[4,0])!=0.0"
+    )