From 121016a50c0c035cb1a91b84397532eed1dffcf5 Mon Sep 17 00:00:00 2001
From: Jakob Maljaars <jakob.maljaars@stcorp.nl>
Date: Mon, 13 Jul 2020 09:43:08 +0200
Subject: [PATCH] Resolve "Include LOFAR GriddedResponse"
---
cpp/common/system.h | 33 ++++++
cpp/coords/coord_utils.h | 18 ++-
cpp/gridded_response/griddedresponse.h | 76 +++++++------
cpp/gridded_response/lofargrid.cc | 150 +++++++++++++++++++++++++
cpp/gridded_response/lofargrid.h | 60 +++++-----
cpp/load.cc | 14 +--
cpp/load.h | 6 +-
cpp/options.h | 10 +-
cpp/telescope/lofar.cc | 46 +++++++-
cpp/telescope/lofar.h | 33 ++++--
cpp/telescope/telescope.h | 43 +++----
cpp/test/tload_lofar.cc | 7 +-
demo/beam-helper.cpp | 6 +-
external/README.md | 2 +-
external/aocommon | 2 +-
15 files changed, 380 insertions(+), 126 deletions(-)
create mode 100644 cpp/common/system.h
diff --git a/cpp/common/system.h b/cpp/common/system.h
new file mode 100644
index 00000000..8742f1c0
--- /dev/null
+++ b/cpp/common/system.h
@@ -0,0 +1,33 @@
+#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
diff --git a/cpp/coords/coord_utils.h b/cpp/coords/coord_utils.h
index 913abce8..0ce31f9b 100644
--- a/cpp/coords/coord_utils.h
+++ b/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
diff --git a/cpp/gridded_response/griddedresponse.h b/cpp/gridded_response/griddedresponse.h
index 2910f6c7..35abf9d5 100644
--- a/cpp/gridded_response/griddedresponse.h
+++ b/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
diff --git a/cpp/gridded_response/lofargrid.cc b/cpp/gridded_response/lofargrid.cc
index e69de29b..e42bfc40 100644
--- a/cpp/gridded_response/lofargrid.cc
+++ b/cpp/gridded_response/lofargrid.cc
@@ -0,0 +1,150 @@
+#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
diff --git a/cpp/gridded_response/lofargrid.h b/cpp/gridded_response/lofargrid.h
index 5fc79d12..ea524aa6 100644
--- a/cpp/gridded_response/lofargrid.h
+++ b/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
diff --git a/cpp/load.cc b/cpp/load.cc
index ab4c5ee8..afa4cf63 100644
--- a/cpp/load.cc
+++ b/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:
diff --git a/cpp/load.h b/cpp/load.h
index 2d8dd411..19102da3 100644
--- a/cpp/load.h
+++ b/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
diff --git a/cpp/options.h b/cpp/options.h
index 4f291bea..afd7629c 100644
--- a/cpp/options.h
+++ b/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
diff --git a/cpp/telescope/lofar.cc b/cpp/telescope/lofar.cc
index 0947cd9d..e1b02a6b 100644
--- a/cpp/telescope/lofar.cc
+++ b/cpp/telescope/lofar.cc
@@ -1,7 +1,11 @@
#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));
diff --git a/cpp/telescope/lofar.h b/cpp/telescope/lofar.h
index bf68d8bb..860f8566 100644
--- a/cpp/telescope/lofar.h
+++ b/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
diff --git a/cpp/telescope/telescope.h b/cpp/telescope/telescope.h
index 1782383f..c73eb508 100644
--- a/cpp/telescope/telescope.h
+++ b/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
diff --git a/cpp/test/tload_lofar.cc b/cpp/test/tload_lofar.cc
index 4b9de6ee..6babadd5 100644
--- a/cpp/test/tload_lofar.cc
+++ b/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
diff --git a/demo/beam-helper.cpp b/demo/beam-helper.cpp
index f2e5c444..26141de2 100644
--- a/demo/beam-helper.cpp
+++ b/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;
diff --git a/external/README.md b/external/README.md
index 0932f0e6..7cd736ab 100644
--- a/external/README.md
+++ b/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)
diff --git a/external/aocommon b/external/aocommon
index be607249..1f77c0c5 160000
--- a/external/aocommon
+++ b/external/aocommon
@@ -1 +1 @@
-Subproject commit be6072490469977f513722b87f3032812438e968
+Subproject commit 1f77c0c5e0d70507f227892353df8c5410cd87a4
--
GitLab