Merge branch '50-mwa-primarybeam' into 'master'

Resolve "Speed up MWA primarybeam calculation"

Closes #50

See merge request !64

CMakeLists.txt

@@ -12,6 +12,7 @@ string(TOLOWER ${CMAKE_PROJECT_NAME} projectname )
 set(CMAKE_MODULE_PATH ${CMAKE_SOURCE_DIR}/CMake)
 
 # Configure directory for data files
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -Wall -DNDEBUG")
 set(CMAKE_INSTALL_DATA_DIR "${CMAKE_INSTALL_PREFIX}/share/${projectname}")
 message("Storing data in: " ${CMAKE_INSTALL_DATA_DIR})
 

cpp/griddedresponse/griddedresponse.h

@@ -148,9 +148,10 @@ class GriddedResponse {
   double ra_, dec_, dl_, dm_, phase_centre_dl_, phase_centre_dm_;
 
  private:
-  void MakeIntegratedSnapshot(std::vector<aocommon::HMC4x4>& matrices,
-                              double time, double frequency, size_t field_id,
-                              const double* baseline_weights_interval);
+  virtual void MakeIntegratedSnapshot(std::vector<aocommon::HMC4x4>& matrices,
+                                      double time, double frequency,
+                                      size_t field_id,
+                                      const double* baseline_weights_interval);
 };
 }  // namespace griddedresponse
 }  // namespace everybeam

cpp/griddedresponse/mwagrid.cc

@@ -9,6 +9,9 @@
 
 #include <memory>
 
+using aocommon::HMC4x4;
+using aocommon::MC2x2;
+using aocommon::UVector;
 using everybeam::griddedresponse::MWAGrid;
 using everybeam::mwabeam::TileBeam2016;
 
@@ -66,4 +69,36 @@ void MWAGrid::CalculateAllStations(std::complex<float>* buffer, double time,
     std::copy_n(buffer, width_ * height_ * 4,
                 buffer + i * width_ * height_ * 4);
   }
-};
\ No newline at end of file
+};
+
+void MWAGrid::MakeIntegratedSnapshot(std::vector<aocommon::HMC4x4>& matrices,
+                                     double time, double frequency,
+                                     size_t field_id,
+                                     const double* baseline_weights_interval) {
+  size_t nstations = telescope_->GetNrStations();
+  UVector<std::complex<float>> buffer_undersampled(
+      GetStationBufferSize(nstations));
+  CalculateAllStations(buffer_undersampled.data(), time, frequency, field_id);
+
+  // For MWA, we can simply weight a (time) snapshot with the accumulated
+  // baseline weights
+  size_t nbaselines = nstations * (nstations + 1) / 2;
+  double snapshot_weight = 0.;
+  for (size_t index = 0; index != nbaselines; ++index) {
+    snapshot_weight += baseline_weights_interval[index];
+  }
+
+  for (size_t y = 0; y != height_; ++y) {
+    for (size_t x = 0; x != width_; ++x) {
+      size_t offset = (y * width_ + x) * 4;
+      MC2x2 A(buffer_undersampled[offset], buffer_undersampled[offset + 1],
+              buffer_undersampled[offset + 2], buffer_undersampled[offset + 3]);
+
+      // Mueller matrix constant for all baselines, so just compute once for
+      // each individual pixel
+      matrices[y * width_ + x] =
+          HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A) *
+          snapshot_weight;
+    }
+  }
+}
\ No newline at end of file

cpp/griddedresponse/mwagrid.h

@@ -46,6 +46,11 @@ class MWAGrid final : public GriddedResponse {
 
  private:
   std::unique_ptr<everybeam::mwabeam::TileBeam2016> tile_beam_;
+
+  // Override MakeIntegrated snapshot for efficiency
+  virtual void MakeIntegratedSnapshot(
+      std::vector<aocommon::HMC4x4>& matrices, double time, double frequency,
+      size_t field_id, const double* baseline_weights_interval) override;
 };
 }  // namespace griddedresponse
 }  // namespace everybeam

cpp/mwabeam/tilebeam2016.cc

@@ -14,27 +14,6 @@ TileBeam2016::TileBeam2016(const double* delays, bool frequency_interpolation,
       mwa_longitude_(116.67081),
       mwa_height_(377.0) {}
 
-void TileBeam2016::ArrayResponse(casacore::MEpoch& time,
-                                 casacore::MPosition& array_position, double ra,
-                                 double dec, double frequency,
-                                 std::complex<double>* gain) {
-  casacore::MeasFrame frame(array_position, time);
-  const casacore::MDirection::Ref hadec_ref(casacore::MDirection::HADEC, frame);
-  const casacore::MDirection::Ref azelgeo_ref(casacore::MDirection::AZELGEO,
-                                              frame);
-  const casacore::MDirection::Ref j2000_ref(casacore::MDirection::J2000, frame);
-  casacore::MPosition wgs = casacore::MPosition::Convert(
-      array_position, casacore::MPosition::WGS84)();
-  double arr_lattitude =
-      wgs.getValue().getLat();  // ant1Pos.getValue().getLat();
-
-  casacore::MDirection::Convert j2000_to_hadec(j2000_ref, hadec_ref),
-      j2000_to_azelgeo(j2000_ref, azelgeo_ref);
-
-  ArrayResponse(ra, dec, j2000_ref, j2000_to_hadec, j2000_to_azelgeo,
-                arr_lattitude, frequency, gain);
-}
-
 void TileBeam2016::ArrayResponse(
     double ra, double dec, const casacore::MDirection::Ref& j2000_ref,
     casacore::MDirection::Convert& j2000_to_hadec,

cpp/mwabeam/tilebeam2016.h

@@ -21,28 +21,35 @@ class TileBeam2016 : public Beam2016Implementation {
                const std::string &coeff_path);
 
   /**
-   * @brief Array response
+   * @brief API method for computing MWA array response
    *
-   * @param time Time (MJD)
-   * @param array_position
    * @param ra right ascension (rad)
    * @param dec declination (rad)
+   * @param j2000_ref J2000 ref coordinates
+   * @param j2000_to_hadecref HADEC coordinates
+   * @param j2000_to_azelgeoref AZELGEO coordinates
+   * @param arr_lattitude Lattitude
    * @param frequency Frequency (Hz)
-   * @param gain
+   * @param gain Gain matrix
    */
-  void ArrayResponse(casacore::MEpoch &time,
-                     casacore::MPosition &array_position, double ra, double dec,
-                     double frequency, std::complex<double> *gain);
-
   void ArrayResponse(double ra, double dec,
                      const casacore::MDirection::Ref &j2000_ref,
                      casacore::MDirection::Convert &j2000_to_hadecref,
                      casacore::MDirection::Convert &j2000_to_azelgeoref,
-                     double arrLatitude, double frequency,
+                     double arr_lattitude, double frequency,
                      std::complex<double> *gain);
 
+  /**
+   * @brief Compute MWA array response in given zenith/azimuth direction
+   *
+   * @param zenith_angle Zenith angle (rad)
+   * @param azimuth Azimuthal angle (rad)
+   * @param frequency Frequency (Hz)
+   * @param gain Gain matrix
+   */
   void ArrayResponse(double zenith_angle, double azimuth, double frequency,
                      std::complex<double> *gain) {
+    // As yet, this conditional is effectively redundant
     if (frequency_interpolation_)
       GetInterpolatedResponse(azimuth, zenith_angle, frequency, gain);
     else
@@ -69,6 +76,7 @@ class TileBeam2016 : public Beam2016Implementation {
     GetTabulatedResponse(az, za, freq, result);
   }
 
+  // TODO: variables seem to be unused. Remove?
   const double mwa_lattitude_;  // Array latitude. degrees North
   const double mwa_longitude_;  // Array longitude. degrees East
   const double mwa_height_;     // Array altitude. meters above sea level