Resolve "Include LOFAR GriddedResponse"

cpp/common/system.h

+#ifndef EVERYBEAM_COMMON_SYSTEM_H_
+#define EVERYBEAM_COMMON_SYSTEM_H_
+
+#include <string>
+#include <sched.h>
+
+namespace everybeam {
+namespace common {
+
+class System {
+ public:
+  static std::size_t ProcessorCount() {
+#ifdef __APPLE__
+    return sysconf(_SC_NPROCESSORS_ONLN);
+#else
+    cpu_set_t cs;
+    CPU_ZERO(&cs);
+    sched_getaffinity(0, sizeof cs, &cs);
+
+    int count = 0;
+    for (int i = 0; i < CPU_SETSIZE; i++) {
+      if (CPU_ISSET(i, &cs)) count++;
+    }
+
+    return count;
+#endif
+  }
+};
+
+}  // namespace common
+}  // namespace everybeam
+
+#endif  // EVERYBEAM_COMMON_SYSTEM_H_
\ No newline at end of file

cpp/coords/coord_utils.h

@@ -7,17 +7,23 @@
 namespace everybeam {
 namespace coords {
 
+struct CoordinateSystem {
+  std::size_t width, height;
+  double ra, dec, dl, dm, phase_centre_dl, phase_centre_dm;
+};
+
 /**
  * @brief Convert Casacore itrfDir to vector3r_t
  *
- * @param itrfDir
+ * @param itrf_dir
  * @param itrf
  */
-void setITRFVector(const casacore::MDirection& itrfDir, vector3r_t& itrf) {
-  const casacore::Vector<double>& itrfVal = itrfDir.getValue().getValue();
-  itrf[0] = itrfVal[0];
-  itrf[1] = itrfVal[1];
-  itrf[2] = itrfVal[2];
+inline void SetITRFVector(const casacore::MDirection& itrf_dir,
+                          vector3r_t& itrf) {
+  const casacore::Vector<double>& itrf_val = itrf_dir.getValue().getValue();
+  itrf[0] = itrf_val[0];
+  itrf[1] = itrf_val[1];
+  itrf[2] = itrf_val[2];
 }
 }  // namespace coords
 }  // namespace everybeam

cpp/gridded_response/griddedresponse.h

@@ -24,7 +24,15 @@
 #ifndef EVERYBEAM_GRIDDEDRESPONSE_GRIDDEDRESPONSE_H_
 #define EVERYBEAM_GRIDDEDRESPONSE_GRIDDEDRESPONSE_H_
 
+#include "./../coords/coord_utils.h"
+#include "./../coords/ITRFDirection.h"
+#include "./../coords/ITRFConverter.h"
+
 #include <memory>
+#include <vector>
+#include <thread>
+#include <aocommon/lane.h>
+#include <casacore/measures/Measures/MDirection.h>
 
 namespace everybeam {
 
@@ -32,12 +40,6 @@ namespace telescope {
 class Telescope;
 }
 
-// TODO: just temporary location!
-struct CoordinateSystem {
-  std::size_t width, height;
-  double ra, dec, dl, dm, phaseCentreDL, phaseCentreDM;
-};
-
 namespace gridded_response {
 
 /**
@@ -46,42 +48,50 @@ namespace gridded_response {
  */
 class GriddedResponse {
  public:
-  typedef std::unique_ptr<GriddedResponse> Ptr;
-
-  // TODO: can be deprecated in a later stage
-  virtual void CalculateStation(std::size_t station_id) = 0;
-  virtual void CalculateStation() = 0;
-  virtual void CalculateAllStations() = 0;
+  /**
+   * @brief Compute the Beam response for a single station
+   *
+   * @param buffer Output buffer
+   * @param station_idx Station index, must be smaller than number of stations
+   * in the Telescope
+   * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+   * @param frequency Frequency (Hz)
+   */
+  virtual bool CalculateStation(std::complex<float>* buffer, double time,
+                                double freq, const size_t station_idx) = 0;
 
-  // TODO: complete!
-  virtual void CalculateStation(std::complex<float>* buffer, size_t station_id,
-                                double time, double freq) = 0;
-  // Repeated call of calculate single?
-  virtual void CalculateAllStations(std::complex<float>* buffer, double time,
-                                    double freq) = 0;
+  /**
+   * @brief Compute the Beam response for all stations in a Telescope
+   *
+   * @param buffer Output buffer
+   * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+   * @param frequency Frequency (Hz)
+   */
+  virtual bool CalculateAllStations(std::complex<float>* buffer, double time,
+                                    double frequency) = 0;
 
  protected:
   /**
    * @brief Construct a new Gridded Response object
    *
    * @param telescope_ptr Pointer to telescope::Telescope object
-   * @param coordinateSystem CoordinateSystem struct
+   * @param coordinate_system CoordinateSystem struct
    */
-  GriddedResponse(const telescope::Telescope* telescope_ptr,
-                  const CoordinateSystem& coordinateSystem)
-      : _telescope(telescope_ptr),
-        _width(coordinateSystem.width),
-        _height(coordinateSystem.height),
-        _ra(coordinateSystem.ra),
-        _dec(coordinateSystem.dec),
-        _dl(coordinateSystem.dl),
-        _dm(coordinateSystem.dm),
-        _phaseCentreDL(coordinateSystem.phaseCentreDL),
-        _phaseCentreDM(coordinateSystem.phaseCentreDM){};
+  GriddedResponse(const telescope::Telescope* const telescope_ptr,
+                  const coords::CoordinateSystem& coordinate_system)
+      : telescope_(telescope_ptr),
+        width_(coordinate_system.width),
+        height_(coordinate_system.height),
+        ra_(coordinate_system.ra),
+        dec_(coordinate_system.dec),
+        dl_(coordinate_system.dl),
+        dm_(coordinate_system.dm),
+        phase_centre_dl_(coordinate_system.phase_centre_dl),
+        phase_centre_dm_(coordinate_system.phase_centre_dm){};
 
-  const telescope::Telescope* _telescope;
-  size_t _width, _height;
-  double _ra, _dec, _dl, _dm, _phaseCentreDL, _phaseCentreDM;
+  const telescope::Telescope* const telescope_;
+  size_t width_, height_;
+  double ra_, dec_, dl_, dm_, phase_centre_dl_, phase_centre_dm_;
 };
 }  // namespace gridded_response
 }  // namespace everybeam

cpp/gridded_response/lofargrid.cc

+#include "lofargrid.h"
+#include "./../telescope/lofar.h"
+#include "./../common/system.h"
+
+#include <aocommon/imagecoordinates.h>
+#include <cmath>
+// #include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
+
+using namespace everybeam::gridded_response;
+
+LOFARGrid::LOFARGrid(telescope::Telescope* telescope_ptr,
+                     const coords::CoordinateSystem& coordinate_system)
+    : GriddedResponse(telescope_ptr, coordinate_system) {
+  const telescope::LOFAR& lofartelescope =
+      static_cast<const telescope::LOFAR&>(*telescope_);
+  size_t ncpus = common::System::ProcessorCount();
+
+  // Set private members
+  nthreads_ = std::min(ncpus, lofartelescope.nstations_);
+  threads_.resize(nthreads_);
+  delay_dir_ = lofartelescope.ms_properties_.delay_dir;
+  tile_beam_dir_ = lofartelescope.ms_properties_.tile_beam_dir;
+  subband_frequency_ = lofartelescope.ms_properties_.subband_freq;
+  use_channel_frequency_ = lofartelescope.options_.useChannelFrequency;
+};
+
+bool LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
+                                 double frequency, const size_t station_idx) {
+  aocommon::Lane<Job> lane(nthreads_);
+  lane_ = &lane;
+
+  SetITRFVectors(time);
+  double sb_freq = use_channel_frequency_ ? frequency : subband_frequency_;
+  // Dummy calculation of gain matrix, needed for multi-threading
+  telescope_->GetStation(station_idx)
+      ->Response(time, frequency, station0_, sb_freq, station0_, tile0_);
+
+  for (size_t i = 0; i != nthreads_; ++i) {
+    threads_[i] =
+        std::thread(&LOFARGrid::CalcThread, this, buffer, time, frequency);
+  }
+
+  for (size_t y = 0; y != height_; ++y) {
+    lane.write(Job{.y = y, .antenna_idx = station_idx});
+  }
+
+  lane.write_end();
+  for (size_t i = 0; i != nthreads_; ++i) threads_[i].join();
+  return true;
+}
+
+bool LOFARGrid::CalculateAllStations(std::complex<float>* buffer, double time,
+                                     double frequency) {
+  aocommon::Lane<Job> lane(nthreads_);
+  lane_ = &lane;
+
+  SetITRFVectors(time);
+
+  double sb_freq = use_channel_frequency_ ? frequency : subband_frequency_;
+
+  // Dummy loop, needed for multi-threading
+  for (size_t i = 0; i != telescope_->GetNrStations(); ++i) {
+    telescope_->GetStation(i)->Response(time, frequency, station0_, sb_freq,
+                                        station0_, tile0_);
+  }
+
+  // Prepare threads
+  for (size_t i = 0; i != nthreads_; ++i) {
+    threads_[i] =
+        std::thread(&LOFARGrid::CalcThread, this, buffer, time, frequency);
+  }
+
+  for (size_t y = 0; y != height_; ++y) {
+    for (size_t antenna_idx = 0; antenna_idx != telescope_->GetNrStations();
+         ++antenna_idx) {
+      lane.write(Job{.y = y, .antenna_idx = antenna_idx});
+    }
+  }
+
+  lane.write_end();
+  for (size_t i = 0; i != nthreads_; ++i) threads_[i].join();
+  return true;
+}
+
+void LOFARGrid::SetITRFVectors(double time) {
+  coords::ITRFConverter itrf_converter(time);
+  coords::SetITRFVector(itrf_converter.toDirection(delay_dir_), station0_);
+  coords::SetITRFVector(itrf_converter.toDirection(tile_beam_dir_), tile0_);
+
+  const casacore::Unit rad_unit("rad");
+
+  casacore::MDirection l_dir(
+      casacore::MVDirection(casacore::Quantity(ra_ + M_PI / 2, rad_unit),
+                            casacore::Quantity(0, rad_unit)),
+      casacore::MDirection::J2000);
+  coords::SetITRFVector(itrf_converter.toDirection(l_dir), l_vector_itrf_);
+
+  casacore::MDirection m_dir(
+      casacore::MVDirection(casacore::Quantity(ra_, rad_unit),
+                            casacore::Quantity(dec_ + M_PI / 2, rad_unit)),
+      casacore::MDirection::J2000);
+  coords::SetITRFVector(itrf_converter.toDirection(m_dir), m_vector_itrf_);
+
+  casacore::MDirection n_dir(
+      casacore::MVDirection(casacore::Quantity(ra_, rad_unit),
+                            casacore::Quantity(dec_, rad_unit)),
+      casacore::MDirection::J2000);
+  coords::SetITRFVector(itrf_converter.toDirection(n_dir), n_vector_itrf_);
+}
+
+void LOFARGrid::CalcThread(std::complex<float>* buffer, double time,
+                           double frequency) {
+  const size_t values_per_ant = width_ * height_ * 4;
+  double sb_freq = use_channel_frequency_ ? frequency : subband_frequency_;
+
+  Job job;
+  while (lane_->read(job)) {
+    for (size_t x = 0; x != width_; ++x) {
+      double l, m, n;
+      aocommon::ImageCoordinates::XYToLM(x, job.y, dl_, dm_, width_, height_, l,
+                                         m);
+      l += phase_centre_dl_;
+      m += phase_centre_dm_;
+      n = sqrt(1.0 - l * l - m * m);
+
+      vector3r_t itrf_direction;
+
+      itrf_direction[0] =
+          l * l_vector_itrf_[0] + m * m_vector_itrf_[0] + n * n_vector_itrf_[0];
+      itrf_direction[1] =
+          l * l_vector_itrf_[1] + m * m_vector_itrf_[1] + n * n_vector_itrf_[1];
+      itrf_direction[2] =
+          l * l_vector_itrf_[2] + m * m_vector_itrf_[2] + n * n_vector_itrf_[2];
+
+      std::complex<float>* base_buffer = buffer + (x + job.y * height_) * 4;
+
+      std::complex<float>* ant_buffer_ptr =
+          base_buffer + job.antenna_idx * values_per_ant;
+      matrix22c_t gain_matrix = telescope_->GetStation(job.antenna_idx)
+                                    ->Response(time, frequency, itrf_direction,
+                                               sb_freq, station0_, tile0_);
+
+      // (Optional) differential beam logic is handled at WSClean level
+      ant_buffer_ptr[0] = gain_matrix[0][0];
+      ant_buffer_ptr[1] = gain_matrix[0][1];
+      ant_buffer_ptr[2] = gain_matrix[1][0];
+      ant_buffer_ptr[3] = gain_matrix[1][1];
+    }
+  }
+}
\ No newline at end of file

cpp/gridded_response/lofargrid.h

@@ -26,7 +26,9 @@
 
 #include "griddedresponse.h"
 #include <iostream>
+#include <aocommon/matrix2x2.h>
 #include <complex>
+#include <limits>
 
 namespace everybeam {
 namespace gridded_response {
@@ -37,41 +39,17 @@ namespace gridded_response {
  */
 class LOFARGrid final : public GriddedResponse {
  public:
-  typedef std::unique_ptr<LOFARGrid> Ptr;
+  // typedef std::unique_ptr<LOFARGrid> Ptr;
 
   /**
    * @brief Construct a new LOFARGrid object
    *
    * @param telescope_ptr Pointer to telescope::LOFAR object
-   * @param coordinateSystem CoordinateSystem struct
+   * @param coordinate_system CoordinateSystem struct
    */
   LOFARGrid(telescope::Telescope* telescope_ptr,
-            const CoordinateSystem& coordinateSystem)
-      : GriddedResponse(telescope_ptr, coordinateSystem){};
+            const coords::CoordinateSystem& coordinate_system);
 
-  // TODO: remove this placeholders
-  void CalculateStation() { std::cout << "One is good" << std::endl; };
-
-  // TODO: remove this placeholder
-  void CalculateStation(std::size_t station_id) {
-    auto station_tmp = _telescope->GetStation(station_id);
-    std::cout << "Station name for index " << station_id << " is "
-              << station_tmp->name() << std::endl;
-  };
-
-  // TODO: remove this placeholder
-  void CalculateAllStations() {
-    std::size_t val = 0;
-    for (std::size_t station_id = 0; station_id < _telescope->stations.size();
-         ++station_id) {
-      val++;
-      // Repeated call to CalculateStation?
-      auto station_tmp = _telescope->GetStation(station_id);
-    };
-    std::cout << "I just read " << val << " stations" << std::endl;
-  };
-
-  // Geared towards implementation
   /**
    * @brief Compute the Beam response for a single station
    *
@@ -81,9 +59,9 @@ class LOFARGrid final : public GriddedResponse {
    * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
    * @param freq Frequency (Hz)
    */
-  void CalculateStation(std::complex<float>* buffer, size_t stationIdx,
-                        double time, double freq) override{};
-  // Repeated call of calculate single?
+  bool CalculateStation(std::complex<float>* buffer, double time,
+                        double frequency, const size_t station_idx) override;
+
   /**
    * @brief Compute the Beam response for all stations in a Telescope
    *
@@ -91,8 +69,26 @@ class LOFARGrid final : public GriddedResponse {
    * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
    * @param freq Frequency (Hz)
    */
-  void CalculateAllStations(std::complex<float>* buffer, double time,
-                            double freq) override{};
+  bool CalculateAllStations(std::complex<float>* buffer, double time,
+                            double frequency) override;
+
+ private:
+  casacore::MDirection delay_dir_, tile_beam_dir_;
+  vector3r_t station0_, tile0_, l_vector_itrf_, m_vector_itrf_, n_vector_itrf_;
+  bool use_channel_frequency_;
+  double subband_frequency_;
+
+  std::size_t nthreads_;
+  std::vector<std::thread> threads_;
+
+  struct Job {
+    size_t y, antenna_idx;
+  };
+  aocommon::Lane<Job>* lane_;
+
+  void SetITRFVectors(double time);
+
+  void CalcThread(std::complex<float>* buffer, double time, double frequency);
 };
 }  // namespace gridded_response
 }  // namespace everybeam

cpp/load.cc

@@ -15,14 +15,11 @@ TelescopeType Convert(const std::string &str) {
   else
     return UNKNOWN_TELESCOPE;
 }
-
-// Default options
-// Options defaultOptions;
 }  // namespace
 
-telescope::Telescope::Ptr Load(const casacore::MeasurementSet &ms,
-                               const ElementResponseModel model,
-                               const Options &options) {
+std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
+                                           const ElementResponseModel model,
+                                           const Options &options) {
   // Read Telescope name and convert to enum
   casacore::ScalarColumn<casacore::String> telescope_name_col(ms.observation(),
                                                               "TELESCOPE_NAME");
@@ -30,8 +27,9 @@ telescope::Telescope::Ptr Load(const casacore::MeasurementSet &ms,
   TelescopeType telescope_name = Convert(telescope_name_col(0));
   switch (telescope_name) {
     case LOFAR_TELESCOPE: {
-      auto telescope =
-          telescope::Telescope::Ptr(new telescope::LOFAR(ms, model, options));
+      std::unique_ptr<telescope::Telescope> telescope =
+          std::unique_ptr<telescope::Telescope>(
+              new telescope::LOFAR(ms, model, options));
       return telescope;
     }
     default:

cpp/load.h

@@ -39,8 +39,8 @@ namespace everybeam {
  * @param options Options
  * @return telescope::Telescope::Ptr
  */
-telescope::Telescope::Ptr Load(const casacore::MeasurementSet &ms,
-                               const ElementResponseModel model,
-                               const Options &options = Options::GetDefault());
+std::unique_ptr<telescope::Telescope> Load(
+    casacore::MeasurementSet &ms, const ElementResponseModel model,
+    const Options &options = Options::GetDefault());
 }  // namespace everybeam
 #endif  // EVERYBEAM_LOAD_H_
\ No newline at end of file

cpp/options.h

@@ -24,6 +24,8 @@
 #ifndef EVERYBEAM_OPTIONS_H_
 #define EVERYBEAM_OPTIONS_H_
 
+#include <string>
+
 namespace everybeam {
 
 /**
@@ -33,7 +35,7 @@ namespace everybeam {
  */
 class Options {
  public:
-  Options(){};
+  Options() : useDifferentialBeam(false), useChannelFrequency(true){};
 
   //! Default - empty - options class
   static Options GetDefault() { return Options(); };
@@ -41,8 +43,10 @@ class Options {
   // Scratch list of potential options
   // Specify path to element response coefficients file
   // std::string coeff_path;
-  // bool useDifferentialBeam
-  // bool useChannelFrequency
+  bool useDifferentialBeam;
+  bool useChannelFrequency;
+  std::string data_column_name;
+  vector3r_t diff_beam_centre;
   //
 };
 }  // namespace everybeam

cpp/telescope/lofar.cc

 #include "lofar.h"
+#include "./../gridded_response/lofargrid.h"
 #include "./../common/math_utils.h"
 #include "./../common/casa_utils.h"
+
+#include <aocommon/banddata.h>
 #include <cassert>
+#include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
 
 using namespace everybeam;
 using namespace everybeam::telescope;
@@ -113,16 +117,54 @@ vector3r_t readStationPhaseReference(const Table &table, unsigned int id) {
 }
 }  // namespace
 
-LOFAR::LOFAR(const MeasurementSet &ms, const ElementResponseModel model,
+LOFAR::LOFAR(MeasurementSet &ms, const ElementResponseModel model,
              const Options &options)
     : Telescope(ms, model, options) {
   ReadAllStations(ms, model);
+
+  // Populate MeasurementSet properties struct
+  aocommon::BandData band(ms.spectralWindow());
+  MSAntenna antenna(ms.antenna());
+  MPosition::ScalarColumn antenna_pos_col(
+      antenna, antenna.columnName(MSAntennaEnums::POSITION));
+  MEpoch::ScalarColumn time_column(ms,
+                                   ms.columnName(casacore::MSMainEnums::TIME));
+
+  // Following is ms.field() related, first check whether field complies with
+  // LOFAR field
+  if (ms.field().nrow() != 1)
+    throw std::runtime_error("LOFAR MeasurementSet has multiple fields");
+
+  if (!ms.field().tableDesc().isColumn("LOFAR_TILE_BEAM_DIR")) {
+    throw std::runtime_error("LOFAR_TILE_BEAM_DIR column not found");
+  }
+
+  casacore::ScalarMeasColumn<casacore::MDirection> delay_dir_col(
+      ms.field(),
+      casacore::MSField::columnName(casacore::MSFieldEnums::DELAY_DIR));
+
+  casacore::ArrayMeasColumn<casacore::MDirection> tile_beam_dir_col(
+      ms.field(), "LOFAR_TILE_BEAM_DIR");
+
+  // Populate struct
+  ms_properties_ = {.subband_freq = band.CentreFrequency(),
+                    .delay_dir = delay_dir_col(0),
+                    .tile_beam_dir = *(tile_beam_dir_col(0).data())};
 }
 
+std::unique_ptr<gridded_response::GriddedResponse> LOFAR::GetGriddedResponse(
+    const coords::CoordinateSystem &coordinate_system) {
+  // Get and return GriddedResponse ptr
+  std::unique_ptr<gridded_response::GriddedResponse> grid(
+      new gridded_response::LOFARGrid(this, coordinate_system));
+  // gridded_response::GriddedResponse grid(LOFARGrid(this, coordinate_system));
+  return grid;
+};
+
 Station::Ptr LOFAR::ReadStation(const MeasurementSet &ms, std::size_t id,
                                 const ElementResponseModel model) const {
   ROMSAntennaColumns antenna(ms.antenna());
-  assert(this->_nstations > id && !antenna.flagRow()(id));
+  assert(nstations_ > id && !antenna.flagRow()(id));
 
   // Get station name.
   const string name(antenna.name()(id));

cpp/telescope/lofar.h

@@ -26,16 +26,26 @@
 #define EVERYBEAM_TELESCOPE_LOFAR_H_
 
 #include "telescope.h"
-#include "./../gridded_response/lofargrid.h"
+
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MEpoch.h>
+#include <memory>
 
 namespace everybeam {
+
+namespace gridded_response {
+class LOFARGrid;
+class GriddedResponse;
+}  // namespace gridded_response
+
 namespace telescope {
 
 //! LOFAR telescope class
 class LOFAR final : public Telescope {
- public:
-  typedef std::unique_ptr<LOFAR> Ptr;
+  friend class gridded_response::LOFARGrid;
 
+ public:
   /**
    * @brief Construct a new LOFAR object
    *
@@ -43,22 +53,21 @@ class LOFAR final : public Telescope {
    * @param model Element Response model
    * @param options telescope options
    */
-  LOFAR(const casacore::MeasurementSet &ms, const ElementResponseModel model,
+  LOFAR(casacore::MeasurementSet &ms, const ElementResponseModel model,
         const Options &options);
 
-  // Docstrings will be inherited from telescope::Telescope
-  gridded_response::GriddedResponse::Ptr GetGriddedResponse(
-      const CoordinateSystem &coordinate_system) override {
-    // Get and return GriddedResponse ptr
-    gridded_response::GriddedResponse::Ptr grid(
-        new gridded_response::LOFARGrid(this, coordinate_system));
-    return grid;
-  };
+  std::unique_ptr<gridded_response::GriddedResponse> GetGriddedResponse(
+      const coords::CoordinateSystem &coordinate_system) override;
 
  private:
   Station::Ptr ReadStation(const casacore::MeasurementSet &ms,
                            const std::size_t id,
                            const ElementResponseModel model) const override;
+  struct MSProperties {
+    double subband_freq;
+    casacore::MDirection delay_dir, tile_beam_dir;
+  };
+  MSProperties ms_properties_;
 };
 }  // namespace telescope
 }  // namespace everybeam

cpp/telescope/telescope.h

@@ -27,15 +27,22 @@
 #include "./../station.h"
 #include "./../options.h"
 #include "./../element_response.h"
-#include "./../gridded_response/griddedresponse.h"
 
 #include <vector>
+#include <memory>
+#include <cassert>
 #include <casacore/ms/MeasurementSets/MeasurementSet.h>
 #include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
 
 namespace everybeam {
 
+namespace gridded_response {
+class GriddedResponse;
+}
+
+namespace coords {
 struct CoordinateSystem;
+}
 
 namespace telescope {
 
@@ -45,19 +52,14 @@ namespace telescope {
  */
 class Telescope {
  public:
-  typedef std::unique_ptr<Telescope> Ptr;
-
-  // Will be filled by the private read_all_stations() method
-  std::vector<Station::Ptr> stations;
-
   /**
    * @brief Return the gridded response object
    *
    * @param coordinate_system Coordinate system struct
    * @return GriddedResponse::Ptr
    */
-  virtual gridded_response::GriddedResponse::Ptr GetGriddedResponse(
-      const CoordinateSystem &coordinate_system) = 0;
+  virtual std::unique_ptr<gridded_response::GriddedResponse> GetGriddedResponse(
+      const coords::CoordinateSystem &coordinate_system) = 0;
 
   /**
    * @brief Get station by index
@@ -65,12 +67,14 @@ class Telescope {
    * @param station_id Station index to retrieve
    * @return Station::Ptr
    */
-  Station::Ptr GetStation(std::size_t station_id) const {
-    // TODO: throw exception if idx >_nstations
-
-    return stations[station_id];
+  Station::Ptr GetStation(std::size_t station_idx) const {
+    // Assert only in DEBUG mode
+    assert(station_idx < nstations_);
+    return stations_[station_idx];
   }
 
+  std::size_t GetNrStations() const { return nstations_; };
+
  protected:
   /**
    * @brief Construct a new Telescope object
@@ -79,10 +83,10 @@ class Telescope {
    * @param model ElementResponse model
    * @param options telescope options
    */
-  Telescope(const casacore::MeasurementSet &ms,
-            const ElementResponseModel model, const Options &options)
-      : _nstations(ms.antenna().nrow()), _options(options) {
-    stations.resize(_nstations);
+  Telescope(casacore::MeasurementSet &ms, const ElementResponseModel model,
+            const Options &options)
+      : nstations_(ms.antenna().nrow()), options_(options) {
+    stations_.resize(nstations_);
   };
 
   /**
@@ -97,12 +101,13 @@ class Telescope {
     casacore::ROMSAntennaColumns antenna(ms.antenna());
 
     for (std::size_t i = 0; i < antenna.nrow(); ++i) {
-      this->stations[i] = ReadStation(ms, i, model);
+      stations_[i] = ReadStation(ms, i, model);
     }
   };
 
-  std::size_t _nstations;
-  Options _options;
+  std::size_t nstations_;
+  Options options_;
+  std::vector<Station::Ptr> stations_;
 
  private:
   // Virtual method for reading a single telescope station

cpp/test/tload_lofar.cc

@@ -2,6 +2,7 @@
 
 #include "./../load.h"
 #include "./../options.h"
+#include "./../gridded_response/lofargrid.h"
 #include "./../element_response.h"
 
 #include "config.h"
@@ -21,19 +22,19 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
 
   // Assert if correct number of stations
   std::size_t nstations = 70;
-  BOOST_CHECK_EQUAL(telescope->stations.size(), nstations);
+  BOOST_CHECK_EQUAL(telescope->GetNrStations(), nstations);
 
   // Assert if GetStation(stationd_id) behaves properly
   BOOST_CHECK_EQUAL(telescope->GetStation(0)->name(), "CS001HBA0");
 
   // Get gridded response
-  CoordinateSystem coord_system;
+  coords::CoordinateSystem coord_system;
   auto grrp = telescope->GetGriddedResponse(coord_system);
   BOOST_CHECK(nullptr !=
               dynamic_cast<gridded_response::LOFARGrid*>(grrp.get()));
 
   // TODO: add test
-  //   grrp->CalculateStation(0);
+  // grrp->CalculateStation(0);
   //   grrp->CalculateStation();
   //   grrp->CalculateAllStations();
 }
\ No newline at end of file

demo/beam-helper.cpp

@@ -94,19 +94,19 @@ void GetITRFDirections(
                 casacore::Quantity(ra + M_PI/2, radUnit),
                 casacore::Quantity(0, radUnit)),
             	casacore::MDirection::J2000);
-            everybeam::coords::setITRFVector(itrfConverter.toDirection(lDir), _l_vector_itrf);
+            everybeam::coords::SetITRFVector(itrfConverter.toDirection(lDir), _l_vector_itrf);
 
             casacore::MDirection mDir(casacore::MVDirection(
                 casacore::Quantity(ra, radUnit),
                 casacore::Quantity(dec + M_PI/2, radUnit)),
             	casacore::MDirection::J2000);
-            everybeam::coords::setITRFVector(itrfConverter.toDirection(mDir), _m_vector_itrf);
+            everybeam::coords::SetITRFVector(itrfConverter.toDirection(mDir), _m_vector_itrf);
 
             casacore::MDirection nDir(casacore::MVDirection(
                 casacore::Quantity(ra, radUnit),
                 casacore::Quantity(dec, radUnit)),
             	casacore::MDirection::J2000);
-            everybeam::coords::setITRFVector(itrfConverter.toDirection(nDir), _n_vector_itrf);
+            everybeam::coords::SetITRFVector(itrfConverter.toDirection(nDir), _n_vector_itrf);
 
             vector3r_t itrfDirection;
 

external/README.md

@@ -3,7 +3,7 @@ External dependencies are included via git submodules, see .gitmodules in the ro
 Please note that the dependencies are fixed to a specific commit, to make sure tests do not break due to
 a silent updateing of the submodules. Dependencies are checked out on the following commits:
 
-- `aocommon`: be6072490469977f513722b87f3032812438e968
+- `aocommon`: 1f77c0c5e0d70507f227892353df8c5410cd87a4
 - `eigen`: `3.3.7` (https://gitlab.com/libeigen/eigen/-/tags/3.3.7)
 - `pybind11`: `v2.5.0` (https://github.com/pybind/pybind11/releases/tag/v2.5.0)
 

aocommon @ 1f77c0c5

-Subproject commit be6072490469977f513722b87f3032812438e968
+Subproject commit 1f77c0c5e0d70507f227892353df8c5410cd87a4