diff --git a/.gitignore b/.gitignore index bfa1a8d198bfb1d8df520af5b72dce63d60b5e4d..8c766b7bfe82d765e40d6f7b751e932ce24ecf1e 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ # Directories .vscode/ build/ +test_data/ diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index fe889e18e124518b0317a6928f0abd37f5900d43..ea054eb4f50d693eb7dee67f46e8d0311303dcaf 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -46,8 +46,13 @@ test-and-coverage: - mkdir -p /opt/everybeam/build - cd /opt/everybeam && git clone https://git.astron.nl/RD/EveryBeam.git EveryBeam - cd /opt/everybeam/EveryBeam && git checkout ${CI_COMMIT_SHORT_SHA} + # Download LOFAR Mock measurement set + - ./scripts/download_lofar_ms.sh + - LOFAR_MOCK_MS=/opt/everybeam/EveryBeam/test_data/LOFAR_MOCK.ms/ + # Change to build dir - cd /opt/everybeam/build - - cmake -DCMAKE_INSTALL_PREFIX=.. -DCMAKE_BUILD_TYPE=Debug -DCMAKE_CXX_FLAGS="-coverage" -DCMAKE_EXE_LINKER_FLAGS="-coverage" ../EveryBeam/ + # Build in Debug mode and add LOFAR_MOCK_MS + - cmake -DCMAKE_INSTALL_PREFIX=.. -DLOFAR_MOCK_MS=$LOFAR_MOCK_MS -DCMAKE_BUILD_TYPE=Debug -DCMAKE_CXX_FLAGS="-coverage" -DCMAKE_EXE_LINKER_FLAGS="-coverage" ../EveryBeam/ - make install -j8 - ctest -T test - gcovr -r .. -e '.*/external/.*' -e '.*/CompilerIdCXX/.*' -e '.*/test/.*' -e '.*/demo/.*' diff --git a/CMake/config.h.in b/CMake/config.h.in index 5a723152fcc198436161117ddbc88c1a4cfe28c6..7074966eb17c8e2639d9a5b1effdcf887ef385fc 100644 --- a/CMake/config.h.in +++ b/CMake/config.h.in @@ -3,5 +3,6 @@ #define EVERYBEAM_DATA_DIR "@CMAKE_INSTALL_DATA_DIR@" #define TEST_MEASUREMENTSET "@TEST_MEASUREMENTSET@" +#define LOFAR_MOCK_MS "@LOFAR_MOCK_MS@" #endif diff --git a/CMakeLists.txt b/CMakeLists.txt index a710cef4834667af96cf73f6bad50688fa73bcc4..d7819beab3902578098b8cbeb83cf86c10a445c3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -61,7 +61,11 @@ cmake_policy(SET CMP0074 NEW) endif() # Set compile options -add_compile_options(-std=c++11 "${OpenMP_CXX_FLAGS}" -Wall -DNDEBUG -Wl,--no-undefined) +add_compile_options(-std=c++11 "${OpenMP_CXX_FLAGS}" -Wall -Wl,--no-undefined) + +if (NOT CMAKE_BUILD_TYPE MATCHES Debug) + add_compile_options(-DNDEBUG) +endif() #------------------------------------------------------------------------------ # Find/load top level dependencies diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index 81b81be4dc69b0132e8834192608e77b5e40cd4d..e251a377a826d264eb93c7b022a30cbb2c606034 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -13,10 +13,13 @@ add_library(everybeam SHARED ElementResponse.cc BeamFormer.cc Element.cc + load.cc coords/ITRFConverter.cc coords/ITRFDirection.cc LofarMetaDataUtil.cc Station.cc + telescope/lofar.cc + gridded_response/lofargrid.cc ) target_include_directories(everybeam PUBLIC ${CASACORE_INCLUDE_DIR}) @@ -37,6 +40,10 @@ install (FILES ElementResponse.h LofarMetaDataUtil.h Station.h + # Related to new API: + load.h + # Not needed as yet: + # options.h DESTINATION "include/${CMAKE_PROJECT_NAME}") install( diff --git a/cpp/LofarMetaDataUtil.cc b/cpp/LofarMetaDataUtil.cc index 1a5572a3c350eefd954d9afc367a08c75491c109..4cb3ff8a0ea1b922cf1be99364abea2332b473e5 100644 --- a/cpp/LofarMetaDataUtil.cc +++ b/cpp/LofarMetaDataUtil.cc @@ -76,14 +76,17 @@ using namespace casacore; typedef std::array<vector3r_t, 16> TileConfig; +// TODO: can be deprecated in favor of common::hasColumn from CasaUtil.h bool hasColumn(const Table &table, const string &column) { return table.tableDesc().isColumn(column); } +// TODO: utility is not used at all bool hasSubTable(const Table &table, const string &name) { return table.keywordSet().isDefined(name); } +// TODO: can be deprecated in favor of common::getSubTable from CasaUtil.h Table getSubTable(const Table &table, const string &name) { return table.keywordSet().asTable(name); } diff --git a/cpp/Station.h b/cpp/Station.h index 978efe5e1f7b718cac1eea675774d06e9aecadb4..88ff762f284c720a3e35396abf567a487a821c3c 100644 --- a/cpp/Station.h +++ b/cpp/Station.h @@ -4,19 +4,20 @@ // ASTRON (Netherlands Institute for Radio Astronomy) // P.O.Box 2, 7990 AA Dwingeloo, The Netherlands // -// This file is part of the LOFAR software suite. -// The LOFAR software suite is free software: you can redistribute it and/or +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or // modify it under the terms of the GNU General Public License as published // by the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // -// The LOFAR software suite is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU General Public License for more details. // // You should have received a copy of the GNU General Public License along -// with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>. +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. // // $Id$ diff --git a/cpp/common/CMakeLists.txt b/cpp/common/CMakeLists.txt index 656f29ab61c6286a304b495f6ac8a13e5bcfc575..f115fb1aaefd58f78db90bc34bd99d11bf4d85cc 100644 --- a/cpp/common/CMakeLists.txt +++ b/cpp/common/CMakeLists.txt @@ -1,6 +1,7 @@ #------------------------------------------------------------------------------ # Install headers install (FILES + casa_utils.h Constants.h MathUtil.h MutablePtr.h diff --git a/cpp/common/casa_utils.h b/cpp/common/casa_utils.h new file mode 100644 index 0000000000000000000000000000000000000000..b0390d2940456c1655093664378ad1e0a02fd0f5 --- /dev/null +++ b/cpp/common/casa_utils.h @@ -0,0 +1,107 @@ +// casa_utils.h: CasaCore utilities. +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_COMMON_CASAUTIL_H_ +#define EVERYBEAM_COMMON_CASAUTIL_H_ + +#include "Types.h" +#include "./../Antenna.h" + +#include <cassert> + +#include <casacore/ms/MeasurementSets/MeasurementSet.h> +#include <casacore/ms/MeasurementSets/MSAntennaColumns.h> +#include <casacore/measures/Measures/MPosition.h> +#include <casacore/measures/Measures/MDirection.h> +#include <casacore/ms/MeasurementSets/MSAntennaColumns.h> +#include <casacore/measures/TableMeasures/ScalarMeasColumn.h> +#include <casacore/measures/Measures/MeasTable.h> +#include <casacore/ms/MSSel/MSSelection.h> +#include <casacore/tables/Tables/TableRecord.h> + +namespace everybeam { +namespace common { + +/** + * @brief Read coordinate system from MeasurementSet + * + * @param table Measurement set (casacore::Table) + * @param id Id of the antenna field in the station (int) + * @return Antenna::CoordinateSystem + */ +Antenna::CoordinateSystem readCoordinateSystem(const casacore::Table &table, + unsigned int id) { + casacore::ArrayQuantColumn<casacore::Double> c_position(table, "POSITION", + "m"); + casacore::ArrayQuantColumn<casacore::Double> c_axes(table, "COORDINATE_AXES", + "m"); + + // Read antenna field center (ITRF). + casacore::Vector<casacore::Quantity> aips_position = c_position(id); + assert(aips_position.size() == 3); + + vector3r_t position = {{aips_position(0).getValue(), + aips_position(1).getValue(), + aips_position(2).getValue()}}; + + // Read antenna field coordinate axes (ITRF). + casacore::Matrix<casacore::Quantity> aips_axes = c_axes(id); + assert(aips_axes.shape().isEqual(casacore::IPosition(2, 3, 3))); + + vector3r_t p = {{aips_axes(0, 0).getValue(), aips_axes(1, 0).getValue(), + aips_axes(2, 0).getValue()}}; + vector3r_t q = {{aips_axes(0, 1).getValue(), aips_axes(1, 1).getValue(), + aips_axes(2, 1).getValue()}}; + vector3r_t r = {{aips_axes(0, 2).getValue(), aips_axes(1, 2).getValue(), + aips_axes(2, 2).getValue()}}; + + Antenna::CoordinateSystem coordinate_system = {position, {p, q, r}}; + return coordinate_system; +} + +/** + * @brief Check if the specified column exists as a column of the + * specified table. + * + * @param table Measurement set (casacore::Table) + * @param column Column name (str) + * @return true If column present + * @return false If column not present + */ +bool hasColumn(const casacore::Table &table, const string &column) { + return table.tableDesc().isColumn(column); +} + +/** + * @brief Provide access to a sub-table by name. + * + * @param table Measurment set (casacore::Table) + * @param name Name of sub table (str) + * @return Table (casacore::Table) + */ +casacore::Table getSubTable(const casacore::Table &table, const string &name) { + return table.keywordSet().asTable(name); +} +} // namespace common +} // namespace everybeam +#endif // EVERYBEAM_COMMON_CASAUTIL_H_ \ No newline at end of file diff --git a/cpp/coords/CMakeLists.txt b/cpp/coords/CMakeLists.txt index 6022220e85a991796410b1c5c269e11a53dcf83b..34f2e0c3003883e8ea02ce7b8a26a2d85f0912f2 100644 --- a/cpp/coords/CMakeLists.txt +++ b/cpp/coords/CMakeLists.txt @@ -3,4 +3,5 @@ install (FILES ITRFConverter.h ITRFDirection.h + coord_utils.h DESTINATION "include/${CMAKE_PROJECT_NAME}/coords") diff --git a/cpp/coords/coord_utils.h b/cpp/coords/coord_utils.h new file mode 100644 index 0000000000000000000000000000000000000000..34d6955374b9203c18d11bf688333064dc09f8ad --- /dev/null +++ b/cpp/coords/coord_utils.h @@ -0,0 +1,24 @@ +#ifndef EVERYBEAM_COORDS_COORDUTILS_H_ +#define EVERYBEAM_COORDS_COORDUTILS_H_ + +#include "./../common/Types.h" +#include <casacore/measures/Measures/MCDirection.h> + +namespace everybeam { +namespace coords { + +/** + * @brief Convert Casacore itrfDir to vector3r_t + * + * @param itrfDir + * @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]; +} +} // namespace coords +} // namespace everybeam +#endif // EVERYBEAM_COORDS_COORDUTILS_H_ \ No newline at end of file diff --git a/cpp/gridded_response/griddedresponse.h b/cpp/gridded_response/griddedresponse.h new file mode 100644 index 0000000000000000000000000000000000000000..2910f6c7226ba31efca7bb776c4922fbb9d3b350 --- /dev/null +++ b/cpp/gridded_response/griddedresponse.h @@ -0,0 +1,88 @@ +// GriddedResponse.h: Base class for computing the (gridded) telescope response. +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_GRIDDEDRESPONSE_GRIDDEDRESPONSE_H_ +#define EVERYBEAM_GRIDDEDRESPONSE_GRIDDEDRESPONSE_H_ + +#include <memory> + +namespace everybeam { + +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 { + +/** + * @brief Virtual base class to compute the 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; + + // 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; + + protected: + /** + * @brief Construct a new Gridded Response object + * + * @param telescope_ptr Pointer to telescope::Telescope object + * @param coordinateSystem 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){}; + + const telescope::Telescope* _telescope; + size_t _width, _height; + double _ra, _dec, _dl, _dm, _phaseCentreDL, _phaseCentreDM; +}; +} // namespace gridded_response +} // namespace everybeam +#endif // EVERYBEAM_GRIDDEDRESPONSE_GRIDDEDRESPONSE_H_ \ No newline at end of file diff --git a/cpp/gridded_response/lofargrid.cc b/cpp/gridded_response/lofargrid.cc new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/cpp/gridded_response/lofargrid.h b/cpp/gridded_response/lofargrid.h new file mode 100644 index 0000000000000000000000000000000000000000..5fc79d1223a7009a6b26e18d105bb66537a569da --- /dev/null +++ b/cpp/gridded_response/lofargrid.h @@ -0,0 +1,99 @@ +// LOFARGrid.h: Class for computing the LOFAR (gridded) response. +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_GRIDDEDRESPONSE_LOFARGRID_H_ +#define EVERYBEAM_GRIDDEDRESPONSE_LOFARGRID_H_ + +#include "griddedresponse.h" +#include <iostream> +#include <complex> + +namespace everybeam { +namespace gridded_response { + +/** + * @brief Class for computing the LOFAR gridded response + * + */ +class LOFARGrid final : public GriddedResponse { + public: + typedef std::unique_ptr<LOFARGrid> Ptr; + + /** + * @brief Construct a new LOFARGrid object + * + * @param telescope_ptr Pointer to telescope::LOFAR object + * @param coordinateSystem CoordinateSystem struct + */ + LOFARGrid(telescope::Telescope* telescope_ptr, + const CoordinateSystem& coordinateSystem) + : GriddedResponse(telescope_ptr, coordinateSystem){}; + + // 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 + * + * @param buffer Output buffer + * @param stationIdx 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 freq Frequency (Hz) + */ + void CalculateStation(std::complex<float>* buffer, size_t stationIdx, + double time, double freq) override{}; + // Repeated call of calculate single? + /** + * @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 freq Frequency (Hz) + */ + void CalculateAllStations(std::complex<float>* buffer, double time, + double freq) override{}; +}; +} // namespace gridded_response +} // namespace everybeam +#endif // EVERYBEAM_GRIDDEDRESPONSE_LOFARGRID_H_ \ No newline at end of file diff --git a/cpp/load.cc b/cpp/load.cc new file mode 100644 index 0000000000000000000000000000000000000000..ab4c5ee87a2a5614f73519ee715560ab3316b5d9 --- /dev/null +++ b/cpp/load.cc @@ -0,0 +1,44 @@ +#include "load.h" +#include "options.h" + +#include <casacore/ms/MeasurementSets/MeasurementSet.h> +#include <casacore/tables/Tables/ScalarColumn.h> + +using namespace everybeam; +namespace everybeam { +namespace { +enum TelescopeType { UNKNOWN_TELESCOPE, LOFAR_TELESCOPE }; + +TelescopeType Convert(const std::string &str) { + if (str == "LOFAR") + return LOFAR_TELESCOPE; + else + return UNKNOWN_TELESCOPE; +} + +// Default options +// Options defaultOptions; +} // namespace + +telescope::Telescope::Ptr Load(const 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"); + + 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)); + return telescope; + } + default: + std::stringstream message; + message << "The requested telescope type " << telescope_name_col(0) + << " is not implemented."; + throw std::runtime_error(message.str()); + } +}; +} // namespace everybeam \ No newline at end of file diff --git a/cpp/load.h b/cpp/load.h new file mode 100644 index 0000000000000000000000000000000000000000..2d8dd4116a63e8055e856589dc6c4305f1810755 --- /dev/null +++ b/cpp/load.h @@ -0,0 +1,46 @@ +// load_telescope.h: Main interface function for loading a telescope +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_LOAD_H_ +#define EVERYBEAM_LOAD_H_ + +#include "./telescope/telescope.h" +#include "./telescope/lofar.h" +#include "options.h" + +namespace everybeam { + +/** + * @brief Load telescope given a measurement set. Telescope is determined + * from MeasurementSet meta-data. + * + * @param ms MeasurementSet + * @param model Element response model + * @param options Options + * @return telescope::Telescope::Ptr + */ +telescope::Telescope::Ptr Load(const 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 new file mode 100644 index 0000000000000000000000000000000000000000..4f291bea376dbd05bea3e94b747064664de8c139 --- /dev/null +++ b/cpp/options.h @@ -0,0 +1,49 @@ +// Options.h: Class for specifying the telescope response options. +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_OPTIONS_H_ +#define EVERYBEAM_OPTIONS_H_ + +namespace everybeam { + +/** + * @brief Class/Struct specifying everybeam Options. Needs further + * implementation! + * + */ +class Options { + public: + Options(){}; + + //! Default - empty - options class + static Options GetDefault() { return Options(); }; + + // Scratch list of potential options + // Specify path to element response coefficients file + // std::string coeff_path; + // bool useDifferentialBeam + // bool useChannelFrequency + // +}; +} // namespace everybeam +#endif // EVERYBEAM_OPTIONS_H_ \ No newline at end of file diff --git a/cpp/telescope/lofar.cc b/cpp/telescope/lofar.cc new file mode 100644 index 0000000000000000000000000000000000000000..87900f056f929595fa13615bcb00c8bbb5c588fe --- /dev/null +++ b/cpp/telescope/lofar.cc @@ -0,0 +1,172 @@ +#include "lofar.h" +#include "./../common/MathUtil.h" +#include "./../common/casa_utils.h" +#include <cassert> + +using namespace everybeam; +using namespace everybeam::telescope; +using namespace casacore; + +// LOFAR Telescope specific utils +namespace { + +constexpr Antenna::CoordinateSystem::Axes lofar_antenna_orientation = { + { + -std::sqrt(.5), + -std::sqrt(.5), + 0.0, + }, + { + std::sqrt(.5), + -std::sqrt(.5), + 0.0, + }, + {0.0, 0.0, 1.0}, +}; + +// TODO: move to common utils? +vector3r_t transformToFieldCoordinates( + const vector3r_t &position, const Antenna::CoordinateSystem::Axes &axes) { + const vector3r_t result{dot(position, axes.p), dot(position, axes.q), + dot(position, axes.r)}; + return result; +} + +// // Seems LOFAR specific? +// void transformToFieldCoordinates(TileConfig &config, +// const Antenna::CoordinateSystem::Axes &axes) +// { +// for (unsigned int i = 0; i < config.size(); ++i) { +// const vector3r_t position = config[i]; +// config[i][0] = dot(position, axes.p); +// config[i][1] = dot(position, axes.q); +// config[i][2] = dot(position, axes.r); +// } +// } + +// LOFAR specific? Or is this generic for each telescope? +BeamFormer::Ptr readAntennaField(const Table &table, std::size_t id, + ElementResponse::Ptr element_response) { + Antenna::CoordinateSystem coordinate_system = + common::readCoordinateSystem(table, id); + + BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system)); + + ROScalarColumn<String> c_name(table, "NAME"); + ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m"); + ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG"); + + const string &name = c_name(id); + + // Read element offsets and flags. + Matrix<Quantity> aips_offset = c_offset(id); + assert(aips_offset.shape().isEqual(IPosition(2, 3, aips_offset.ncolumn()))); + + Matrix<Bool> aips_flag = c_flag(id); + assert(aips_flag.shape().isEqual(IPosition(2, 2, aips_offset.ncolumn()))); + + // TileConfig tile_config; + // if(name != "LBA") readTileConfig(table, id); + // transformToFieldCoordinates(tile_config, coordinate_system.axes); + + for (size_t i = 0; i < aips_offset.ncolumn(); ++i) { + vector3r_t antenna_position = {aips_offset(0, i).getValue(), + aips_offset(1, i).getValue(), + aips_offset(2, i).getValue()}; + antenna_position = + transformToFieldCoordinates(antenna_position, coordinate_system.axes); + + Antenna::Ptr antenna; + Antenna::CoordinateSystem antenna_coordinate_system{ + antenna_position, lofar_antenna_orientation}; + if (name == "LBA") { + antenna = Element::Ptr( + new Element(antenna_coordinate_system, element_response, id)); + } else { + antenna = Element::Ptr( + new Element(antenna_coordinate_system, element_response, id)); + + // TODO + // HBA, HBA0, HBA1 + // antenna = make_tile(name, coordinate_system, tile_config, + // element_response); + } + + antenna->m_enabled[0] = !aips_flag(0, i); + antenna->m_enabled[1] = !aips_flag(1, i); + beam_former->add_antenna(antenna); + } + return beam_former; +} + +vector3r_t readStationPhaseReference(const Table &table, unsigned int id) { + vector3r_t phase_reference = {0.0, 0.0, 0.0}; + const string columnName("LOFAR_PHASE_REFERENCE"); + if (common::hasColumn(table, columnName)) { + ROScalarMeasColumn<MPosition> c_reference(table, columnName); + MPosition mReference = + MPosition::Convert(c_reference(id), MPosition::ITRF)(); + MVPosition mvReference = mReference.getValue(); + phase_reference = {mvReference(0), mvReference(1), mvReference(2)}; + } + return phase_reference; +} +} // namespace + +LOFAR::LOFAR(const MeasurementSet &ms, const ElementResponseModel model, + const Options &options) + : Telescope(ms, model, options) { + ReadAllStations(ms, model); +} + +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)); + + // Get station name. + const string name(antenna.name()(id)); + + // Get station position (ITRF). + MPosition mPosition = + MPosition::Convert(antenna.positionMeas()(id), MPosition::ITRF)(); + MVPosition mvPosition = mPosition.getValue(); + const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}}; + + // Create station. + Station::Ptr station(new Station(name, position, model)); + + // Read phase reference position (if available). + station->setPhaseReference(readStationPhaseReference(ms.antenna(), id)); + + Table tab_field = common::getSubTable(ms, "LOFAR_ANTENNA_FIELD"); + tab_field = tab_field(tab_field.col("ANTENNA_ID") == static_cast<Int>(id)); + + // The Station will consist of a BeamFormer that combines the fields + // coordinate system is ITRF + // phase reference is station position + auto beam_former = BeamFormer::Ptr(new BeamFormer( + Antenna::identity_coordinate_system, station->phaseReference())); + + for (std::size_t i = 0; i < tab_field.nrow(); ++i) { + beam_former->add_antenna( + readAntennaField(tab_field, i, station->get_element_response())); + } + + // TODO + // If There is only one field, the top level beamformer is not needed + // and the station antenna can be set the the beamformer of the field + station->set_antenna(beam_former); + + size_t field_id = 0; + size_t element_id = 0; + Antenna::CoordinateSystem coordinate_system = + common::readCoordinateSystem(tab_field, field_id); + + // TODO: rotate coordinate system for antenna + auto element = Element::Ptr(new Element( + coordinate_system, station->get_element_response(), element_id)); + station->set_element(element); + + return station; +} \ No newline at end of file diff --git a/cpp/telescope/lofar.h b/cpp/telescope/lofar.h new file mode 100644 index 0000000000000000000000000000000000000000..bf68d8bbb751606d0b4aa479850ce94bfd996ffe --- /dev/null +++ b/cpp/telescope/lofar.h @@ -0,0 +1,66 @@ +// LOFARTelescope.h: Base class for computing the response for the LOFAR +// telescope. +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_TELESCOPE_LOFAR_H_ +#define EVERYBEAM_TELESCOPE_LOFAR_H_ + +#include "telescope.h" +#include "./../gridded_response/lofargrid.h" + +namespace everybeam { +namespace telescope { + +//! LOFAR telescope class +class LOFAR final : public Telescope { + public: + typedef std::unique_ptr<LOFAR> Ptr; + + /** + * @brief Construct a new LOFAR object + * + * @param ms MeasurementSet + * @param model Element Response model + * @param options telescope options + */ + LOFAR(const 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; + }; + + private: + Station::Ptr ReadStation(const casacore::MeasurementSet &ms, + const std::size_t id, + const ElementResponseModel model) const override; +}; +} // namespace telescope +} // namespace everybeam + +#endif // EVERYBEAM_TELESCOPE_LOFAR_H_ \ No newline at end of file diff --git a/cpp/telescope/telescope.h b/cpp/telescope/telescope.h new file mode 100644 index 0000000000000000000000000000000000000000..cf023f55022d2c4d0cc5c64c89c3259cda750753 --- /dev/null +++ b/cpp/telescope/telescope.h @@ -0,0 +1,116 @@ +// Telescope.h: Base class for computing the Telescope response. +// +// Copyright (C) 2020 +// ASTRON (Netherlands Institute for Radio Astronomy) +// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands +// +// This file is part of the EveryBeam software suite. +// The EveryBeam software suite is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as published +// by the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// The EveryBeam software suite is distributed in the hope that it will be +// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License along +// with the EveryBeam software suite. If not, see +// <http://www.gnu.org/licenses/>. +// +// $Id$ + +#ifndef EVERYBEAM_TELESCOPE_TELESCOPE_H_ +#define EVERYBEAM_TELESCOPE_TELESCOPE_H_ + +#include "./../Station.h" +#include "./../options.h" +#include "./../ElementResponse.h" +#include "./../gridded_response/griddedresponse.h" + +#include <vector> +#include <casacore/ms/MeasurementSets/MeasurementSet.h> +#include <casacore/ms/MeasurementSets/MSAntennaColumns.h> + +namespace everybeam { + +struct CoordinateSystem; + +namespace telescope { + +/** + * @brief Telescope class, forming the base class for specific telescopes. + * + */ +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; + + /** + * @brief Get station by index + * + * @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]; + } + + protected: + /** + * @brief Construct a new Telescope object + * + * @param ms MeasurementSet + * @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); + }; + + /** + * @brief Read stations into vector + * + * @param out_it std::vector of stations + * @param ms measurement set + * @param model model + */ + void ReadAllStations(const casacore::MeasurementSet &ms, + const ElementResponseModel model) { + casacore::ROMSAntennaColumns antenna(ms.antenna()); + + for (std::size_t i = 0; i < antenna.nrow(); ++i) { + this->stations[i] = ReadStation(ms, i, model); + } + }; + + std::size_t _nstations; + Options _options; + + private: + // Virtual method for reading a single telescope station + virtual Station::Ptr ReadStation(const casacore::MeasurementSet &ms, + const std::size_t id, + const ElementResponseModel model) const = 0; +}; +} // namespace telescope +} // namespace everybeam + +#endif // EVERYBEAM_TELESCOPE_TELESCOPE_H_ \ No newline at end of file diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt index 25798968771317928e9d7159659a2f20cbeca2ac..bc662495c34c726ee1cecf36c84be7b0443857b6 100644 --- a/cpp/test/CMakeLists.txt +++ b/cpp/test/CMakeLists.txt @@ -1,16 +1,40 @@ +set(LOFAR_MOCK_MS CACHE STRING "LOFAR measurement set used for testing") + #------------------------------------------------------------------------------ -# Add executable(s) -add_executable(tstation tstation.cc) -target_link_libraries(tstation PUBLIC everybeam) -target_link_libraries(tstation PUBLIC OpenMP::OpenMP_CXX) +# Find boost +find_package(Boost COMPONENTS unit_test_framework) add_executable(taocommon taocommon.cc) add_executable(teigen teigen.cc) add_executable(tpybind11 tpybind11.cc) target_link_libraries(tpybind11 PRIVATE pybind11::embed) -#------------------------------------------------------------------------------ + # Add test -add_test(station-tests tstation) add_test(submodule-tests taocommon) add_test(submodule-tests teigen) -add_test(submodule-tests tpybind11) \ No newline at end of file +add_test(submodule-tests tpybind11) + +#------------------------------------------------------------------------------ +# Add tests if Boost found +if(Boost_FOUND) + include_directories(${Boost_INCLUDE_DIR}) + set(TEST_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/runtests.cc) + + # Test loading of LOFAR + if( LOFAR_MOCK_MS ) + set(TEST_SOURCES + ${TEST_SOURCES} + ${CMAKE_CURRENT_SOURCE_DIR}/tload_lofar.cc + ) + endif() + + # Add test executable + add_executable(runtests ${TEST_SOURCES}) + target_link_libraries(runtests everybeam ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) + # Required to get the config.h header + target_include_directories(runtests PRIVATE "${CMAKE_BINARY_DIR}") + + add_test(runtests runtests) +endif() + + diff --git a/cpp/test/runtests.cc b/cpp/test/runtests.cc new file mode 100644 index 0000000000000000000000000000000000000000..e93c2e3ae49574285f7ef558ef195a6b6bda042e --- /dev/null +++ b/cpp/test/runtests.cc @@ -0,0 +1,4 @@ +#define BOOST_TEST_MODULE everybeam +#define BOOST_TEST_DYN_LINK + +#include <boost/test/unit_test.hpp> \ No newline at end of file diff --git a/cpp/test/tload_lofar.cc b/cpp/test/tload_lofar.cc new file mode 100644 index 0000000000000000000000000000000000000000..23898aef60b78d622db53500d363bec99d882f6c --- /dev/null +++ b/cpp/test/tload_lofar.cc @@ -0,0 +1,39 @@ +#include <boost/test/unit_test.hpp> + +#include "./../load.h" +#include "./../options.h" +#include "./../ElementResponse.h" + +#include "config.h" + +using namespace everybeam; + +BOOST_AUTO_TEST_CASE(load_lofar) { + ElementResponseModel response_model = ElementResponseModel::Hamaker; + Options options; + casacore::MeasurementSet ms(LOFAR_MOCK_MS); + + // Load LOFAR Telescope + auto telescope = Load(ms, response_model, options); + + // Assert if we indeed have a LOFAR pointer + BOOST_CHECK(nullptr != dynamic_cast<telescope::LOFAR*>(telescope.get())); + + // Assert if correct number of stations + std::size_t nstations = 70; + BOOST_CHECK_EQUAL(telescope->stations.size(), nstations); + + // Assert if GetStation(stationd_id) behaves properly + BOOST_CHECK_EQUAL(telescope->GetStation(0)->name(), "CS001HBA0"); + + // Get gridded response + 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(); + // grrp->CalculateAllStations(); +} \ No newline at end of file diff --git a/cpp/test/tstation.cc b/cpp/test/tstation.cc index f9ced6450b06eb1af514de48b562589708243fc2..4d1198100b79532776e530672111015fe076ad31 100644 --- a/cpp/test/tstation.cc +++ b/cpp/test/tstation.cc @@ -1,5 +1,6 @@ #include "./../Station.h" +// TODO: make a test out of this int main() { const everybeam::vector3r_t position = {{1.0, 2.0, 3.0}}; diff --git a/demo/beam-helper.cpp b/demo/beam-helper.cpp index f5b804b51f00e811ad4ab8171bfeead093a1c5c9..45eeae300e4b961821a5a1764d91548e0a7296d3 100644 --- a/demo/beam-helper.cpp +++ b/demo/beam-helper.cpp @@ -1,5 +1,6 @@ #include "beam-helper.h" +#include <aocommon/imagecoordinates.h> #include <fitsio.h> #include "./../cpp/common/MathUtil.h" @@ -75,7 +76,7 @@ void GetITRFDirections( double l, m, n, ra, dec; - XYToLM<double>(x, y, dl, dm, subgrid_size, subgrid_size, l, m); + aocommon::ImageCoordinates::XYToLM<double>(x, y, dl, dm, subgrid_size, subgrid_size, l, m); l += pdl; m += pdm; @@ -93,19 +94,19 @@ void GetITRFDirections( casacore::Quantity(ra + M_PI/2, radUnit), casacore::Quantity(0, radUnit)), casacore::MDirection::J2000); - 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); - 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); - setITRFVector(itrfConverter.toDirection(nDir), _n_vector_itrf); + everybeam::coords::setITRFVector(itrfConverter.toDirection(nDir), _n_vector_itrf); vector3r_t itrfDirection; @@ -200,16 +201,6 @@ void StoreBeam( WriteFits<double>(filename, img.data(), nx*width, ny*height); } -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]; -} - void GetRaDecZenith( vector3r_t position, double time, diff --git a/demo/beam-helper.h b/demo/beam-helper.h index 4bfbab9e13b6694d2bf44515d2c160c94e7991ce..20b6fd93c9b7a5fa96e0ab1f2cf3f1c58aecde94 100644 --- a/demo/beam-helper.h +++ b/demo/beam-helper.h @@ -1,6 +1,7 @@ #include "./../cpp/ElementResponse.h" #include "./../cpp/Station.h" #include "./../cpp/LofarMetaDataUtil.h" +#include "./../cpp/coords/coord_utils.h" #include <casacore/measures/Measures/MPosition.h> #include <casacore/measures/Measures/MEpoch.h> @@ -27,8 +28,6 @@ void GetITRFDirections(vector3r_t* itrfDirections, size_t subgrid_size, void StoreBeam(const std::string& filename, const std::complex<float>* buffer, size_t nStations, size_t width, size_t height); -void setITRFVector(const casacore::MDirection& itrfDir, vector3r_t& itrf); - inline matrix22c_t operator*(const matrix22c_t& arg0, const matrix22c_t& arg1) { matrix22c_t result; result[0][0] = arg0[0][0] * arg1[0][0] + arg0[0][1] * arg1[1][0]; @@ -38,14 +37,6 @@ inline matrix22c_t operator*(const matrix22c_t& arg0, const matrix22c_t& arg1) { return result; } -template <typename T> -void XYToLM(size_t x, size_t y, T pixelSizeX, T pixelSizeY, size_t width, - size_t height, T& l, T& m) { - T midX = (T)width / 2.0, midY = (T)height / 2.0; - l = (midX - (T)x) * pixelSizeX; - m = ((T)y - midY) * pixelSizeY; -} - void GetRaDecZenith(vector3r_t position, double time, double& ra, double& dec); std::string GetFieldName(casacore::MeasurementSet& ms, diff --git a/demo/tStation.cc b/demo/tStation.cc index c6b6e3d311d1de19ab2083f284de97775e417337..e63cf65e72e7fa6c77c34d3ef4cd103917e7f132 100644 --- a/demo/tStation.cc +++ b/demo/tStation.cc @@ -1,10 +1,8 @@ #include <iostream> -#include "Station.h" - -// #include "BeamFormer.h" -#include "LofarMetaDataUtil.h" -#include "MathUtil.h" +#include "./../cpp/Station.h" +#include "./../cpp/LofarMetaDataUtil.h" +#include "./../cpp/common/MathUtil.h" #include "config.h" diff --git a/scripts/download_lofar_ms.sh b/scripts/download_lofar_ms.sh new file mode 100755 index 0000000000000000000000000000000000000000..487f878315e6d87e5cccf707ca75dafd993e5fd3 --- /dev/null +++ b/scripts/download_lofar_ms.sh @@ -0,0 +1,33 @@ +#!/bin/sh +# +# Author: Jakob Maljaars +# Email: jakob.maljaars_@_stcorp.nl + +# Script for downloading a mock LOFAR Measurement set + +set -e + +SCRIPT_PATH=$(dirname "$0") +cd $SCRIPT_PATH + +# Move up to parent folder which contains the source +cd .. +mkdir -p test_data +cd test_data/ + +LOFAR_MOCK_ARCHIVE=LOFAR_ARCHIVE.tar.bz2 +LOFAR_MOCK_MS=LOFAR_MOCK.ms + +if [ ! -f "$LOFAR_MOCK_ARCHIVE" ]; then + wget https://www.astron.nl/citt/EveryBeam/L258627-one-timestep.tar.bz2 -O $LOFAR_MOCK_ARCHIVE +fi + +if [ -d $LOFAR_MOCK_MS ] +then + echo "Directory already exists" +else + mkdir $LOFAR_MOCK_MS +fi + +tar -xf $LOFAR_MOCK_ARCHIVE -C $LOFAR_MOCK_MS --strip-components=1 +rm $LOFAR_MOCK_ARCHIVE \ No newline at end of file diff --git a/scripts/run-clang-format.sh b/scripts/run-clang-format.sh index fd7998fd6cb413012a2a736fbe84ef670af6599a..4191ca4ca414d11e978337aba83a6520beaee78f 100755 --- a/scripts/run-clang-format.sh +++ b/scripts/run-clang-format.sh @@ -1,4 +1,3 @@ - #!/bin/sh # # Author: Jakob Maljaars