diff --git a/CMake/config.h.in b/CMake/config.h.in
index e2270e4ff5aabf245fde9cdd2c413e182ed4bcc8..10c2186ab94fb341d734c8a0aca6d37d74a1b0ce 100644
--- a/CMake/config.h.in
+++ b/CMake/config.h.in
@@ -14,6 +14,7 @@
#define GMRT_MOCK_MS "@GMRT_MOCK_MS@"
#define KL_SCREEN_PATH "@KL_SCREEN_PATH@"
#define LWA_MOCK_PATH "@LWA_MOCK_PATH@"
+#define DISH_MOCK_PATH "@DISH_MOCK_PATH@"
#define H5_SOLUTIONS_PATH "@H5_SOLUTIONS_PATH@"
#define SCREEN_FITTING_MS "@SCREEN_FITTING_MS@"
#define OSKAR_MOCK_MS "@OSKAR_MOCK_MS@"
diff --git a/cpp/circularsymmetric/CMakeLists.txt b/cpp/circularsymmetric/CMakeLists.txt
index 894a2fa8424fcc0a507ee0b140682338b713e601..f5d6611b5823d345f5d6d9173bcac49d92f2caec 100644
--- a/cpp/circularsymmetric/CMakeLists.txt
+++ b/cpp/circularsymmetric/CMakeLists.txt
@@ -1,6 +1,6 @@
# Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
# SPDX-License-Identifier: GPL-3.0-or-later
-install(FILES atcacoefficients.h gmrtcoefficients.h meerkatcoefficients.h
- vlacoefficients.h voltagepattern.h
+install(FILES atcacoefficients.h coefficients.h gmrtcoefficients.h
+ meerkatcoefficients.h vlacoefficients.h voltagepattern.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/circularsymmetric")
diff --git a/cpp/pointresponse/dishpoint.cc b/cpp/pointresponse/dishpoint.cc
index 114e2ac74eb4412b189b060a109335bfc11eba59..98c8a9baa4bba18b7184be8aac499d2c209891bc 100644
--- a/cpp/pointresponse/dishpoint.cc
+++ b/cpp/pointresponse/dishpoint.cc
@@ -1,38 +1,102 @@
-// Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
+// Copyright (C) 2024 ASTRON (Netherlands Institute for Radio Astronomy)
// SPDX-License-Identifier: GPL-3.0-or-later
#include "dishpoint.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include <aocommon/uvector.h>
+
#include "../telescope/dish.h"
#include "../circularsymmetric/voltagepattern.h"
#include "../circularsymmetric/vlacoefficients.h"
-
-#include <aocommon/uvector.h>
+#include "./../coords/itrfdirection.h"
+#include "./../coords/itrfconverter.h"
+#include "../common/mathutils.h"
namespace everybeam {
namespace pointresponse {
+
+DishPoint::DishPoint(const telescope::Dish& dish, double time)
+ : PointResponse(&dish, time),
+ dish_(dish),
+ pointing_(dish_.GetFieldPointingMDirection().at(0)){};
+
void DishPoint::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
double ra, double dec, double freq,
size_t /* station_idx */, size_t field_id) {
- const telescope::Dish& dish_telescope =
- static_cast<const telescope::Dish&>(GetTelescope());
-
double pdir_ra;
double pdir_dec;
- std::tie(pdir_ra, pdir_dec) = dish_telescope.GetFieldPointing()[field_id];
+ std::tie(pdir_ra, pdir_dec) = dish_.GetFieldPointing()[field_id];
const double max_radius_arc_min =
- dish_telescope.GetDishCoefficients()->MaxRadiusInArcMin();
+ dish_.GetDishCoefficients()->MaxRadiusInArcMin();
const double reference_frequency =
- dish_telescope.GetDishCoefficients()->ReferenceFrequency();
+ dish_.GetDishCoefficients()->ReferenceFrequency();
+
+ // The computation below is ineffecient just a single point,
+ // see https://git.astron.nl/RD/EveryBeam/-/merge_requests/335#note_88457
circularsymmetric::VoltagePattern vp(
- dish_telescope.GetDishCoefficients()->GetFrequencies(freq),
- max_radius_arc_min);
+ dish_.GetDishCoefficients()->GetFrequencies(freq), max_radius_arc_min);
const aocommon::UVector<double> coefs_vec =
- dish_telescope.GetDishCoefficients()->GetCoefficients(freq);
+ dish_.GetDishCoefficients()->GetCoefficients(freq);
vp.EvaluatePolynomial(coefs_vec, reference_frequency,
- dish_telescope.GetDishCoefficients()->AreInverted());
+ dish_.GetDishCoefficients()->AreInverted());
vp.Render(buffer, ra, dec, pdir_ra, pdir_dec, freq);
}
+aocommon::MC2x2 DishPoint::Response(BeamMode beam_mode, size_t station_idx,
+ double freq, const vector3r_t& direction,
+ std::mutex* mutex) {
+ std::complex<float> buffer[4];
+
+ if (HasTimeUpdate()) {
+ // Compute ITRF pointing
+ {
+ std::unique_lock<std::mutex> lock(*(mutex ? mutex : &mutex_));
+ const coords::ItrfConverter itrf_converter(GetIntervalMidPoint());
+ pointing_itrf_ = itrf_converter.ToItrf(pointing_);
+ }
+ ClearTimeUpdate();
+ }
+
+ // Compute pointing_direction_angle
+ double pointing_direction_angle =
+ acos(std::max(std::min(dot(pointing_itrf_, direction), 1.0), -1.0));
+
+ const double max_radius_arc_min =
+ dish_.GetDishCoefficients()->MaxRadiusInArcMin();
+ const double reference_frequency =
+ dish_.GetDishCoefficients()->ReferenceFrequency();
+
+ // The computation below is ineffecient just a single point,
+ // see https://git.astron.nl/RD/EveryBeam/-/merge_requests/335#note_88457
+ circularsymmetric::VoltagePattern vp(
+ dish_.GetDishCoefficients()->GetFrequencies(freq), max_radius_arc_min);
+ const aocommon::UVector<double> coefs_vec =
+ dish_.GetDishCoefficients()->GetCoefficients(freq);
+ vp.EvaluatePolynomial(coefs_vec, reference_frequency,
+ dish_.GetDishCoefficients()->AreInverted());
+
+ // Because the voltage pattern is circular symmetric, the result
+ // depends only on the angle between the pointing and the requested direction
+ // Calling Render with direction_ra=pointing_direction_angle,
+ // direction_dec=0.0 and pointing_ra=0.0, pointing_ra=dec=0.0 yields the
+ // correct result
+
+ vp.Render(buffer, pointing_direction_angle, 0.0, 0.0, 0.0, freq);
+
+ // Because of the conversion from std::complex<float> to std::complex<double>
+ // It is not possible to directly return the buffer by a
+ // return aocommon::MC2x2(buffer);
+ // Instead, cast the buffer element by element before returning
+ aocommon::MC2x2 result;
+ for (int i = 0; i < 4; i++) {
+ result[i] = buffer[i];
+ }
+ return result;
+}
+
void DishPoint::ResponseAllStations(BeamMode beam_mode,
std::complex<float>* buffer, double ra,
double dec, double freq, size_t field_id) {
diff --git a/cpp/pointresponse/dishpoint.h b/cpp/pointresponse/dishpoint.h
index 634c37487aca4660da16ba1802270dfb1818d890..41a8d382361f71b4c65601814d2ebf3b2f00b4b5 100644
--- a/cpp/pointresponse/dishpoint.h
+++ b/cpp/pointresponse/dishpoint.h
@@ -7,7 +7,10 @@
#ifndef EVERYBEAM_POINTRESPONSE_DISHPOINT_H_
#define EVERYBEAM_POINTRESPONSE_DISHPOINT_H_
+#include <casacore/measures/Measures/MDirection.h>
+
#include "pointresponse.h"
+#include "../telescope/dish.h"
namespace everybeam {
namespace pointresponse {
@@ -19,16 +22,25 @@ namespace pointresponse {
*/
class [[gnu::visibility("default")]] DishPoint : public PointResponse {
public:
- DishPoint(const telescope::Telescope* telescope_ptr, double time)
- : PointResponse(telescope_ptr, time){};
+ DishPoint(const telescope::Dish& dish, double time);
void Response(BeamMode beam_mode, std::complex<float> * buffer, double ra,
double dec, double freq, size_t station_idx, size_t field_id)
override;
+ aocommon::MC2x2 Response(BeamMode beam_mode, size_t station_idx, double freq,
+ const vector3r_t& direction,
+ std::mutex* mutex = nullptr) override;
+
void ResponseAllStations(BeamMode beam_mode, std::complex<float> * buffer,
double ra, double dec, double freq, size_t field_id)
final override;
+
+ private:
+ const telescope::Dish& dish_;
+ const casacore::MDirection& pointing_;
+ vector3r_t pointing_itrf_;
+ std::mutex mutex_;
};
} // namespace pointresponse
} // namespace everybeam
diff --git a/cpp/pointresponse/skamidpoint.cc b/cpp/pointresponse/skamidpoint.cc
index 47067cbb9ff7def2b95cc771651b5dd568dde71c..240ff1fd127327c7e271a45b39b359de4eee2a18 100644
--- a/cpp/pointresponse/skamidpoint.cc
+++ b/cpp/pointresponse/skamidpoint.cc
@@ -14,7 +14,7 @@ namespace pointresponse {
SkaMidPoint::SkaMidPoint(const telescope::Telescope* telescope_ptr, double time,
ElementResponseModel element_response_model)
- : DishPoint(telescope_ptr, time),
+ : PointResponse(telescope_ptr, time),
element_response_model_(element_response_model) {
const telescope::SkaMid& ska_mid =
static_cast<const telescope::SkaMid&>(GetTelescope());
diff --git a/cpp/pointresponse/skamidpoint.h b/cpp/pointresponse/skamidpoint.h
index 23a08f34ed9b00e7c1657bf2a1a8b980c22ac88b..96dd84e468b6991b55da3704434296b9d963efa3 100644
--- a/cpp/pointresponse/skamidpoint.h
+++ b/cpp/pointresponse/skamidpoint.h
@@ -7,7 +7,7 @@
#ifndef EVERYBEAM_POINTRESPONSE_SKAMIDPOINT_H_
#define EVERYBEAM_POINTRESPONSE_SKAMIDPOINT_H_
-#include "dishpoint.h"
+#include "pointresponse.h"
#include "../elementresponse.h"
namespace everybeam {
@@ -17,7 +17,7 @@ class SkaMidResponse;
namespace pointresponse {
-class SkaMidPoint final : public DishPoint {
+class SkaMidPoint final : public PointResponse {
public:
SkaMidPoint(const telescope::Telescope* telescope_ptr, double time,
ElementResponseModel element_response_model);
diff --git a/cpp/telescope/dish.cc b/cpp/telescope/dish.cc
index 92f37884db6ec327520030b7a889e42c23356059..0f4160a74b54a3b2c4329dacb5afa811c2bf3381 100644
--- a/cpp/telescope/dish.cc
+++ b/cpp/telescope/dish.cc
@@ -20,12 +20,16 @@ Dish::Dish(const casacore::MeasurementSet& ms,
casacore::MSField field_table = ms.field();
casacore::ArrayColumn<double> pointing_dir_col(
field_table, casacore::MSField::columnName(casacore::MSField::DELAY_DIR));
+ casacore::ScalarMeasColumn<casacore::MDirection> pointing_dir_measure_col(
+ field_table,
+ casacore::MSField::columnName(casacore::MSFieldEnums::DELAY_DIR));
for (std::size_t field_id = 0; field_id != field_table.nrow(); ++field_id) {
casacore::Array<double> pdir = pointing_dir_col(field_id);
double pdir_ra = *pdir.cbegin();
double pdir_dec = *(pdir.cbegin() + 1);
field_pointing_.emplace_back(pdir_ra, pdir_dec);
+ field_pointing_mdirection_.emplace_back(pointing_dir_measure_col(field_id));
}
SetIsTimeRelevant(false);
}
@@ -36,5 +40,5 @@ std::unique_ptr<GriddedResponse> Dish::GetGriddedResponse(
}
std::unique_ptr<PointResponse> Dish::GetPointResponse(double time) const {
- return std::make_unique<DishPoint>(this, time);
+ return std::make_unique<DishPoint>(*this, time);
}
diff --git a/cpp/telescope/dish.h b/cpp/telescope/dish.h
index cf55c0e3b9e3e7c9715b922c3b9c7892f7cf8a75..65cb866978e950c9d40f62cf2ae998056622099b 100644
--- a/cpp/telescope/dish.h
+++ b/cpp/telescope/dish.h
@@ -51,6 +51,10 @@ class [[gnu::visibility("default")]] Dish final : public Telescope {
return field_pointing_;
}
+ const std::vector<casacore::MDirection>& GetFieldPointingMDirection() const {
+ return field_pointing_mdirection_;
+ }
+
const circularsymmetric::Coefficients* GetDishCoefficients() const {
return dish_coefficients_.get();
}
@@ -59,6 +63,7 @@ class [[gnu::visibility("default")]] Dish final : public Telescope {
std::unique_ptr<circularsymmetric::Coefficients> dish_coefficients_;
/// Store ra, dec pointing per field id from measurement set
std::vector<std::pair<double, double>> field_pointing_;
+ std::vector<casacore::MDirection> field_pointing_mdirection_;
};
} // namespace telescope
} // namespace everybeam
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index b875f7094c4a542182b5d8ddcb1cda737199ef02..617ceb1d1a45df8aea1291aeddb1bf9ce61da13d 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -40,6 +40,9 @@ set(KL_SCREEN_PATH
set(LWA_MOCK_PATH
${DATA_DIR}/LWA_OVRO_MOCK.ms
CACHE PATH "")
+set(DISH_MOCK_PATH
+ ${DATA_DIR}/DISH_MOCK.ms
+ CACHE PATH "")
set(LWA_COEFF_PATH
${CMAKE_BINARY_DIR}/coeffs/lwa
CACHE PATH "")
@@ -51,6 +54,7 @@ set(LOBES_COEFF_PATH
set(TEST_FILENAMES
runtests.cc
tatermconfig.cc
+ tdish.cc
telementresponsefixeddirection.cc
tsphericalharmonicsresponse.cc
tgmrt.cc
@@ -106,6 +110,8 @@ add_custom_target(
SCREEN_FITTING.ms
COMMAND ${CMAKE_SOURCE_DIR}/scripts/download_ms.sh lwa-ovro.MS.tar.bz2
LWA_OVRO_MOCK.ms
+ COMMAND ${CMAKE_SOURCE_DIR}/scripts/download_ms.sh meerkat_test.tar.gz
+ DISH_MOCK.ms
COMMAND ${CMAKE_SOURCE_DIR}/scripts/download_compressed_file.sh
OVRO_LWA_FEKO_TEST.tar.gz OVRO_LWA_FEKO_TEST.h5
COMMAND ${CMAKE_SOURCE_DIR}/scripts/download_h5_solutions.sh
diff --git a/cpp/test/tdish.cc b/cpp/test/tdish.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c4738af64b221baed6a7f3a07355983b496bc7ec
--- /dev/null
+++ b/cpp/test/tdish.cc
@@ -0,0 +1,92 @@
+// Copyright (C) 2024 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include <boost/test/unit_test.hpp>
+
+#include "config.h"
+#include "../load.h"
+#include "../options.h"
+#include "../pointresponse/dishpoint.h"
+#include "../telescope/dish.h"
+#include "../coords/itrfconverter.h"
+
+namespace everybeam {
+namespace {
+// First time stamp in mock ms"
+const double kTime = 5068498314.005126;
+const double kFrequency = 8.56313e+08;
+const double kRa = 0.90848969;
+const double kDec = -0.48149271;
+
+} // namespace
+
+BOOST_AUTO_TEST_SUITE(dish)
+
+BOOST_AUTO_TEST_CASE(load_dish) {
+ everybeam::Options options;
+
+ casacore::MeasurementSet ms(DISH_MOCK_PATH);
+
+ std::unique_ptr<everybeam::telescope::Telescope> telescope =
+ everybeam::Load(ms, options);
+
+ // Check that we have an Dish pointer.
+ const everybeam::telescope::Dish* dish_telescope =
+ dynamic_cast<const everybeam::telescope::Dish*>(telescope.get());
+ BOOST_REQUIRE(dish_telescope);
+
+ // Assert if correct number of stations
+ BOOST_CHECK_EQUAL(dish_telescope->GetNrStations(), size_t{62});
+
+ // Assert that we have a dish point response
+ std::unique_ptr<everybeam::pointresponse::PointResponse> point_response =
+ dish_telescope->GetPointResponse(kTime);
+ everybeam::pointresponse::DishPoint* dish_point_response =
+ dynamic_cast<everybeam::pointresponse::DishPoint*>(point_response.get());
+ BOOST_REQUIRE(dish_point_response);
+
+ std::vector<std::vector<std::complex<float>>> reference_response = {
+ {{1.0, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {1.0, 0.0}},
+ {{0.382599, 0.0}, {0.0, 0.0}, {0.0, 0.0}, {0.382599, 0.0}}};
+
+ std::vector<std::pair<double, double>> offsets = {{0.0, 0.0}, {0.01, -0.02}};
+
+ for (size_t j = 0; j < offsets.size(); j++) {
+ // Check the two response functions
+ double ra = kRa + offsets[j].first;
+ double dec = kDec + offsets[j].second;
+ size_t station_id = 0;
+ size_t field_id = 0;
+
+ std::array<std::complex<float>, 4> point_response_buffer;
+ dish_point_response->Response(everybeam::BeamMode::kFull,
+ point_response_buffer.data(), ra, dec,
+ kFrequency, station_id, field_id);
+
+ for (std::size_t i = 0; i < 4; ++i) {
+ std::cout << point_response_buffer[i] << std::endl;
+ BOOST_CHECK(
+ std::abs(point_response_buffer[i] - reference_response[j][i]) < 1e-6);
+ }
+
+ casacore::MDirection pointing(casacore::Quantity(ra, "rad"),
+ casacore::Quantity(dec, "rad"),
+ casacore::MDirection::J2000);
+
+ const coords::ItrfConverter itrf_converter(kTime);
+ vector3r_t direction = itrf_converter.ToItrf(pointing);
+
+ aocommon::MC2x2 response = dish_point_response->Response(
+ everybeam::BeamMode::kFull, station_id, kFrequency, direction);
+
+ for (std::size_t i = 0; i < 4; ++i) {
+ std::cout << response[i] << " " << reference_response[j][i] << std::endl;
+ BOOST_CHECK(std::abs(response[i] - std::complex<double>(
+ reference_response[j][i])) < 1e-6);
+ }
+ }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+} // namespace everybeam