diff --git a/CMakeLists.txt b/CMakeLists.txt index 7f980325817019769fbdfc340c65ca107d34c816..f1f759986fc21548d73ec9e7ffb5580c71694875 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -700,6 +700,7 @@ if(BUILD_TESTING) include(CTest) set(TEST_FILENAMES + aartfaacreader/test/unit/tAntennaConfig.cc antennaflagger/test/unit/tFlagger.cc base/test/runtests.cc base/test/unit/tAartfaacSubtableWriter.cc diff --git a/aartfaacreader/AntennaConfig.h b/aartfaacreader/AntennaConfig.h new file mode 100644 index 0000000000000000000000000000000000000000..b11247645bf357008cb4b3c57f500e61bbe61a57 --- /dev/null +++ b/aartfaacreader/AntennaConfig.h @@ -0,0 +1,267 @@ +// Copyright (C) 2023 ASTRON (Netherlands Institute for Radio Astronomy) +// SPDX-License-Identifier: GPL-3.0-or-later + +/// This class is responsible to parse the LOFAR antenna configuration file. +/// Such configuration contains the antenna positions for a specific aartfaac +/// configration but also the orientation matrix of the common reference field + +/// This code was taken from +/// https://git.astron.nl/RD/aartfaac-tools/-/raw/master/lib/aartfaac2ms/antennaconfig.h?ref_type=heads +/// and adapted to fit into the DP3 codebase. + +#ifndef AARTFAACREADER_ANTENNACONFIG_H_ +#define AARTFAACREADER_ANTENNACONFIG_H_ + +#include <algorithm> +#include <array> +#include <cctype> +#include <fstream> +#include <locale> +#include <map> +#include <numeric> +#include <sstream> +#include <string> +#include <utility> +#include <vector> + +#include <base/RcuMode.h> +#include <casacore/measures/Measures/MPosition.h> + +namespace { + +// trim from start (in place) +void LTrim(std::string &s) { + s.erase(s.begin(), std::find_if(s.begin(), s.end(), + [](int ch) { return !std::isspace(ch); })); +} + +// trim from end (in place) +void RTrim(std::string &s) { + s.erase(std::find_if(s.rbegin(), s.rend(), + [](int ch) { return !std::isspace(ch); }) + .base(), + s.end()); +} + +// trim from both ends (in place) +void Trim(std::string &s) { + LTrim(s); + RTrim(s); +} +} // namespace + +namespace dp3::aartfaacreader { + +class AntennaConfig { + public: + explicit AntennaConfig(const char *filename) : file_(filename) { + ParseFile(); + } + + std::vector<casacore::MPosition> GetLBAPositions() const { + return GetPositions("LBA"); + } + + std::vector<casacore::MPosition> GetHBAPositions() const { + return GetPositions("HBA"); + } + + std::array<double, 9> GetLBAAxes() const { + return GetAxes("LBA_ROTATION_MATRIX"); + } + + std::array<double, 9> GetHBA0Axes() const { + return GetAxes("HBA0_ROTATION_MATRIX"); + } + + std::array<double, 9> GetHBA1Axes() const { + return GetAxes("HBA1_ROTATION_MATRIX"); + } + + std::array<double, 9> GetAxesFromMode(base::RcuMode mode) { + switch (mode.mode) { + case base::RcuMode::LBAInner10_90: + case base::RcuMode::LBAInner30_90: + case base::RcuMode::LBAOuter10_90: + case base::RcuMode::LBAOuter30_90: + return GetLBAAxes(); + case base::RcuMode::HBA110_190: + case base::RcuMode::HBA170_230: + case base::RcuMode::HBA210_270: + return GetHBA0Axes(); + break; + default: + throw std::runtime_error("Wrong RCU mode"); + } + } + + std::vector<casacore::MPosition> GetArrayFromMode(base::RcuMode mode) { + switch (mode.mode) { + case base::RcuMode::LBAInner10_90: + case base::RcuMode::LBAInner30_90: + case base::RcuMode::LBAOuter10_90: + case base::RcuMode::LBAOuter30_90: + return GetLBAPositions(); + break; + case base::RcuMode::HBA110_190: + case base::RcuMode::HBA170_230: + case base::RcuMode::HBA210_270: + return GetHBAPositions(); + default: + throw std::runtime_error("Wrong RCU mode"); + } + } + + private: + struct Array { + std::string name; + std::string band; + std::vector<double> data; + }; + + void ParseFile() { + Next(); + struct Array a; + while (ReadArray(a.name, a.band, a.data)) { + std::string key; + if (a.band.empty()) + key = a.name; + else + key = a.band + "_" + a.name; + values_.insert(std::make_pair(key, a)); + } + } + const std::vector<double> &GetArray(const std::string &name) const { + return values_.find(name)->second.data; + } + std::vector<casacore::MPosition> GetPositions( + const std::string &arrayName) const { + const std::vector<double> &arr = GetArray(arrayName); + std::vector<casacore::MPosition> position; + for (size_t index = 0; index < arr.size(); index += 6) { + position.emplace_back(casacore::MPosition{ + casacore::MVPosition{arr[index], arr[index + 1], arr[index + 2]}, + casacore::MPosition::ITRF}); + } + return position; + } + + std::array<double, 9> GetAxes(const std::string &arrayName) const { + const std::vector<double> &arr = GetArray(arrayName); + if (arr.size() != 9) + throw std::runtime_error( + "The array for coordinate axes in the antenna " + "config file had an incorrect size"); + + std::array<double, 9> axes; + for (size_t index = 0; index < 9; ++index) axes[index] = arr[index]; + return axes; + } + + bool Next() { + if (line_.empty() || line_position_ >= line_.size()) { + do { + std::getline(file_, line_); + if (!file_) { + token_.clear(); + line_.clear(); + return false; + } + Trim(line_); + } while (line_.empty() || line_[0] == '#'); + line_position_ = 0; + } + size_t pos = line_.find_first_of(" \t", line_position_); + if (pos == std::string::npos) { + token_ = line_.substr(line_position_); + line_.clear(); + line_position_ = 0; + return true; + } else { + token_ = line_.substr(line_position_, pos - line_position_); + line_position_ = pos + 1; + while (line_position_ < line_.size() && + (line_[line_position_] == ' ' || line_[line_position_] == '\t')) + ++line_position_; + Trim(token_); + if (token_.empty()) + return Next(); + else + return true; + } + } + + std::vector<int> ReadDimensions() { + std::vector<int> dimensions = {std::atoi(token_.c_str())}; + do { + if (!Next()) + throw std::runtime_error( + "Antenna config file has bad format: expected dimensions"); + if (token_ == "x") { + if (!Next()) + throw std::runtime_error( + "Antenna config file has bad format: " + "expected another dimension after x"); + int dimension_value = std::atoi(token_.c_str()); + dimensions.push_back(dimension_value); + } else if (token_ != "[") { + throw std::runtime_error("Antenna config file has bad format"); + } + } while (token_ != "["); + return dimensions; + } + + const std::vector<double> ReadData(const std::vector<int> &dimensions) { + int count = std::accumulate(dimensions.begin(), dimensions.end(), 1, + [](int a, int b) { return a * b; }); + std::vector<double> values(count); + for (int i = 0; i != count; ++i) { + if (!Next()) { + throw std::runtime_error("Missing numbers"); + } + values[i] = std::atof(token_.c_str()); + } + Next(); // move TO ']' + return values; + } + + bool ReadArray(std::string &name, std::string &band, + std::vector<double> &values) { + values.clear(); + if (!std::isalpha(token_[0])) return false; + name = token_; + Next(); + + if (std::isalpha(token_[0])) { + band = token_; + Next(); + } else { + band.clear(); + } + + std::vector<int> dimensions1 = ReadDimensions(); + std::vector<double> data1 = ReadData(dimensions1); + if (Next() && token_[0] >= '0' && token_[0] <= '9') { + std::vector<int> dimensions2 = ReadDimensions(); + std::vector<double> data2 = ReadData(dimensions2); + Next(); // skip ']' + + values = std::move(data2); + for (size_t i = 0; i != values.size(); ++i) + values[i] += data1[i % data1.size()]; + } else { + values = std::move(data1); + } + return true; + } + + std::map<std::string, Array> values_; + + std::ifstream file_; + std::string line_; + size_t line_position_; + std::string token_; +}; +} // namespace dp3::aartfaacreader + +#endif // AARTFAACREADER_ANTENNACONFIG_H_ diff --git a/aartfaacreader/test/unit/tAntennaConfig.cc b/aartfaacreader/test/unit/tAntennaConfig.cc new file mode 100644 index 0000000000000000000000000000000000000000..0b76b3252e12974f0ec87dad549ce9a5987aae1c --- /dev/null +++ b/aartfaacreader/test/unit/tAntennaConfig.cc @@ -0,0 +1,233 @@ +#include <aartfaacreader/AntennaConfig.h> +#include <base/RcuMode.h> +#include <casacore/measures/Measures/MPosition.h> +#include <filesystem> +#include <fstream> +#include <array> +#include <boost/filesystem.hpp> // for the unique_path generation +#include <boost/test/unit_test.hpp> +#include <boost/test/data/test_case.hpp> + +namespace { +bool ComparePositions(const casacore::MPosition &left, + const casacore::MPosition &right, + const casacore::MPosition &reference) { + return left.getValue() == (right.getValue() + reference.getValue()); +} + +casacore::MPosition PositionFromITRFCoordinates(double x, double y, double z) { + return casacore::MPosition(casacore::MVPosition{x, y, z}, + casacore::MPosition::ITRF); +} + +void TestPosition(const casacore::MPosition &left, + const casacore::MPosition &right, + const casacore::MPosition &reference) { + BOOST_CHECK(ComparePositions(left, right, reference)); +} + +void TestArrayPositions( + const std::vector<casacore::MPosition> &positions, + const std::vector<casacore::MPosition> &expected_position, + const casacore::MPosition &reference_antenna) { + BOOST_CHECK_EQUAL(positions.size(), expected_position.size()); + for (size_t i = 0; i < expected_position.size(); i++) { + TestPosition(positions[i], expected_position[i], reference_antenna); + } +} + +const int kAxesDimension = 3 * 3; +void TestAxes(const std::array<double, kAxesDimension> &axes, + const std::array<double, kAxesDimension> &expected_axes) { + for (int i = 0; i < kAxesDimension; i++) { + BOOST_CHECK_EQUAL(axes[i], expected_axes[i]); + } +} +const std::string configuration_example = + "\ +#\n\ +# AntennaPositions for AARTFAAC-12 LBA_OUTER antennas\n\ +# ITRF2005 target_date = 2015.5\n\ +# Created: 2023-10-16 14:35:38\n\ +#\n\ +\n\ +NORMAL_VECTOR LBA\n\ +3 [ 0.598753 0.072099 0.797682 ]\n\ +\n\ +ROTATION_MATRIX LBA\n\ +3 x 3 [\n\ + -0.1195950000 -0.7919540000 0.5987530000 \n\ + 0.9928230000 -0.0954190000 0.0720990000 \n\ + 0.0000330000 0.6030780000 0.7976820000 \n\ +]\n\ +\n\ +LBA\n\ +3 [ 3826577.462000000 461022.624000000 5064892.526 ]\n\ +4 x 2 x 3 [\n\ + 2.38300 -17.79500 -0.18000 2.38300 -17.79500 -0.18000\n\ + 0.95600 -20.19400 1.10800 0.95600 -20.19400 1.10800\n\ + -10.83100 -14.47100 9.43800 -10.83100 -14.47100 9.43800\n\ + -15.87100 0.64500 11.85500 -15.87100 0.64500 11.85500\n\ +]\n\ +\n\ +HBA\n\ +3 [ 3826577.462000000 461022.624000000 5064892.526 ]\n\ +4 x 2 x 3 [\n\ + -2.38300 -17.79500 -0.18000 -2.38300 -17.79500 -0.18000\n\ + -0.95600 -20.19400 1.10800 -0.95600 -20.19400 1.10800\n\ + +10.83100 -14.47100 9.43800 +10.83100 -14.47100 9.43800\n\ + +15.87100 0.64500 11.85500 +15.87100 0.64500 11.85500\n\ +]\n\ +\n\ +NORMAL_VECTOR HBA0\n\ +3 [ 0.598753 0.072099 0.797682 ]\n\ +\n\ +ROTATION_MATRIX HBA0\n\ +3 x 3 [\n\ + -0.1195950000 -0.7919540000 0.5987530000 \n\ + 0.9928230000 -0.0954190000 0.0720990000 \n\ + 0.0000330000 0.6030780000 0.7976820000 \n\ +]\n\ +\n\ +HBA0\n\ +3 [ 0.000000000 0.000000000 0.000 ]\n\ +\n\ +NORMAL_VECTOR HBA1\n\ +3 [ 0.598753 0.072099 0.797682 ]\n\ +\n\ +ROTATION_MATRIX HBA1\n\ +3 x 3 [\n\ + -0.1195950000 -0.7919540000 0.5987530000 \n\ + 0.9928230000 -0.0954190000 0.0720990000 \n\ + 0.0000330000 0.6030780000 0.7976820000 \n\ +]\n\ +\n\ +HBA1\n\ +3 [ 0.000000000 0.000000000 0.000 ]\n\ +\n"; + +struct AntennaConfigFixture { + AntennaConfigFixture() { + const std::filesystem::path tmp_path = + std::filesystem::temp_directory_path() / + boost::filesystem::unique_path("tmp%%%%%%%.config").string(); + path = tmp_path.string(); + std::ofstream config; + config.open(path); + config << configuration_example; + config.close(); + }; + + ~AntennaConfigFixture() { std::filesystem::remove(path); }; + + std::string path; +}; + +const casacore::MPosition kLbaReference{PositionFromITRFCoordinates( + 3826577.462000000, 461022.62400000, 5064892.526)}; + +const std::vector<casacore::MPosition> kLbaShifts = { + PositionFromITRFCoordinates(2.38300, -17.79500, -0.18000), + PositionFromITRFCoordinates(0.95600, -20.19400, 1.10800), + PositionFromITRFCoordinates(-10.83100, -14.47100, 9.43800), + PositionFromITRFCoordinates(-15.87100, 0.64500, 11.85500)}; + +const casacore::MPosition kHbaReference{PositionFromITRFCoordinates( + 3826577.462000000, 461022.62400000, 5064892.526)}; + +const std::vector<casacore::MPosition> kHbaShifts = { + PositionFromITRFCoordinates(-2.38300, -17.79500, -0.18000), + PositionFromITRFCoordinates(-0.95600, -20.19400, 1.10800), + PositionFromITRFCoordinates(+10.83100, -14.47100, 9.43800), + PositionFromITRFCoordinates(+15.87100, 0.64500, 11.85500)}; + +const std::array<double, kAxesDimension> kLbaAxes{ + -0.1195950000, -0.7919540000, 0.5987530000, 0.9928230000, -0.0954190000, + 0.0720990000, 0.0000330000, 0.6030780000, 0.7976820000}; + +const std::array<double, kAxesDimension> kHba0Axes{ + -0.1195950000, -0.7919540000, 0.5987530000, 0.9928230000, -0.0954190000, + 0.0720990000, 0.0000330000, 0.6030780000, 0.7976820000}; + +const std::array<double, kAxesDimension> kHba1Axes{ + -0.1195950000, -0.7919540000, 0.5987530000, 0.9928230000, -0.0954190000, + 0.0720990000, 0.0000330000, 0.6030780000, 0.7976820000}; + +const int kFirstLbaMode = 1; +const int kLastLbaMode = 4; +const int kFirstHbaMode = 5; +const int kLastHbaMode = 7; +} // namespace +BOOST_AUTO_TEST_SUITE(aartfaacantennaconfig) + +BOOST_FIXTURE_TEST_CASE(constructor_and_parsing, AntennaConfigFixture) { + BOOST_REQUIRE_NO_THROW( + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str())); +}; + +BOOST_FIXTURE_TEST_CASE(get_lba_positions, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + std::vector<casacore::MPosition> positions = antennaConfig.GetLBAPositions(); + TestArrayPositions(positions, kLbaShifts, kLbaReference); +} + +BOOST_FIXTURE_TEST_CASE(get_hba_positions, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + std::vector<casacore::MPosition> positions = antennaConfig.GetHBAPositions(); + TestArrayPositions(positions, kHbaShifts, kHbaReference); +} +BOOST_FIXTURE_TEST_CASE(get_lba_axes, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + std::array<double, kAxesDimension> lba_axes = antennaConfig.GetLBAAxes(); + TestAxes(lba_axes, kLbaAxes); +} +BOOST_FIXTURE_TEST_CASE(get_hba0_axes, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + std::array<double, kAxesDimension> hba0_axes = antennaConfig.GetHBA0Axes(); + TestAxes(hba0_axes, kHba0Axes); +} +BOOST_FIXTURE_TEST_CASE(get_hba1_axes, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + std::array<double, kAxesDimension> hba1_axes = antennaConfig.GetHBA1Axes(); + TestAxes(hba1_axes, kHba1Axes); +} + +BOOST_FIXTURE_TEST_CASE(get_axes_from_mode, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + for (int i = kFirstLbaMode; i <= kLastLbaMode; i++) { + TestAxes(antennaConfig.GetAxesFromMode(dp3::base::RcuMode::FromNumber(i)), + kLbaAxes); + } + for (int i = kFirstHbaMode; i <= kLastHbaMode; i++) { + TestAxes(antennaConfig.GetAxesFromMode(dp3::base::RcuMode::FromNumber(i)), + kHba0Axes); + } + + BOOST_CHECK_THROW( + antennaConfig.GetAxesFromMode(dp3::base::RcuMode::FromNumber(0)), + std::runtime_error); +} +BOOST_FIXTURE_TEST_CASE(get_array_from_mode, AntennaConfigFixture) { + dp3::aartfaacreader::AntennaConfig antennaConfig(path.c_str()); + + for (int i = kFirstLbaMode; i <= kLastLbaMode; i++) { + TestArrayPositions( + antennaConfig.GetArrayFromMode(dp3::base::RcuMode::FromNumber(i)), + kLbaShifts, kLbaReference); + } + for (int i = kFirstHbaMode; i <= kLastHbaMode; i++) { + TestArrayPositions( + antennaConfig.GetArrayFromMode(dp3::base::RcuMode::FromNumber(i)), + kHbaShifts, kHbaReference); + } + BOOST_CHECK_THROW( + antennaConfig.GetArrayFromMode(dp3::base::RcuMode::FromNumber(0)), + std::runtime_error); +} +BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file