diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index bcc721e533894259cb89b211d209e063b84d33c8..78e2f786bea6b334bcd804092d11599d9b665258 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -18,7 +18,7 @@ clang-format:
- apt-get -y install python3-pip
- pip3 install clang-format
script:
- - ./scripts/clang-format-check.sh
+ - ./scripts/run-clang-format.sh
# Build a debug version of EveryBeam from the base image
test-and-coverage:
@@ -30,16 +30,24 @@ test-and-coverage:
- apt-get -y install python3-pip
- pip3 install gcovr
script:
- - mkdir build
- # Download LOFAR/VLA Mock measurement set
- - ./scripts/download_lofar_ms.sh && cd test_data/LOFAR_MOCK.ms/ && export LOFAR_MOCK_MS=$PWD && cd ../../
- - ./scripts/download_vla_ms.sh && cd test_data/VLA_MOCK.ms/ && export VLA_MOCK_MS=$PWD && cd ../../build
+ # Download casacore wsrt measures
+ - wget -q ftp://ftp.astron.nl/outgoing/Measures/WSRT_Measures.ztar && tar -xf WSRT_Measures.ztar -C /var/lib/casacore/data/ && rm -f WSRT_Measures.ztar
+ # Download LOFAR/VLA/MWA Mock measurement set
+ - ./scripts/download_lofar_ms.sh && export LOFAR_MOCK_MS=$PWD/test_data/LOFAR_MOCK.ms/
+ - ./scripts/download_vla_ms.sh && export VLA_MOCK_MS=$PWD/test_data/VLA_MOCK.ms/
+ - ./scripts/download_mwa_ms.sh && export MWA_MOCK_MS=$PWD/test_data/MWA_MOCK.ms/
+ - ./scripts/download_mwa_coeff.sh && export MWA_COEFF_PATH=$PWD/test_data/mwa_full_embedded_element_pattern.h5
# Build in Debug mode and add LOFAR_MOCK_MS
- - cmake -DCMAKE_INSTALL_PREFIX=.. -DLOFAR_MOCK_MS=$LOFAR_MOCK_MS -DVLA_MOCK_MS=$VLA_MOCK_MS -DCMAKE_BUILD_TYPE=Debug -DCMAKE_CXX_FLAGS="-coverage" -DCMAKE_EXE_LINKER_FLAGS="-coverage" ..
+ - mkdir build && cd build
+ - |
+ cmake \
+ -DCMAKE_INSTALL_PREFIX=.. -DLOFAR_MOCK_MS=$LOFAR_MOCK_MS \
+ -DVLA_MOCK_MS=$VLA_MOCK_MS \
+ -DMWA_MOCK_MS=$MWA_MOCK_MS -DMWA_COEFF_PATH=$MWA_COEFF_PATH \
+ -DCMAKE_BUILD_TYPE=Debug -DCMAKE_CXX_FLAGS="-coverage" -DCMAKE_EXE_LINKER_FLAGS="-coverage" ..
- make install -j8
- ctest -T test
- # Explicit test that vla is checked
- - ./cpp/test/runtests --log_level=test_suite --run_test=load_vla
+ # Capture coverage
- gcovr -r .. -e '.*/external/.*' -e '.*/CompilerIdCXX/.*' -e '.*/test/.*' -e '.*/demo/.*'
- gcovr -r .. -e '.*/external/.*' -e '.*/CompilerIdCXX/.*' -e '.*/test/.*' -e '.*/demo/.*' --xml > coverage.xml
artifacts:
diff --git a/CMake/config.h.in b/CMake/config.h.in
index 29fe68a0093c1fd592cba6628021bd23fe03b328..2599abdbdd646366501ef8f809c4434627e59edf 100644
--- a/CMake/config.h.in
+++ b/CMake/config.h.in
@@ -5,5 +5,7 @@
#define TEST_MEASUREMENTSET "@TEST_MEASUREMENTSET@"
#define LOFAR_MOCK_MS "@LOFAR_MOCK_MS@"
#define VLA_MOCK_MS "@VLA_MOCK_MS@"
+#define MWA_MOCK_MS "@MWA_MOCK_MS@"
+#define MWA_COEFF_PATH "@MWA_COEFF_PATH@"
#endif
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d7819beab3902578098b8cbeb83cf86c10a445c3..551ba66a7b72dada77c568ffc794eb683c06b226 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -50,6 +50,10 @@ include_directories(${CASACORE_INCLUDE_DIR})
# Find and include OpenMP
find_package(OpenMP REQUIRED)
+# Find and include Boost headers (boost::math required for MWA beam)
+find_package(Boost REQUIRED)
+include_directories(${Boost_INCLUDE_DIR})
+
#------------------------------------------------------------------------------
# Set CMake and compiler options
if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME)
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index c1746ca68264a01002eb35586498c5c3adc1057f..9ce65a85d13b659028a540c95d386c63478731e2 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -21,10 +21,15 @@ add_library(everybeam SHARED
station.cc
telescope/lofar.cc
telescope/dish.cc
+ telescope/mwa.cc
griddedresponse/lofargrid.cc
griddedresponse/dishgrid.cc
+ griddedresponse/mwagrid.cc
circularsymmetric/voltagepattern.cc
circularsymmetric/vlabeam.cc
+ # mwabeam related
+ mwabeam/tilebeam2016.cc
+ mwabeam/beam2016implementation.cc
)
target_include_directories(everybeam PUBLIC ${CASACORE_INCLUDE_DIR})
diff --git a/cpp/circularsymmetric/vlabeam.cc b/cpp/circularsymmetric/vlabeam.cc
index d8bd3d3cec846d2fe189f157087b11ad7a006504..27de521560cd5c34e8f5451950cf41c74f119447 100644
--- a/cpp/circularsymmetric/vlabeam.cc
+++ b/cpp/circularsymmetric/vlabeam.cc
@@ -3,7 +3,7 @@
#include <cmath>
#include <vector>
-using namespace everybeam::circularsymmetric;
+using everybeam::circularsymmetric::VLABeam;
std::array<double, 5> VLABeam::GetCoefficients(const std::string& bandName,
double freq) {
diff --git a/cpp/circularsymmetric/voltagepattern.cc b/cpp/circularsymmetric/voltagepattern.cc
index bb7874cd36b97848633a2afb6e22a958bdbc2ad6..1983c3c743a04e687afd11627de86c9bbfdb14bb 100644
--- a/cpp/circularsymmetric/voltagepattern.cc
+++ b/cpp/circularsymmetric/voltagepattern.cc
@@ -1,15 +1,12 @@
#include "voltagepattern.h"
-// #include "../wsclean/logger.h"
-// #include "../wsclean/primarybeamimageset.h"
-
#include <aocommon/imagecoordinates.h>
#include <cmath>
-using namespace everybeam::circularsymmetric;
using aocommon::ImageCoordinates;
using aocommon::UVector;
+using everybeam::circularsymmetric::VoltagePattern;
void VoltagePattern::EvaluatePolynomial(const UVector<double>& coefficients,
bool invert) {
diff --git a/cpp/common/system.h b/cpp/common/system.h
deleted file mode 100644
index 8742f1c06d4c15a9319766d0bfa07c64a6600250..0000000000000000000000000000000000000000
--- a/cpp/common/system.h
+++ /dev/null
@@ -1,33 +0,0 @@
-#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/itrfconverter.cc b/cpp/coords/itrfconverter.cc
index afe533ed15e833fb874ad0cb768e585e462e2bb1..79f956943986df1b1ca8f60ab34b2f186ec5b348 100644
--- a/cpp/coords/itrfconverter.cc
+++ b/cpp/coords/itrfconverter.cc
@@ -7,9 +7,8 @@
#include <casacore/measures/Measures/MDirection.h>
#include <casacore/measures/Measures/MEpoch.h>
-using namespace everybeam;
-using namespace everybeam::coords;
-
+namespace everybeam {
+namespace coords {
// TODO: Initialize converter with a time (and fixed position) and convert
// specific directions.
// Needed for wslean as well as for the makeeverybeam executable.
@@ -79,3 +78,5 @@ casacore::MDirection ITRFConverter::ToDirection(
const casacore::MDirection &direction) const {
return converter_(direction);
}
+} // namespace coords
+} // namespace everybeam
\ No newline at end of file
diff --git a/cpp/coords/itrfdirection.cc b/cpp/coords/itrfdirection.cc
index c2dbd485bbfb5b2022d99d79abbadbfc446b6518..fb996d6f9d48bf84760b7724da4aa1c2a6cdf4fb 100644
--- a/cpp/coords/itrfdirection.cc
+++ b/cpp/coords/itrfdirection.cc
@@ -26,8 +26,8 @@
#include <casacore/measures/Measures/MDirection.h>
#include <casacore/measures/Measures/MEpoch.h>
-using namespace everybeam;
-using namespace everybeam::coords;
+namespace everybeam {
+namespace coords {
// ITRF position of CS002LBA, just to use a fixed reference
const vector3r_t ITRFDirection::lofar_position_ = {
@@ -90,3 +90,5 @@ vector3r_t ITRFDirection::at(real_t time) const {
vector3r_t itrf = {{mvITRF(0), mvITRF(1), mvITRF(2)}};
return itrf;
}
+} // namespace coords
+} // namespace everybeam
\ No newline at end of file
diff --git a/cpp/griddedresponse/CMakeLists.txt b/cpp/griddedresponse/CMakeLists.txt
index 300981f2e38b7342f5b910fd834b843d95c4f028..2ed7a567bb480cccd7bbd82d10f99053089d5e3e 100644
--- a/cpp/griddedresponse/CMakeLists.txt
+++ b/cpp/griddedresponse/CMakeLists.txt
@@ -2,4 +2,5 @@ install (FILES
griddedresponse.h
lofargrid.h
dishgrid.h
+ mwagrid.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/griddedresponse")
diff --git a/cpp/griddedresponse/dishgrid.cc b/cpp/griddedresponse/dishgrid.cc
index d340730d46ef2669c993048468acd8bb5cce3c27..c5bd5a74ff2deb964c0ab67e6e57ebf9e54e114f 100644
--- a/cpp/griddedresponse/dishgrid.cc
+++ b/cpp/griddedresponse/dishgrid.cc
@@ -7,12 +7,10 @@
#include <algorithm>
-using namespace everybeam;
-using namespace everybeam::griddedresponse;
+using everybeam::griddedresponse::DishGrid;
void DishGrid::CalculateStation(std::complex<float>* buffer, double,
- double frequency, const size_t,
- const size_t field_id) {
+ double frequency, size_t, size_t field_id) {
const telescope::Dish& dishtelescope =
static_cast<const telescope::Dish&>(*telescope_);
@@ -29,7 +27,7 @@ void DishGrid::CalculateStation(std::complex<float>* buffer, double,
};
void DishGrid::CalculateAllStations(std::complex<float>* buffer, double,
- double frequency, const size_t field_id) {
+ double frequency, size_t field_id) {
CalculateStation(buffer, 0., frequency, 0, field_id);
// Just repeat nstations times
@@ -37,4 +35,4 @@ void DishGrid::CalculateAllStations(std::complex<float>* buffer, double,
std::copy_n(buffer, width_ * height_ * 4,
buffer + i * width_ * height_ * 4);
}
-}
+}
\ No newline at end of file
diff --git a/cpp/griddedresponse/dishgrid.h b/cpp/griddedresponse/dishgrid.h
index 8df4f45c87cacec07eae2d1a4e6baab28aa049cd..e0a9bf5a681fbf774572c0824218b86c7140f259 100644
--- a/cpp/griddedresponse/dishgrid.h
+++ b/cpp/griddedresponse/dishgrid.h
@@ -1,4 +1,4 @@
-// vlagrid.h: Class for computing the VLA circular symmetric (gridded) response.
+// dishgrid.h: Class for computing the circular symmetric (gridded) response.
//
// Copyright (C) 2020
// ASTRON (Netherlands Institute for Radio Astronomy)
@@ -36,11 +36,11 @@ class DishGrid final : public GriddedResponse {
: GriddedResponse(telescope_ptr, coordinate_system){};
void CalculateStation(std::complex<float>* buffer, double time,
- double frequency, const size_t field_id,
- const size_t station_idx) override;
+ double frequency, size_t station_idx,
+ size_t field_id) override;
void CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency, const size_t field_id) override;
+ double frequency, size_t field_id) override;
};
} // namespace griddedresponse
} // namespace everybeam
diff --git a/cpp/griddedresponse/griddedresponse.h b/cpp/griddedresponse/griddedresponse.h
index ca158f6078e7ac94691aa7fa6ad2975fdbe7cb42..ce070ec3bb57ffdb748817ca265a957fade34ed9 100644
--- a/cpp/griddedresponse/griddedresponse.h
+++ b/cpp/griddedresponse/griddedresponse.h
@@ -59,8 +59,8 @@ class GriddedResponse {
* @param frequency Frequency (Hz)
*/
virtual void CalculateStation(std::complex<float>* buffer, double time,
- double freq, const size_t station_idx,
- const size_t field_id) = 0;
+ double freq, size_t station_idx,
+ size_t field_id) = 0;
/**
* @brief Compute the Beam response for all stations in a Telescope
@@ -71,10 +71,9 @@ class GriddedResponse {
* @param frequency Frequency (Hz)
*/
virtual void CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency,
- const size_t field_id) = 0;
+ double frequency, size_t field_id) = 0;
- std::size_t GetBufferSize(std::size_t nstations) {
+ std::size_t GetBufferSize(std::size_t nstations) const {
return std::size_t(nstations * width_ * height_ * 2 * 2);
}
diff --git a/cpp/griddedresponse/lofargrid.cc b/cpp/griddedresponse/lofargrid.cc
index 03abdc5c25cc4ce9d02fabd56b5481d4c824656c..51d3700614506f4578ccf8c6e9ad79fff0284ea0 100644
--- a/cpp/griddedresponse/lofargrid.cc
+++ b/cpp/griddedresponse/lofargrid.cc
@@ -1,19 +1,19 @@
#include "lofargrid.h"
#include "./../telescope/lofar.h"
-#include "./../common/system.h"
#include <aocommon/imagecoordinates.h>
+#include <aocommon/threadpool.h>
#include <cmath>
#include <iostream>
-using namespace everybeam::griddedresponse;
+using everybeam::griddedresponse::LOFARGrid;
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();
+ size_t ncpus = aocommon::ThreadPool::NCPUs();
// Set private members
nthreads_ = std::min(ncpus, lofartelescope.nstations_);
@@ -26,8 +26,7 @@ LOFARGrid::LOFARGrid(telescope::Telescope* telescope_ptr,
};
void LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
- double frequency, const size_t station_idx,
- const size_t) {
+ double frequency, size_t station_idx, size_t) {
const telescope::LOFAR& lofartelescope =
static_cast<const telescope::LOFAR&>(*telescope_);
aocommon::Lane<Job> lane(nthreads_);
@@ -63,7 +62,7 @@ void LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
}
void LOFARGrid::CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency, const size_t) {
+ double frequency, size_t) {
const telescope::LOFAR& lofartelescope =
static_cast<const telescope::LOFAR&>(*telescope_);
aocommon::Lane<Job> lane(nthreads_);
diff --git a/cpp/griddedresponse/lofargrid.h b/cpp/griddedresponse/lofargrid.h
index 50f9a515b3ba1e65bfe973f9bb623ff352f49cd7..a032993475f49d3b27b2c97c26a2814d64cb01e5 100644
--- a/cpp/griddedresponse/lofargrid.h
+++ b/cpp/griddedresponse/lofargrid.h
@@ -61,8 +61,8 @@ class LOFARGrid final : public GriddedResponse {
* @param freq Frequency (Hz)
*/
void CalculateStation(std::complex<float>* buffer, double time,
- double frequency, const size_t station_idx,
- const size_t field_id) override;
+ double frequency, size_t station_idx,
+ size_t field_id) override;
/**
* @brief Compute the Beam response for all stations in a Telescope
@@ -73,7 +73,7 @@ class LOFARGrid final : public GriddedResponse {
* @param freq Frequency (Hz)
*/
void CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency, const size_t field_id) override;
+ double frequency, size_t field_id) override;
private:
casacore::MDirection delay_dir_, tile_beam_dir_, preapplied_beam_dir_;
diff --git a/cpp/griddedresponse/mwagrid.cc b/cpp/griddedresponse/mwagrid.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d7c3cbec8857613858cd2d58e8c8b21017792c25
--- /dev/null
+++ b/cpp/griddedresponse/mwagrid.cc
@@ -0,0 +1,69 @@
+#include "mwagrid.h"
+#include "../telescope/mwa.h"
+
+#include <aocommon/imagecoordinates.h>
+
+#include <casacore/measures/Measures/MEpoch.h>
+#include <casacore/measures/Measures/MEpoch.h>
+#include <casacore/measures/Measures/MDirection.h>
+
+#include <memory>
+
+using everybeam::griddedresponse::MWAGrid;
+using everybeam::mwabeam::TileBeam2016;
+
+void MWAGrid::CalculateStation(std::complex<float>* buffer, double time,
+ double frequency, size_t station_idx, size_t) {
+ const telescope::MWA& mwatelescope =
+ static_cast<const telescope::MWA&>(*telescope_);
+ casacore::MEpoch time_epoch(casacore::Quantity(time, "s"));
+ casacore::MeasFrame frame(mwatelescope.ms_properties_.array_position,
+ time_epoch);
+
+ const casacore::MDirection::Ref hadec_ref(casacore::MDirection::HADEC, frame);
+ const casacore::MDirection::Ref azelgeo_ref(casacore::MDirection::AZELGEO,
+ frame);
+ const casacore::MDirection::Ref j2000_ref(casacore::MDirection::J2000, frame);
+ casacore::MDirection::Convert j2000_to_hadecref(j2000_ref, hadec_ref),
+ j2000_to_azelgeoref(j2000_ref, azelgeo_ref);
+ casacore::MPosition wgs = casacore::MPosition::Convert(
+ mwatelescope.ms_properties_.array_position, casacore::MPosition::WGS84)();
+ double arrLatitude = wgs.getValue().getLat();
+
+ if (!tile_beam_) {
+ tile_beam_.reset(
+ new TileBeam2016(mwatelescope.ms_properties_.delays,
+ mwatelescope.GetOptions().frequency_interpolation,
+ mwatelescope.GetOptions().coeff_path));
+ }
+ std::complex<float>* buffer_ptr = buffer;
+ for (size_t y = 0; y != height_; ++y) {
+ for (size_t x = 0; x != width_; ++x) {
+ double l, m, ra, dec;
+ aocommon::ImageCoordinates::XYToLM(x, y, dl_, dm_, width_, height_, l, m);
+ l += phase_centre_dl_;
+ m += phase_centre_dm_;
+ aocommon::ImageCoordinates::LMToRaDec(l, m, ra_, dec_, ra, dec);
+
+ std::complex<double> gain[4];
+ tile_beam_->ArrayResponse(ra, dec, j2000_ref, j2000_to_hadecref,
+ j2000_to_azelgeoref, arrLatitude, frequency,
+ gain);
+
+ for (size_t i = 0; i != 4; ++i) {
+ *buffer_ptr = gain[i];
+ ++buffer_ptr;
+ }
+ }
+ }
+};
+
+void MWAGrid::CalculateAllStations(std::complex<float>* buffer, double time,
+ double frequency, size_t) {
+ CalculateStation(buffer, time, frequency, 0, 0);
+ // Repeated copy for nstations
+ for (size_t i = 1; i != telescope_->GetNrStations(); ++i) {
+ std::copy_n(buffer, width_ * height_ * 4,
+ buffer + i * width_ * height_ * 4);
+ }
+};
\ No newline at end of file
diff --git a/cpp/griddedresponse/mwagrid.h b/cpp/griddedresponse/mwagrid.h
new file mode 100644
index 0000000000000000000000000000000000000000..601f16de962350e24da721e55041fbbd770be173
--- /dev/null
+++ b/cpp/griddedresponse/mwagrid.h
@@ -0,0 +1,52 @@
+// mwagrid.h: Class for computing the circular symmetric (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_MWAGRID_H_
+#define EVERYBEAM_GRIDDEDRESPONSE_MWAGRID_H_
+
+#include "griddedresponse.h"
+#include "../mwabeam/tilebeam2016.h"
+
+#include <memory>
+
+namespace everybeam {
+namespace griddedresponse {
+class MWAGrid final : public GriddedResponse {
+ public:
+ MWAGrid(telescope::Telescope* telescope_ptr,
+ const coords::CoordinateSystem coordinate_system)
+ : GriddedResponse(telescope_ptr, coordinate_system){};
+
+ void CalculateStation(std::complex<float>* buffer, double time,
+ double frequency, size_t station_idx,
+ size_t field_id) override;
+
+ void CalculateAllStations(std::complex<float>* buffer, double time,
+ double frequency, const size_t field_id) override;
+
+ private:
+ std::unique_ptr<everybeam::mwabeam::TileBeam2016> tile_beam_;
+};
+} // namespace griddedresponse
+} // namespace everybeam
+#endif // EVERYBEAM_GRIDDEDRESPONSE_MWAGRID_H_
diff --git a/cpp/load.cc b/cpp/load.cc
index a91810d40063bf90370396652bc4e5853bff9a06..93e8a33c27e36e56a4fa42ab762c0185745ecb06 100644
--- a/cpp/load.cc
+++ b/cpp/load.cc
@@ -13,7 +13,8 @@ enum TelescopeType {
kUnknownTelescope,
kLofarTelescope,
kVLATelescope,
- kATCATelescope
+ kATCATelescope,
+ kMWATelescope
};
TelescopeType Convert(const std::string &telescope_name) {
@@ -23,6 +24,8 @@ TelescopeType Convert(const std::string &telescope_name) {
return kVLATelescope;
else if (telescope_name.compare(0, 4, "ATCA") == 0)
return kATCATelescope;
+ else if (telescope_name == "MWA")
+ return kMWATelescope;
else
return kUnknownTelescope;
}
@@ -55,6 +58,12 @@ std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
new telescope::Dish(ms, options));
return telescope;
}
+ case kMWATelescope: {
+ std::unique_ptr<telescope::Telescope> telescope =
+ std::unique_ptr<telescope::Telescope>(
+ new telescope::MWA(ms, options));
+ return telescope;
+ }
default:
std::stringstream message;
message << "The requested telescope type " << telescope_name_col(0)
diff --git a/cpp/load.h b/cpp/load.h
index ff9ed40467e3d51eb07de2c766abc4f18777b8d9..50748b2430ee6e58d160db4e75eee3791620114a 100644
--- a/cpp/load.h
+++ b/cpp/load.h
@@ -27,6 +27,7 @@
#include "./telescope/telescope.h"
#include "./telescope/lofar.h"
#include "./telescope/dish.h"
+#include "./telescope/mwa.h"
#include "options.h"
namespace everybeam {
diff --git a/cpp/mwabeam/beam2016implementation.cc b/cpp/mwabeam/beam2016implementation.cc
new file mode 100644
index 0000000000000000000000000000000000000000..34f4e82de961cb6c4a0556e68781a2e8bbe35004
--- /dev/null
+++ b/cpp/mwabeam/beam2016implementation.cc
@@ -0,0 +1,731 @@
+/**
+ * C++ implementation of Full Embeded Element beam model for MWA based on
+ * beam_full_EE.py script and Sokolowski et al (2016) paper Implemented by
+ * Marcin Sokolowski (May 2017) - marcin.sokolowski@curtin.edu.au 20 May 2017 :
+ * Somewhat optimized by André Offringa.
+ */
+
+#include <algorithm>
+#include <cmath>
+#include <complex>
+#include <cstdlib>
+#include <cstring>
+#include <iostream>
+#include <string>
+#include <vector>
+
+#include <boost/math/special_functions/legendre.hpp>
+
+#include <H5Cpp.h>
+
+#include "beam2016implementation.h"
+
+using namespace std;
+using namespace H5;
+using namespace everybeam::mwabeam;
+
+// constants :
+static const double deg2rad = M_PI / 180.00;
+
+// beamformer step in pico-seconds
+#define DELAY_STEP 435.0e-12
+#define N_ANT_COUNT 16
+
+bool Beam2016Implementation::has_freq(int freq_hz) const {
+ return std::find(m_freq_list.begin(), m_freq_list.end(), freq_hz) !=
+ m_freq_list.end();
+}
+
+int Beam2016Implementation::find_closest_freq(int freq_hz) const {
+ double min_diff = 1e20;
+ int best_idx = -1;
+
+ for (size_t i = 0; i < m_freq_list.size(); i++) {
+ double diff = abs(m_freq_list[i] - freq_hz);
+ if (diff < min_diff) {
+ min_diff = diff;
+ best_idx = i;
+ }
+ }
+
+ if (best_idx >= 0) {
+ return m_freq_list[best_idx];
+ }
+
+ return m_freq_list[0];
+}
+
+Beam2016Implementation::Beam2016Implementation(const double* delays,
+ const double* amps,
+ const std::string& searchPath)
+ : _calcModesLastFreqHz(-1),
+ _calcModesLastDelays(),
+ _calcModesLastAmps(),
+ _h5File(),
+ _searchPath(searchPath),
+ _factorial(100) {
+ if (delays == nullptr)
+ std::fill(_delays, _delays + N_ANT_COUNT, 0.0);
+ else
+ std::copy_n(delays, N_ANT_COUNT, _delays);
+
+ if (amps == nullptr)
+ std::fill(_amps, _amps + N_ANT_COUNT, 1.0);
+ else
+ std::copy_n(amps, N_ANT_COUNT, _amps);
+
+ Read();
+ _jPowerTable[0] = std::complex<double>(1.0, 0.0);
+ _jPowerTable[1] = std::complex<double>(0, 1.0);
+ _jPowerTable[2] = std::complex<double>(-1.0, 0.0);
+ _jPowerTable[3] = std::complex<double>(0.0, -1.0);
+}
+
+//--------------------------------------------------------------------
+// Calculation of Jones matrix
+//------------------------------------------------------------------------------------------------------------------
+void Beam2016Implementation::CalcJonesArray(vector<vector<double>>& azim_arr,
+ vector<vector<double>>& za_arr,
+ vector<vector<JonesMatrix>>& jones,
+ int freq_hz_param,
+ bool bZenithNorm) {
+ // convert AZIM -> FEKO PHI phi=90-azim or azim=90-phi :
+ // python : phi_arr=math.pi/2-phi_arr #Convert to East through North (FEKO
+ // coords)
+ for (size_t y = 0; y < azim_arr.size(); y++) {
+ vector<double>& image_row = azim_arr[y];
+
+ for (size_t x = 0; x < image_row.size(); x++) {
+ image_row[x] = M_PI / 2.00 - image_row[x];
+
+ // phi_arr[phi_arr < 0] += 2*math.pi #360 wrap
+ if (image_row[x] < 0) {
+ image_row[x] += 2.0 * M_PI;
+ }
+ }
+ }
+
+ JonesMatrix::zeros(jones, azim_arr[0].size(), azim_arr.size());
+
+ for (size_t y = 0; y < azim_arr.size(); y++) {
+ vector<double>& image_row = azim_arr[y];
+ for (size_t x = 0; x < image_row.size(); x++) {
+ double azim_deg = image_row[x];
+ double za_deg = (za_arr[y])[x];
+
+ jones[y][x] = CalcJones(azim_deg, za_deg, freq_hz_param, bZenithNorm);
+ }
+ }
+}
+
+JonesMatrix Beam2016Implementation::CalcJonesDirect(
+ double az_rad, double za_rad, const Coefficients& coefsX,
+ const Coefficients& coefsY) {
+ // convert AZIM -> FEKO PHI phi=90-azim or azim=90-phi :
+ // python : phi_arr=math.pi/2-phi_arr #Convert to East through North (FEKO
+ // coords)
+ JonesMatrix jones;
+ double phi_rad = M_PI / 2.00 - az_rad;
+
+ CalcSigmas(phi_rad, za_rad, coefsX, 'X', jones);
+ CalcSigmas(phi_rad, za_rad, coefsY, 'Y', jones);
+
+ return jones;
+}
+
+JonesMatrix Beam2016Implementation::CalcJones(double az_deg, double za_deg,
+ int freq_hz_param,
+ bool bZenithNorm) {
+ RecursiveLock<std::mutex> lock(mutex_, std::defer_lock);
+ return CalcJones(az_deg, za_deg, freq_hz_param, _delays, _amps, lock,
+ bZenithNorm);
+}
+
+JonesMatrix Beam2016Implementation::CalcJones(double az_deg, double za_deg,
+ int freq_hz, const double* delays,
+ const double* amps,
+ RecursiveLock<std::mutex>& lock,
+ bool bZenithNorm) {
+ if (!has_freq(freq_hz)) {
+ freq_hz = find_closest_freq(freq_hz);
+ }
+
+ Coefficients coefsX, coefsY;
+ GetModes(freq_hz, N_ANT_COUNT, delays, amps, coefsX, coefsY, lock);
+
+ JonesMatrix result =
+ CalcJonesDirect(az_deg * deg2rad, za_deg * deg2rad, coefsX, coefsY);
+
+ if (bZenithNorm) {
+ JonesMatrix normMatrix;
+
+ std::lock_guard<RecursiveLock<std::mutex>> glock(lock);
+ std::map<int, JonesMatrix>::const_iterator iter =
+ _normJonesCache.find(freq_hz);
+ if (iter == _normJonesCache.end()) {
+ normMatrix = CalcZenithNormMatrix(freq_hz, lock);
+ _normJonesCache.insert(std::make_pair(freq_hz, normMatrix));
+ } else {
+ normMatrix = iter->second;
+ }
+
+ result.j00 = result.j00 / normMatrix.j00;
+ result.j01 = result.j01 / normMatrix.j01;
+ result.j10 = result.j10 / normMatrix.j10;
+ result.j11 = result.j11 / normMatrix.j11;
+ }
+
+ return result;
+}
+
+JonesMatrix Beam2016Implementation::CalcZenithNormMatrix(
+ int freq_hz, RecursiveLock<std::mutex>& lock) {
+ // std::cout << "INFO : calculating Jones matrix for frequency = " << freq_hz
+ // << " Hz\n";
+
+ // Azimuth angles at which Jones components are maximum (see beam_full_EE.py
+ // for comments):
+ // max_phis=[[math.pi/2,math.pi],[0,math.pi/2]] #phi where each Jones vector
+ // is max
+ double j00_max_az = 90.00;
+ double j01_max_az = 180.00;
+ double j10_max_az = 0.00;
+ double j11_max_az = 90.00;
+
+ JonesMatrix tmp_jones = JonesMatrix(1, 1, 1, 1);
+ JonesMatrix jonesMatrix;
+
+ // default delays at zenith
+ const double defaultDelays[N_ANT_COUNT] = {0, 0, 0, 0, 0, 0, 0, 0,
+ 0, 0, 0, 0, 0, 0, 0, 0};
+ const double defaultAmps[N_ANT_COUNT] = {1, 1, 1, 1, 1, 1, 1, 1,
+ 1, 1, 1, 1, 1, 1, 1, 1};
+
+ // j00 :
+ tmp_jones = CalcJones(j00_max_az, 0, freq_hz, defaultDelays, defaultAmps,
+ lock, false);
+ jonesMatrix.j00 = abs(tmp_jones.j00);
+
+ // j01 :
+ tmp_jones = CalcJones(j01_max_az, 0, freq_hz, defaultDelays, defaultAmps,
+ lock, false);
+ jonesMatrix.j01 = abs(tmp_jones.j01);
+
+ // j10 :
+ tmp_jones = CalcJones(j10_max_az, 0, freq_hz, defaultDelays, defaultAmps,
+ lock, false);
+ jonesMatrix.j10 = abs(tmp_jones.j10);
+
+ // j11 :
+ tmp_jones = CalcJones(j11_max_az, 0, freq_hz, defaultDelays, defaultAmps,
+ lock, false);
+ jonesMatrix.j11 = abs(tmp_jones.j11);
+
+ return jonesMatrix;
+}
+
+//--------------------------------------------------------------------
+// Calculation of Jones matrix components for given polarisation (eq. 4 and 5 in
+// the Sokolowski et al (2017) paper --------------------------------
+void Beam2016Implementation::CalcSigmas(double phi, double theta,
+ const Coefficients& coefs, char pol,
+ JonesMatrix& jones_matrix) const {
+ int nmax = int(coefs.Nmax);
+
+ double u = cos(theta);
+
+ vector<double> P1sin_arr, P1_arr;
+
+ P1sin(nmax, theta, P1sin_arr, P1_arr);
+
+ if (coefs.N_accum.size() != coefs.M_accum.size()) {
+ throw std::runtime_error("ERROR : size of N_accnum != M_accum ( " +
+ std::to_string(coefs.N_accum.size()) + " != " +
+ std::to_string(coefs.M_accum.size()) + ")");
+ }
+
+ complex<double> sigma_P(0, 0), sigma_T(0, 0);
+ for (size_t i = 0; i < coefs.N_accum.size(); i++) {
+ double N = coefs.N_accum[i];
+ int n = int(N);
+
+ double M = coefs.M_accum[i];
+ // int m = int(M);
+ double m_abs_m = coefs.MabsM[i];
+
+ double c_mn_sqr =
+ (0.5 * (2 * N + 1) * _factorial(N - abs(M)) / _factorial(N + abs(M)));
+ double c_mn = sqrt(c_mn_sqr);
+
+ complex<double> ejm_phi(cos(M * phi), sin(M * phi));
+ complex<double> phi_comp = (ejm_phi * c_mn) / (sqrt(N * (N + 1))) * m_abs_m;
+
+ complex<double> j_power_n = JPower(n);
+ complex<double> E_theta_mn =
+ j_power_n * (P1sin_arr[i] * (fabs(M) * coefs.Q2_accum[i] * u -
+ M * coefs.Q1_accum[i]) +
+ coefs.Q2_accum[i] * P1_arr[i]);
+
+ complex<double> j_power_np1 = JPower(n + 1);
+ complex<double> E_phi_mn =
+ j_power_np1 * (P1sin_arr[i] * (M * coefs.Q2_accum[i] -
+ fabs(M) * coefs.Q1_accum[i] * u) -
+ coefs.Q1_accum[i] * P1_arr[i]);
+
+ sigma_P = sigma_P + phi_comp * E_phi_mn;
+ sigma_T = sigma_T + phi_comp * E_theta_mn;
+ }
+
+ if (pol == 'X') {
+ jones_matrix.j00 = sigma_T;
+ jones_matrix.j01 =
+ -sigma_P; // as it is now in python code in sign_fix branch
+ } else {
+ jones_matrix.j10 = sigma_T;
+ jones_matrix.j11 =
+ -sigma_P; // as it is now in python code in sign_fix branch
+ }
+}
+
+//--------------------------------------------------------------------
+// Calculation of spherical harmonics coefficients (eq. 3-6 in the Sokolowski et
+// al 2016 paper) --------------------------------
+// function comparing current parameters : frequency, delays and amplitudes with
+// those previously used to calculate spherical waves coefficients (stored in
+// the 3 variables : m_CalcModesLastFreqHz , , m_CalcModesLastAmps )
+bool Beam2016Implementation::IsCalcModesRequired(int freq_hz, int n_ant,
+ const double* delays,
+ const double* amps) {
+ if (freq_hz != _calcModesLastFreqHz) {
+ return true;
+ }
+
+ if (_calcModesLastDelays.empty() || _calcModesLastAmps.empty()) {
+ return true;
+ }
+
+ for (int i = 0; i < n_ant; i++) {
+ if (delays[i] != _calcModesLastDelays[i]) return true;
+ if (amps[i] != _calcModesLastAmps[i]) return true;
+ }
+
+ return false;
+}
+
+// function calculating coefficients for X and Y and storing parameters
+// frequency, delays and amplitudes
+void Beam2016Implementation::GetModes(int freq_hz, size_t n_ant,
+ const double* delays, const double* amps,
+ Coefficients& coefsX,
+ Coefficients& coefsY,
+ RecursiveLock<std::mutex>& lock) {
+ std::unique_lock<RecursiveLock<std::mutex>> glock(lock);
+ if (IsCalcModesRequired(freq_hz, n_ant, delays, amps)) {
+ glock.unlock();
+ coefsX.Nmax = CalcModes(freq_hz, n_ant, delays, amps, 'X', coefsX, lock);
+ coefsY.Nmax = CalcModes(freq_hz, n_ant, delays, amps, 'Y', coefsY, lock);
+
+ glock.lock();
+ _coefX = coefsX;
+ _coefY = coefsY;
+ _calcModesLastFreqHz = freq_hz;
+ _calcModesLastDelays.assign(delays, delays + n_ant);
+ _calcModesLastAmps.assign(amps, amps + n_ant);
+ } else {
+ coefsX = _coefX;
+ coefsY = _coefY;
+ }
+}
+
+// function calculating all coefficients Q1, Q2, N, M and derived MabsM, Nmax
+// for a given polarisation (X or Y)
+double Beam2016Implementation::CalcModes(int freq_hz, size_t n_ant,
+ const double* delays,
+ const double* amp, char pol,
+ Coefficients& coefs,
+ RecursiveLock<std::mutex>& lock) {
+ vector<double> phases(n_ant);
+ double Nmax = 0;
+ coefs.M_accum.clear();
+ coefs.N_accum.clear();
+ coefs.MabsM.clear();
+ coefs.Cmn.clear();
+
+ int modes_size = m_Modes[0].size();
+ coefs.Q1_accum.assign(modes_size, 0.0);
+ coefs.Q2_accum.assign(modes_size, 0.0);
+
+ for (size_t a = 0; a < n_ant; a++) {
+ double phase = 2 * M_PI * freq_hz * (-double(delays[a]) * DELAY_STEP);
+
+ phases[a] = phase;
+
+ // complex excitation voltage:
+ // self.amps[pol]*np.exp(1.0j*phases)
+ complex<double> phase_factor(cos(phase), sin(phase));
+ complex<double> Vcplx = amp[a] * phase_factor;
+
+ DataSetIndex index(pol, a, freq_hz);
+
+ const vector<vector<double>>& Q_all = GetDataSet(index, lock);
+
+ size_t n_ant_coeff = Q_all[0].size();
+ vector<double> Ms1, Ns1, Ms2, Ns2;
+ const vector<double>& modes_Type = m_Modes[0];
+ const vector<double>& modes_M = m_Modes[1];
+ const vector<double>& modes_N = m_Modes[2];
+
+ int bUpdateNAccum = 0;
+
+ // list of indexes where S=1 coefficients seat in array m_Modes ( and
+ // m_Modes_M and m_Modes_N )
+ vector<int> s1_list;
+ // idem for S=2
+ vector<int> s2_list;
+
+ for (size_t coeff = 0; coeff < n_ant_coeff; coeff++) {
+ int mode_type = modes_Type[coeff];
+
+ if (mode_type <= 1) {
+ // s=1 modes :
+ s1_list.push_back(coeff);
+
+ Ms1.push_back(modes_M[coeff]);
+ Ns1.push_back(modes_N[coeff]);
+ if (modes_N[coeff] > Nmax) {
+ Nmax = modes_N[coeff];
+ bUpdateNAccum = 1;
+ }
+ } else {
+ // s=2 modes :
+ s2_list.push_back(coeff);
+
+ Ms2.push_back(modes_M[coeff]);
+ Ns2.push_back(modes_N[coeff]);
+ }
+ }
+
+ if (bUpdateNAccum > 0) {
+ coefs.N_accum = Ns1;
+ coefs.M_accum = Ms1;
+ }
+
+ if (s1_list.size() != s2_list.size() ||
+ s2_list.size() != (n_ant_coeff / 2)) {
+ throw std::runtime_error(
+ "Wrong number of coefficients for s1 and s2 condition " +
+ std::to_string(s1_list.size()) +
+ " = =" + std::to_string(s2_list.size()) +
+ " == " + std::to_string(n_ant_coeff / 2) + " not satisfied");
+ }
+
+ vector<std::complex<double>> Q1, Q2;
+ const vector<double>& Q_all_0 = Q_all[0];
+ const vector<double>& Q_all_1 = Q_all[1];
+ int my_len_half = (n_ant_coeff / 2);
+
+ for (int i = 0; i < my_len_half; i++) {
+ // calculate Q1:
+ int s1_idx = s1_list[i];
+ double s10_coeff = Q_all_0[s1_idx];
+ double s11_coeff = Q_all_1[s1_idx];
+ double arg = s11_coeff * deg2rad;
+ complex<double> tmp(cos(arg), sin(arg));
+ complex<double> q1_val = s10_coeff * tmp;
+ Q1.push_back(q1_val);
+
+ // calculate Q2:
+ int s2_idx = s2_list[i];
+ double s20_coeff = Q_all_0[s2_idx];
+ double s21_coeff = Q_all_1[s2_idx];
+ double arg2 = s21_coeff * deg2rad;
+ complex<double> tmp2(cos(arg2), sin(arg2));
+ complex<double> q2_val = s20_coeff * tmp2;
+ Q2.push_back(q2_val);
+
+ coefs.Q1_accum[i] = coefs.Q1_accum[i] + q1_val * Vcplx;
+ coefs.Q2_accum[i] = coefs.Q2_accum[i] + q2_val * Vcplx;
+ }
+ }
+
+ // Same as tragic python code:
+ // MabsM=-M/np.abs(M)
+ // MabsM[MabsM==np.NaN]=1 #for M=0, replace NaN with MabsM=1;
+ // MabsM=(MabsM)**M
+ for (size_t i = 0; i < coefs.M_accum.size(); i++) {
+ int m = int(coefs.M_accum[i]);
+
+ double m_abs_m = 1;
+ if (m > 0) {
+ if ((m % 2) != 0) {
+ m_abs_m = -1;
+ }
+ }
+
+ coefs.MabsM.push_back(m_abs_m);
+ }
+
+ return Nmax;
+}
+
+//------------------------------------------------------------------------------------------------------
+// maths functions and wrappers
+//---------------------------------------------------------------------------------------
+// OUTPUT : returns list of Legendre polynomial values calculated up to order
+// nmax :
+int Beam2016Implementation::P1sin(int nmax, double theta,
+ vector<double>& p1sin_out,
+ vector<double>& p1_out) {
+ int size = nmax * nmax + 2 * nmax;
+ p1sin_out.resize(size);
+ p1_out.resize(size);
+
+ double sin_th = std::sin(theta);
+ double u = std::cos(theta);
+ double delu = 1e-6;
+
+ vector<double> P, Pm1, Pm_sin, Pu_mdelu, Pm_sin_merged, Pm1_merged;
+ P.reserve(nmax + 1);
+ Pm1.reserve(nmax + 1);
+ Pm_sin.reserve(nmax + 1);
+ Pu_mdelu.reserve(nmax + 1);
+ Pm_sin_merged.reserve(nmax * 2 + 1);
+ Pm1_merged.reserve(nmax * 2 + 1);
+
+ // Create a look-up table for the legendre polynomials
+ // Such that legendre_table[ m * nmax + (n-1) ] = legendre(n, m, u)
+ vector<double> legendre_table(nmax * (nmax + 1));
+ vector<double>::iterator legendre_iter = legendre_table.begin();
+ for (int m = 0; m != nmax + 1; ++m) {
+ double leg0 = boost::math::legendre_p(0, m, u);
+ double leg1 = boost::math::legendre_p(1, m, u);
+ *legendre_iter = leg1;
+ ++legendre_iter;
+ for (int n = 2; n != nmax + 1; ++n) {
+ if (m < n)
+ *legendre_iter = boost::math::legendre_next(n - 1, m, u, leg1, leg0);
+ else if (m == n)
+ *legendre_iter = boost::math::legendre_p(n, m, u);
+ else
+ *legendre_iter = 0.0;
+ leg0 = leg1;
+ leg1 = *legendre_iter;
+ ++legendre_iter;
+ }
+ }
+
+ for (int n = 1; n <= nmax; n++) {
+ P.resize(n + 1);
+ // Assign P[m] to legendre(n, m, u)
+ // This is equal to:
+ // lpmv(P, n , u);
+ for (size_t m = 0; m != size_t(n) + 1; ++m)
+ P[m] = legendre_table[m * nmax + (n - 1)];
+
+ // skip first 1 and build table Pm1 (P minus 1 )
+ Pm1.assign(P.begin() + 1, P.end());
+ Pm1.push_back(0);
+
+ Pm_sin.assign(n + 1, 0.0);
+ if (u == 1 || u == -1) {
+ // In this case we take the easy approach and don't use
+ // precalculated polynomials, since this path does not occur often.
+
+ // TODO This path doesn't make sense.
+ // I think Marcin should look at this; this path only occurs on polar
+ // positions so is rare, but what is done here looks wrong: first
+ // calculate *all* polynomials, then only use the pol for m=1. Pm_sin for
+ // indices 0 and >=2 are not calculated, this seems not right.
+ Pu_mdelu.resize(1);
+ lpmv(Pu_mdelu, n, u - delu);
+
+ // Pm_sin[1,0]=-(P[0]-Pu_mdelu[0])/delu #backward difference
+ if (u == -1)
+ Pm_sin[1] = -(Pu_mdelu[0] - P[0]) / delu; // #forward difference
+ else
+ Pm_sin[1] = -(P[0] - Pu_mdelu[0]) / delu;
+ } else {
+ for (size_t i = 0; i < P.size(); i++) {
+ Pm_sin[i] = P[i] / sin_th;
+ }
+ }
+
+ Pm_sin_merged.assign(Pm_sin.rbegin(), Pm_sin.rend() - 1);
+ Pm_sin_merged.insert(Pm_sin_merged.end(), Pm_sin.begin(), Pm_sin.end());
+
+ int ind_start =
+ (n - 1) * (n - 1) + 2 * (n - 1); // #start index to populate
+ int ind_stop = n * n + 2 * n; //#stop index to populate
+
+ // P_sin[np.arange(ind_start,ind_stop)]=np.append(np.flipud(Pm_sin[1::,0]),Pm_sin)
+ int modified = 0;
+ for (int i = ind_start; i < ind_stop; i++) {
+ p1sin_out[i] = (Pm_sin_merged[modified]);
+ modified++;
+ }
+
+ // P1[np.arange(ind_start,ind_stop)]=np.append(np.flipud(Pm1[1::,0]),Pm1)
+ Pm1_merged.assign(Pm1.rbegin(), Pm1.rend() - 1);
+ Pm1_merged.insert(Pm1_merged.end(), Pm1.begin(), Pm1.end());
+ modified = 0;
+ for (int i = ind_start; i < ind_stop; i++) {
+ p1_out[i] = Pm1_merged[modified];
+ modified++;
+ }
+ }
+
+ return nmax;
+}
+
+// Legendre polynomials :
+void Beam2016Implementation::lpmv(vector<double>& output, int n, double x) {
+ for (size_t order = 0; order < output.size(); order++) {
+ double val = boost::math::legendre_p<double>(n, order, x);
+ output[order] = val;
+ }
+}
+
+//-----------------------------------------------------------------------------------
+// HDF5 File interface and data structures for H5 data
+//-----------------------------------------------------------------------------------
+// This function goes thorugh all dataset names and records them info list of
+// strings : Beam2016Implementation::m_obj_list
+herr_t Beam2016Implementation::list_obj_iterate(hid_t /*loc_id*/,
+ const char* name,
+ const H5O_info_t* info,
+ void* operator_data) {
+ string szTmp;
+ Beam2016Implementation* pBeamModelPtr =
+ (Beam2016Implementation*)operator_data;
+
+ if (pBeamModelPtr == nullptr)
+ throw std::runtime_error(
+ "The pointer to Beam2016Implementation class in "
+ "Beam2016Implementation::list_obj_iterate must not be null");
+
+ if (name[0] == '.') { /* Root group, do not print '.' */
+ } else {
+ switch (info->type) {
+ case H5O_TYPE_GROUP:
+ case H5O_TYPE_NAMED_DATATYPE:
+ default:
+ break;
+ case H5O_TYPE_DATASET:
+ szTmp = name;
+ pBeamModelPtr->m_obj_list.push_back(szTmp);
+ break;
+ }
+ }
+
+ return 0;
+}
+
+const std::vector<std::vector<double>>& Beam2016Implementation::GetDataSet(
+ const DataSetIndex& index, RecursiveLock<std::mutex>& lock) {
+ std::lock_guard<RecursiveLock<std::mutex>> glock(lock);
+ auto iter = _dataSetCache.find(index);
+ if (iter == _dataSetCache.end()) {
+ iter =
+ _dataSetCache.emplace(index, std::vector<std::vector<double>>()).first;
+ ReadDataSet(index.Name(), iter->second, *_h5File);
+ }
+ return iter->second;
+}
+
+// Read dataset_name from H5 file
+void Beam2016Implementation::ReadDataSet(const std::string& dataset_name,
+ vector<vector<double>>& out_vector,
+ H5::H5File& h5File) {
+ DataSet modes = h5File.openDataSet(dataset_name.c_str());
+ DataSpace modes_dataspace = modes.getSpace();
+ int rank = modes_dataspace.getSimpleExtentNdims();
+ hsize_t dims_out[2];
+ modes_dataspace.getSimpleExtentDims(dims_out, NULL);
+ modes_dataspace.selectAll();
+
+ aocommon::UVector<float> data(dims_out[0] * dims_out[1]);
+ aocommon::UVector<float*> modes_data(dims_out[0]);
+ for (size_t i = 0; i < dims_out[0]; i++)
+ modes_data[i] = &data[i * dims_out[1]];
+
+ DataSpace memspace(rank, dims_out);
+ modes.read(data.data(), PredType::NATIVE_FLOAT, memspace, modes_dataspace);
+
+ for (size_t i = 0; i != dims_out[0]; i++) {
+ const float* startElement = &data[i * dims_out[1]];
+ const float* endElement = startElement + dims_out[1];
+ out_vector.emplace_back(startElement, endElement);
+ }
+}
+
+// Interface to HDF5 file format and structures to store H5 data
+// Functions for reading H5 file and its datasets :
+// Read data from H5 file, file name is specified in the object constructor
+void Beam2016Implementation::Read() {
+ std::string h5_path = _searchPath;
+ _h5File.reset(new H5File(h5_path.c_str(), H5F_ACC_RDONLY));
+
+ // hid_t group_id = m_pH5File->getId();
+ hid_t file_id = _h5File->getId();
+
+ /* TODO : not sure how to read attribute with the official HDF5 library ...
+ if( H5Aexists( file_id, "VERSION" ) ){
+ char szVersion[128];
+ strcpy(szVersion,"TEST");
+ //H5std_string szVersion;
+
+ hid_t attr_id = H5Aopen(file_id, "VERSION", H5P_DEFAULT );
+ hid_t attr_type = H5Aget_type(attr_id);
+ herr_t err = H5Aread( attr_id, attr_type, (void*)szVersion );
+ printf("Ids = %d -> %d -> type =
+ %d\n",file_id,attr_id,attr_type); printf("Version of the %s file = %s (err =
+ %d)\n",m_h5file.c_str(),szVersion,err); }else{ printf("ERROR : attribute
+ version does not exist\n");
+ }*/
+
+ m_obj_list.clear();
+ m_freq_list.clear();
+ herr_t status =
+ H5Ovisit(file_id, H5_INDEX_NAME, H5_ITER_NATIVE, list_obj_iterate, this);
+ if (status < 0) {
+ throw std::runtime_error(
+ "H5Ovisit returned with negative value which indicates a critical "
+ "error");
+ }
+
+ int max_ant_idx = -1;
+ for (size_t i = 0; i < m_obj_list.size(); i++) {
+ const char* key = m_obj_list[i].c_str();
+ if (strstr(key, "X1_")) {
+ const char* szFreq = key + 3;
+ m_freq_list.push_back(atol(szFreq));
+ }
+
+ if (key[0] == 'X') {
+ int ant_idx = 0, freq_hz = 0;
+ int scanf_ret = sscanf(key, "X%d_%d", &ant_idx, &freq_hz);
+ if (scanf_ret == 2) {
+ if (ant_idx > max_ant_idx) {
+ max_ant_idx = ant_idx;
+ }
+ }
+ }
+ }
+ // number of antenna is read from the file
+ if (max_ant_idx != N_ANT_COUNT) {
+ throw std::runtime_error(
+ "Number of simulated antennae = " + std::to_string(max_ant_idx) +
+ ", the code is currently implemented for " +
+ std::to_string(N_ANT_COUNT));
+ }
+
+ std::sort(m_freq_list.begin(), m_freq_list.end());
+
+ ReadDataSet("modes", m_Modes, *_h5File);
+}
+
+void JonesMatrix::zeros(vector<vector<JonesMatrix>>& jones, size_t x_size,
+ size_t y_size) {
+ vector<JonesMatrix> zero_vector(x_size, JonesMatrix(0, 0, 0, 0));
+ jones.assign(y_size, zero_vector);
+}
diff --git a/cpp/mwabeam/beam2016implementation.h b/cpp/mwabeam/beam2016implementation.h
new file mode 100644
index 0000000000000000000000000000000000000000..63cb233aacbbed17bb7abacce894465524cde026
--- /dev/null
+++ b/cpp/mwabeam/beam2016implementation.h
@@ -0,0 +1,248 @@
+#ifndef _GET_FF_H__
+#define _GET_FF_H__
+
+/**
+ * C++ implementation of Full Embeded Element beam model for MWA based on
+ * beam_full_EE.py script and Sokolowski et al (2017) paper Implemented by
+ * Marcin Sokolowski (May 2017) - marcin.sokolowski@curtin.edu.au
+ */
+
+#include <algorithm>
+#include <complex>
+#include <cstdlib>
+#include <cstring>
+#include <string>
+#include <vector>
+#include <map>
+#include <memory>
+#include <mutex>
+
+#include <boost/math/special_functions/legendre.hpp>
+
+#include <H5Cpp.h>
+
+#include "factorialtable.h"
+#include "recursivelock.h"
+
+namespace everybeam::mwabeam {
+
+/*
+ Structure for Jones matrix :
+ j00 j01
+ j10 j11
+*/
+class JonesMatrix {
+ public:
+ std::complex<double> j00, j01, j10, j11;
+
+ JonesMatrix(double j00_r = 0.00, double j01_r = 0.00, double j10_r = 0.00,
+ double j11_r = 0.00)
+ : j00(j00_r, 0), j01(j01_r, 0), j10(j10_r, 0), j11(j11_r, 0) {}
+
+ static void zeros(std::vector<std::vector<JonesMatrix>>& jones, size_t x_size,
+ size_t y_size);
+};
+
+class Beam2016Implementation {
+ public:
+ struct DataSetIndex {
+ DataSetIndex(char pol_, int antenna_, int freqHz_)
+ : pol(pol_), antenna(antenna_), frequencyHz(freqHz_) {}
+
+ bool operator<(const DataSetIndex& rhs) const {
+ return std::make_tuple(pol, antenna, frequencyHz) <
+ std::make_tuple(rhs.pol, rhs.antenna, rhs.frequencyHz);
+ }
+
+ std::string Name() const {
+ return pol + std::to_string(antenna + 1) + '_' +
+ std::to_string(frequencyHz);
+ }
+
+ char pol;
+ int antenna, frequencyHz;
+ };
+
+ Beam2016Implementation(const double* delays, const double* amps,
+ const std::string& searchPath);
+
+ //--------------------------------------------------------------------
+ // Calculation of Jones matrix
+ //------------------------------------------------------------------------------------------------------------------
+ // Calculate jones matrix in a specified direction for a given frequency,
+ // delays and amplitudes. Zenith normalisation can also be disabled - but by
+ // default is enabled : This function will :
+ // - calculate coefficients of spherical waves (SPW) if required (meaning
+ // if frequency, delays or amplitudes are different than last set of
+ // coefficients was calculated for)
+ // - calculate electric fields (equations 4 and 5 in Sokolowski et al
+ // paper)
+ // - normalise Jones matrix to one at zenith at the same frequency (if
+ // required by parameter):
+ //
+ // Only thse the functions should be used for external calls. The rest is
+ // protected from external use.
+ //
+ // INPUT : (az_deg,za_deg) - azimuth and zenith angles [in degrees]
+ // freq_hz_param - frequency in Hz
+ // delays - delays in beamformer steps
+ // amps - amplitudes
+ // bZenithNorm - normalise to zenith (>0) or not (<=0)
+ // OUTPUT : Jones matrix (normalised or not - depending on the parameter
+ // bZenithNorm )
+ JonesMatrix CalcJones(double az_deg, double za_deg, int freq_hz_param,
+ bool bZenithNorm = true);
+ JonesMatrix CalcJones(double az_deg, double za_deg, int freq_hz_param,
+ const double* delays, const double* amps,
+ RecursiveLock<std::mutex>& lock,
+ bool bZenithNorm = true);
+
+ // Calculation of Jones matrix for an image passed in the arrays azimuth and
+ // zenith angles maps. This function calls the single direction one (above)
+ // for all pixel in the input image. INPUT :
+ // 2D array of azimuth angles in degrees
+ // 2D array of zenith angles in degrees
+ // freq_hz_param - frequency in Hz
+ // delays - delays in beamformer steps
+ // amps - amplutudes
+ // bZenithNorm - normalise to zenith (>0) or not (<=0)
+ // OUTPUT :
+ // 2D array of JonesMatrix for each pixel
+ void CalcJonesArray(std::vector<std::vector<double>>& azim_arr,
+ std::vector<std::vector<double>>& za_arr,
+ std::vector<std::vector<JonesMatrix>>& jones,
+ int freq_hz_param, bool bZenithNorm = true);
+
+ private:
+ //---------------------------------------------------- Calculations and
+ // variables for spherical waves coefficients (see equations 3-6 in the
+ // Sokolowski et al paper)
+ // ----------------------------------------------------
+ // Coefficients of spherical waves (SPW) - see equations 3-6 in the Sokolowski
+ // et al paper 1 refers to s=1 (transverse electric - TE modes) - eq.5 2
+ // refers to s=2 (transverse magnetic - TM modes) - eq.5 These coefficients
+ // are calculated ones for a given frequency, delays, amplitudes so these
+ // parameters are stored to only calculate new coefficients when they change.
+ // X polarisation :
+ struct Coefficients {
+ std::vector<std::complex<double>> Q1_accum;
+ std::vector<std::complex<double>> Q2_accum;
+ std::vector<double> M_accum;
+ std::vector<double> N_accum;
+ std::vector<double>
+ MabsM; // precalculated m/abs(m) to make it once for all pointings
+ double Nmax; // maximum N coefficient for Y (=max(N_accum_X)) - to avoid
+ // relaculations
+ std::vector<double>
+ Cmn; // coefficient under sumation in equation 3 for X pol.
+ } _coefX, _coefY;
+
+ // Calculation of Jones matrix for a single pointing direction (internal
+ // function): INPUT : (az_rad,za_rad) - azimuth and zenith in [radians]
+ JonesMatrix CalcJonesDirect(double az_rad, double za_rad,
+ const Coefficients& coefsX,
+ const Coefficients& coefsY);
+
+ //--------------------------------------------------------------------
+ // Calculation of Jones matrix components for given polarisation (eq. 4 and 5
+ // in the Sokolowski et al (2017) paper --------------------------------
+ // Internal function to calculate Jones matrix componenets for a given
+ // polarisation (pol). The are calculated as electric field vectors E_theta_mn
+ // (eq.4) and E_phi_mn (eq.5 in the Sokolowski et al 2017 paper). INPUT :
+ // phi,theta - FEKO convention coordinates (phi=90-azim), theta=za in
+ // radians already Q1_accum, Q2_accum, M_accum, N_accum, MabsM, Nmax -
+ // coefficients calculated earlier in CalcModes function for given
+ // frequency, delays and amplitudes
+ // OUTPUT :
+ // Jones matrix filled with components for given polarisation (pol input
+ // parameter)
+ void CalcSigmas(double phi, double theta, const Coefficients& coefficients,
+ char pol, JonesMatrix& jones_matrix) const;
+
+ std::complex<double> JPower(size_t i) const { return _jPowerTable[i % 4]; }
+ std::complex<double> _jPowerTable[4];
+
+ // Information on last modes parameters - not to recalculate the same again
+ // and again !
+ int _calcModesLastFreqHz;
+ std::vector<double> _calcModesLastDelays;
+ std::vector<double> _calcModesLastAmps;
+
+ // function comparing current parameters : frequency, delays and amplitudes
+ // with those previously used to calculate spherical waves coefficients
+ // (stored in the 3 variables above)
+ bool IsCalcModesRequired(int freq_hz, int n_ant, const double* delays,
+ const double* amps);
+
+ // Calculation of modes Q1 and Q2 and coefficients N and M and some derived
+ // variables (MabsM_X,MabsM_Y,Nmax_X and Nmax_Y) to make it once for a given
+ // pointing and then re-use for many different (azim,za) angles:
+
+ // function calculating coefficients for X and Y and storing parameters
+ // frequency, delays and amplitudes
+ void GetModes(int freq_hz, size_t n_ant, const double* delays,
+ const double* amps, Coefficients& coefsX, Coefficients& coefsY,
+ RecursiveLock<std::mutex>& lock);
+
+ // function calculating all coefficients Q1, Q2, N, M and derived MabsM, Nmax
+ // for a given polarisation ("X" or "Y") - perhaps enum should be used here
+ double CalcModes(int freq_hz, size_t n_ant, const double* delays,
+ const double* amp, char pol, Coefficients& coefs,
+ RecursiveLock<std::mutex>& lock);
+
+ // Calculation of normalisation matrix :
+ JonesMatrix CalcZenithNormMatrix(int freq_hz,
+ RecursiveLock<std::mutex>& lock);
+
+ std::map<int, JonesMatrix> _normJonesCache;
+
+ double _delays[16], _amps[16];
+
+ // ----------------------------------------------------------------
+ // Checking frequences included in the H5 file (stored in vector m_freq_list)
+ // :
+ bool has_freq(int freq_hz) const;
+ int find_closest_freq(int freq_hz) const;
+
+ // HDF5 File interface and data structures for H5 data
+ // Interface to HDF5 file format and structures to store H5 data
+ // Functions for reading H5 file and its datasets :
+ // Read data from H5 file, file name is specified in the object constructor
+ void Read();
+
+ const std::vector<std::vector<double>>& GetDataSet(
+ const DataSetIndex& index, RecursiveLock<std::mutex>& lock);
+
+ // Read dataset_name from H5 file
+ void ReadDataSet(const std::string& dataset_name,
+ std::vector<std::vector<double>>& out_vector,
+ H5::H5File& h5File);
+
+ // function for iteration call to H5Ovisit function :
+ static herr_t list_obj_iterate(hid_t loc_id, const char* name,
+ const H5O_info_t* info, void* operator_data);
+
+ std::unique_ptr<H5::H5File> _h5File;
+ std::string _searchPath;
+
+ // Data structures for H5 file data :
+ std::vector<std::string> m_obj_list; // list of datasets in H5 file
+ std::vector<int> m_freq_list; // list of simulated frequencies
+ std::vector<std::vector<double>> m_Modes; // data in Modes DataSet
+
+ FactorialTable _factorial;
+
+ // Calculations of Legendre polynomials :
+ static void lpmv(std::vector<double>& output, int n, double x);
+
+ // OUTPUT : returns list of Legendre polynomial values calculated up to order
+ // nmax :
+ static int P1sin(int nmax, double theta, std::vector<double>& p1sin_out,
+ std::vector<double>& p1_out);
+
+ std::map<DataSetIndex, std::vector<std::vector<double>>> _dataSetCache;
+
+ std::mutex mutex_;
+};
+} // namespace everybeam::mwabeam
+#endif
diff --git a/cpp/mwabeam/factorialtable.h b/cpp/mwabeam/factorialtable.h
new file mode 100644
index 0000000000000000000000000000000000000000..0c6e57496706148ebb941360f0e46d870f93890c
--- /dev/null
+++ b/cpp/mwabeam/factorialtable.h
@@ -0,0 +1,31 @@
+#ifndef EVERYBEAM_MWABEAM_FACTORIALTABLE_H_
+#define EVERYBEAM_MWABEAM_FACTORIALTABLE_H_
+
+#include <boost/math/special_functions/factorials.hpp>
+
+#include <aocommon/uvector.h>
+
+namespace everybeam::mwabeam {
+class FactorialTable {
+ public:
+ FactorialTable(size_t nprecalculated) : table_(nprecalculated) {
+ for (unsigned i = 0; i <= nprecalculated; i++)
+ table_[i] = boost::math::factorial<double>(i);
+ }
+
+ double operator()(unsigned n) const {
+ if (n >= table_.size()) {
+ return boost::math::factorial<double>(n);
+ }
+
+ if (n < table_.size())
+ return table_[n];
+ else
+ return boost::math::factorial<double>(n);
+ }
+
+ private:
+ aocommon::UVector<double> table_;
+};
+} // namespace everybeam::mwabeam
+#endif // EVERYBEAM_MWABEAM_FACTORIALTABLE_H_
diff --git a/cpp/mwabeam/recursivelock.h b/cpp/mwabeam/recursivelock.h
new file mode 100644
index 0000000000000000000000000000000000000000..cdcdcd21e9e7ebdd3e086ebdce125748f49e6ec9
--- /dev/null
+++ b/cpp/mwabeam/recursivelock.h
@@ -0,0 +1,105 @@
+#ifndef EVERYBEAM_MWABEAM_RECURSIVELOCK_H_
+#define EVERYBEAM_MWABEAM_RECURSIVELOCK_H_
+
+#include <cstring>
+#include <mutex>
+#include <system_error>
+
+namespace everybeam::mwabeam {
+
+template <typename Mutex>
+class RecursiveLock {
+ public:
+ RecursiveLock() : mutex_(nullptr), nlocks_(0) {}
+
+ RecursiveLock(RecursiveLock&& other)
+ : mutex_(other.mutex_), nlocks_(other.nlocks_) {
+ other.mutex_ = nullptr;
+ other.nlocks_ = 0;
+ }
+
+ RecursiveLock(Mutex& mutex) : mutex_(&mutex), nlocks_(1) { mutex_->lock(); }
+
+ RecursiveLock(Mutex& mutex, std::defer_lock_t) noexcept
+ : mutex_(&mutex), nlocks_(0) {}
+
+ RecursiveLock(Mutex& mutex, std::try_to_lock_t) : mutex_(&mutex), nlocks_(0) {
+ try_lock();
+ }
+
+ RecursiveLock(Mutex& mutex, std::adopt_lock_t) noexcept
+ : mutex_(&mutex), nlocks_(1) {}
+
+ ~RecursiveLock() {
+ if (nlocks_ != 0) mutex_->unlock();
+ }
+
+ RecursiveLock& operator=(const RecursiveLock& other) {
+ if (nlocks_ != 0) mutex_->unlock();
+ mutex_ = other.mutex_;
+ nlocks_ = other.nLocks;
+ other.mutex_ = nullptr;
+ other.nlocks_ = 0;
+ }
+
+ void lock() {
+ if (mutex_ == nullptr)
+ throw std::system_error(
+ std::make_error_code(std::errc::operation_not_permitted),
+ "lock() called on RecursiveLock without mutex");
+ if (nlocks_ == 0) mutex_->lock();
+ ++nlocks_;
+ }
+
+ bool try_lock() {
+ if (mutex_ == nullptr)
+ throw std::system_error(
+ std::make_error_code(std::errc::operation_not_permitted),
+ "lock() called on RecursiveLock without mutex");
+ if (nlocks_ == 0) {
+ if (mutex_->try_lock()) {
+ ++nlocks_;
+ return true;
+ } else {
+ return false;
+ }
+ } else {
+ ++nlocks_;
+ return true;
+ }
+ }
+
+ void unlock() {
+ nlocks_--;
+ if (mutex_ == nullptr)
+ throw std::system_error(
+ std::make_error_code(std::errc::operation_not_permitted),
+ "unlock() called on RecursiveLock without mutex");
+ else if (nlocks_ == 0)
+ mutex_->unlock();
+ }
+
+ void swap(RecursiveLock& other) noexcept {
+ std::swap(mutex_, other.mutex_);
+ std::swap(nlocks_, other.nlocks_);
+ }
+
+ bool owns_lock() const noexcept { return nlocks_ != 0; }
+
+ Mutex* release() noexcept {
+ Mutex* m = mutex_;
+ mutex_ = nullptr;
+ nlocks_ = 0;
+ return m;
+ }
+
+ Mutex* mutex() const noexcept { return mutex_; }
+
+ operator bool() const noexcept { return owns_lock(); }
+
+ private:
+ Mutex* mutex_;
+ size_t nlocks_;
+};
+} // namespace everybeam::mwabeam
+#endif // EVERYBEAM_MWABEAM_RECURSIVELOCK_H_
diff --git a/cpp/mwabeam/tilebeam2016.cc b/cpp/mwabeam/tilebeam2016.cc
new file mode 100644
index 0000000000000000000000000000000000000000..46c636584da95d649970a590431242038151982a
--- /dev/null
+++ b/cpp/mwabeam/tilebeam2016.cc
@@ -0,0 +1,78 @@
+#include "tilebeam2016.h"
+
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MCPosition.h>
+#include <casacore/measures/Measures/MeasConvert.h>
+
+namespace everybeam::mwabeam {
+
+TileBeam2016::TileBeam2016(const double* delays, bool frequency_interpolation,
+ const std::string& coeff_path)
+ : Beam2016Implementation(delays, nullptr, coeff_path),
+ frequency_interpolation_(frequency_interpolation),
+ mwa_lattitude_(-26.703319),
+ mwa_longitude_(116.67081),
+ mwa_height_(377.0) {}
+
+void TileBeam2016::ArrayResponse(casacore::MEpoch& time,
+ casacore::MPosition& array_position, double ra,
+ double dec, double frequency,
+ std::complex<double>* gain) {
+ casacore::MeasFrame frame(array_position, time);
+ const casacore::MDirection::Ref hadec_ref(casacore::MDirection::HADEC, frame);
+ const casacore::MDirection::Ref azelgeo_ref(casacore::MDirection::AZELGEO,
+ frame);
+ const casacore::MDirection::Ref j2000_ref(casacore::MDirection::J2000, frame);
+ casacore::MPosition wgs = casacore::MPosition::Convert(
+ array_position, casacore::MPosition::WGS84)();
+ double arr_lattitude =
+ wgs.getValue().getLat(); // ant1Pos.getValue().getLat();
+
+ casacore::MDirection::Convert j2000_to_hadec(j2000_ref, hadec_ref),
+ j2000_to_azelgeo(j2000_ref, azelgeo_ref);
+
+ ArrayResponse(ra, dec, j2000_ref, j2000_to_hadec, j2000_to_azelgeo,
+ arr_lattitude, frequency, gain);
+}
+
+void TileBeam2016::ArrayResponse(
+ double ra, double dec, const casacore::MDirection::Ref& j2000_ref,
+ casacore::MDirection::Convert& j2000_to_hadec,
+ casacore::MDirection::Convert& j2000_to_azelgeo, double arr_lattitude,
+ double frequency, std::complex<double>* gain) {
+ static const casacore::Unit rad_unit("rad");
+ casacore::MDirection image_dir(
+ casacore::MVDirection(casacore::Quantity(ra, rad_unit), // RA
+ casacore::Quantity(dec, rad_unit)), // DEC
+ j2000_ref);
+
+ // convert ra, dec to ha
+ casacore::MDirection hadec = j2000_to_hadec(image_dir);
+ double ha = hadec.getValue().get()[0];
+ double sin_lat = std::sin(arr_lattitude), cos_lat = std::cos(arr_lattitude);
+ double sin_dec = std::sin(dec), cos_dec = std::cos(dec);
+ double cos_ha = std::cos(ha);
+ double zenith_distance =
+ std::acos(sin_lat * sin_dec + cos_lat * cos_dec * cos_ha);
+ casacore::MDirection azel = j2000_to_azelgeo(image_dir);
+ double azimuth = azel.getValue().get()[0];
+ ArrayResponse(zenith_distance, azimuth, frequency, gain);
+}
+
+/**
+ * Get the full Jones matrix response of the tile including the dipole
+ * reponse and array factor incorporating any mutual coupling effects
+ * from the impedance matrix. freq in Hz.
+ */
+void TileBeam2016::GetTabulatedResponse(double az, double za, double freq,
+ std::complex<double>* result) {
+ // input are radians -> convert to degrees as implementation class expects :
+ double az_deg = az * (180.00 / M_PI);
+ double za_deg = za * (180.00 / M_PI);
+ JonesMatrix jones = CalcJones(az_deg, za_deg, freq, 1);
+ result[0] = jones.j00;
+ result[1] = jones.j01;
+ result[2] = jones.j10;
+ result[3] = jones.j11;
+}
+} // namespace everybeam::mwabeam
diff --git a/cpp/mwabeam/tilebeam2016.h b/cpp/mwabeam/tilebeam2016.h
new file mode 100644
index 0000000000000000000000000000000000000000..bcb542f649cc9b186a8b56aa024155855a63986d
--- /dev/null
+++ b/cpp/mwabeam/tilebeam2016.h
@@ -0,0 +1,78 @@
+#ifndef EVERYBEAM_MWABEAM_TILEBEAM2016_H_
+#define EVERYBEAM_MWABEAM_TILEBEAM2016_H_
+
+#include <complex>
+#include <map>
+#include <set>
+
+#include "beam2016implementation.h"
+
+#include <casacore/measures/Measures/MEpoch.h>
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MCDirection.h>
+
+namespace everybeam {
+namespace mwabeam {
+
+class TileBeam2016 : public Beam2016Implementation {
+ public:
+ TileBeam2016(const double *delays, bool frequency_interpolation,
+ const std::string &coeff_path);
+
+ /**
+ * @brief Array response
+ *
+ * @param time Time (MJD)
+ * @param array_position
+ * @param ra right ascension (rad)
+ * @param dec declination (rad)
+ * @param frequency Frequency (Hz)
+ * @param gain
+ */
+ void ArrayResponse(casacore::MEpoch &time,
+ casacore::MPosition &array_position, double ra, double dec,
+ double frequency, std::complex<double> *gain);
+
+ void ArrayResponse(double ra, double dec,
+ const casacore::MDirection::Ref &j2000_ref,
+ casacore::MDirection::Convert &j2000_to_hadecref,
+ casacore::MDirection::Convert &j2000_to_azelgeoref,
+ double arrLatitude, double frequency,
+ std::complex<double> *gain);
+
+ void ArrayResponse(double zenith_angle, double azimuth, double frequency,
+ std::complex<double> *gain) {
+ if (frequency_interpolation_)
+ GetInterpolatedResponse(azimuth, zenith_angle, frequency, gain);
+ else
+ GetTabulatedResponse(azimuth, zenith_angle, frequency, gain);
+ }
+
+ private:
+ bool frequency_interpolation_;
+
+ /**
+ * Get the full Jones matrix response of the tile including the dipole
+ * reponse and array factor incorporating any mutual coupling effects
+ * from the impedance matrix. freq in Hz.
+ */
+ void GetTabulatedResponse(double az, double za, double freq,
+ std::complex<double> *result);
+
+ /**
+ * Create a few tabulated responses and interpolated over these.
+ */
+ void GetInterpolatedResponse(double az, double za, double freq,
+ std::complex<double> *result) {
+ // Not implemented yet: just call normal function
+ GetTabulatedResponse(az, za, freq, result);
+ }
+
+ const double mwa_lattitude_; // Array latitude. degrees North
+ const double mwa_longitude_; // Array longitude. degrees East
+ const double mwa_height_; // Array altitude. meters above sea level
+};
+} // namespace mwabeam
+} // namespace everybeam
+#endif // EVERYBEAM_MWABEAM_TILEBEAM2016_H_
diff --git a/cpp/options.h b/cpp/options.h
index af9576dc6e171aacb5f176dbcd607e584b353c5d..f1fbc4b9258cb79bd9270ab3e304ee80d0506b3f 100644
--- a/cpp/options.h
+++ b/cpp/options.h
@@ -36,18 +36,25 @@ namespace everybeam {
class Options {
public:
Options()
- : use_differential_beam(false),
+ : coeff_path("."),
+ use_differential_beam(false),
use_channel_frequency(true),
- data_column_name("DATA"){};
+ data_column_name("DATA"),
+ frequency_interpolation(false){};
//! Default - empty - options class
static Options GetDefault() { return Options(); };
- // TODO? Specify path to element response coefficients file
- // std::string coeff_path;
+ // Path to coefficients file
+ std::string coeff_path;
+
+ // LOFAR specific
bool use_differential_beam;
bool use_channel_frequency;
std::string data_column_name;
+
+ // MWA specific (Lofar probably will follow)
+ bool frequency_interpolation;
};
} // namespace everybeam
#endif // EVERYBEAM_OPTIONS_H_
\ No newline at end of file
diff --git a/cpp/telescope/CMakeLists.txt b/cpp/telescope/CMakeLists.txt
index 23934e1b6a917fae0ae4f2d63a9d2dbf43d6644d..81dbc809e607f7d5fef94c736b663629852bf383 100644
--- a/cpp/telescope/CMakeLists.txt
+++ b/cpp/telescope/CMakeLists.txt
@@ -1,6 +1,7 @@
install (FILES
dish.h
- lofar.h
+ lofar.h
+ mwa.h
telescope.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/telescope")
diff --git a/cpp/telescope/atcatelescope.h b/cpp/telescope/atcatelescope.h
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/cpp/telescope/dish.cc b/cpp/telescope/dish.cc
index 0aeba4e28b53c1cdc9eaf1c685ce68573662ed58..ac5db7209011d4ba9c37f24cc56a84457af80be4 100644
--- a/cpp/telescope/dish.cc
+++ b/cpp/telescope/dish.cc
@@ -3,8 +3,9 @@
#include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
-using namespace everybeam;
-using namespace everybeam::telescope;
+using everybeam::griddedresponse::DishGrid;
+using everybeam::griddedresponse::GriddedResponse;
+using everybeam::telescope::Dish;
Dish::Dish(casacore::MeasurementSet &ms, const Options &options)
: Telescope(ms, options) {
@@ -20,10 +21,9 @@ Dish::Dish(casacore::MeasurementSet &ms, const Options &options)
}
}
-std::unique_ptr<griddedresponse::GriddedResponse> Dish::GetGriddedResponse(
+std::unique_ptr<GriddedResponse> Dish::GetGriddedResponse(
const coords::CoordinateSystem &coordinate_system) {
// Get and return GriddedResponse ptr
- std::unique_ptr<griddedresponse::GriddedResponse> grid(
- new griddedresponse::DishGrid(this, coordinate_system));
+ std::unique_ptr<GriddedResponse> grid(new DishGrid(this, coordinate_system));
return grid;
};
\ No newline at end of file
diff --git a/cpp/telescope/dish.h b/cpp/telescope/dish.h
index efa2a0b9ffcf03e7af7b4a2df9d144cabe1059c4..12d3cce33779b0f69bdb9486b278a356e6969673 100644
--- a/cpp/telescope/dish.h
+++ b/cpp/telescope/dish.h
@@ -1,4 +1,4 @@
-// dishtelescope.h: Base class for dish telescopes (VLA, ATCA, ...).
+// dish.h: Base class for dish telescopes (VLA, ATCA, ...).
// Inherits from Telescope class.
//
// Copyright (C) 2020
diff --git a/cpp/telescope/lofar.cc b/cpp/telescope/lofar.cc
index 263b5e01cd9b7318edaaf8437f43729f010f7561..489d4f915bbf569fd458f248c97c9d2799246b9a 100644
--- a/cpp/telescope/lofar.cc
+++ b/cpp/telescope/lofar.cc
@@ -8,27 +8,28 @@
#include <cassert>
#include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
-using namespace everybeam;
-using namespace everybeam::telescope;
-using namespace casacore;
+using everybeam::Station;
+using everybeam::griddedresponse::GriddedResponse;
+using everybeam::griddedresponse::LOFARGrid;
+using everybeam::telescope::LOFAR;
namespace {
bool GetPreappliedBeamDirection(casacore::MeasurementSet &ms,
- const std::string &dataColumnName,
- bool useDifferentialBeam,
- casacore::MDirection &preappliedBeamDir) {
+ const std::string &data_column_name,
+ bool use_differential_beam,
+ casacore::MDirection &preapplied_beam_dir) {
casacore::ScalarMeasColumn<casacore::MDirection> referenceDirColumn(
ms.field(),
casacore::MSField::columnName(casacore::MSFieldEnums::REFERENCE_DIR));
- preappliedBeamDir = referenceDirColumn(0);
+ preapplied_beam_dir = referenceDirColumn(0);
// Read beam keywords of input datacolumn
- casacore::ArrayColumn<std::complex<float>> dataCol(ms, dataColumnName);
- bool wasBeamApplied = false;
+ casacore::ArrayColumn<std::complex<float>> dataCol(ms, data_column_name);
+ bool was_beam_applied = false;
if (dataCol.keywordSet().isDefined("LOFAR_APPLIED_BEAM_MODE")) {
std::string mode = dataCol.keywordSet().asString("LOFAR_APPLIED_BEAM_MODE");
if (mode == "None")
- wasBeamApplied = false;
+ was_beam_applied = false;
else {
if (mode == "Element" || mode == "ArrayFactor")
throw std::runtime_error(
@@ -36,27 +37,27 @@ bool GetPreappliedBeamDirection(casacore::MeasurementSet &ms,
" beam. WSClean can only handle a full pre-applied beam (both "
"arrayfactor + element).");
else if (mode == "Full") {
- wasBeamApplied = true;
+ was_beam_applied = true;
casacore::String error;
casacore::MeasureHolder mHolder;
if (!mHolder.fromRecord(
error, dataCol.keywordSet().asRecord("LOFAR_APPLIED_BEAM_DIR")))
throw std::runtime_error(
"Error while reading LOFAR_APPLIED_BEAM_DIR keyword: " + error);
- preappliedBeamDir = mHolder.asMDirection();
+ preapplied_beam_dir = mHolder.asMDirection();
} else
throw std::runtime_error(
"Measurement set specifies an unknown beam correction: " + mode);
}
}
- if (wasBeamApplied || useDifferentialBeam) {
- useDifferentialBeam = true;
+ if (was_beam_applied || use_differential_beam) {
+ use_differential_beam = true;
}
- return useDifferentialBeam;
+ return use_differential_beam;
}
} // namespace
-LOFAR::LOFAR(MeasurementSet &ms, const ElementResponseModel model,
+LOFAR::LOFAR(casacore::MeasurementSet &ms, const ElementResponseModel model,
const Options &options)
: Telescope(ms, options) {
stations_.resize(nstations_);
@@ -64,11 +65,11 @@ LOFAR::LOFAR(MeasurementSet &ms, const ElementResponseModel 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));
+ casacore::MSAntenna antenna(ms.antenna());
+ casacore::MPosition::ScalarColumn antenna_pos_col(
+ antenna, antenna.columnName(casacore::MSAntennaEnums::POSITION));
+ casacore::MEpoch::ScalarColumn time_column(
+ ms, ms.columnName(casacore::MSMainEnums::TIME));
// Following is ms.field() related, first check whether field complies with
// LOFAR field
@@ -98,15 +99,15 @@ LOFAR::LOFAR(MeasurementSet &ms, const ElementResponseModel model,
.preapplied_beam_dir = preapplied_beam_dir};
}
-std::unique_ptr<griddedresponse::GriddedResponse> LOFAR::GetGriddedResponse(
+std::unique_ptr<GriddedResponse> LOFAR::GetGriddedResponse(
const coords::CoordinateSystem &coordinate_system) {
// Get and return GriddedResponse ptr
- std::unique_ptr<griddedresponse::GriddedResponse> grid(
- new griddedresponse::LOFARGrid(this, coordinate_system));
+ std::unique_ptr<GriddedResponse> grid(new LOFARGrid(this, coordinate_system));
return grid;
};
-Station::Ptr LOFAR::ReadStation(const MeasurementSet &ms, std::size_t id,
+Station::Ptr LOFAR::ReadStation(const casacore::MeasurementSet &ms,
+ std::size_t id,
const ElementResponseModel model) const {
Station::Ptr station = ReadLofarStation(ms, id, model);
return station;
diff --git a/cpp/telescope/lofar.h b/cpp/telescope/lofar.h
index 6e52d958eddd8661487630a606ac60253fd03c2b..0b63dfd0895597cb098f9f78c8c404e832de2d2d 100644
--- a/cpp/telescope/lofar.h
+++ b/cpp/telescope/lofar.h
@@ -1,4 +1,4 @@
-// LOFARTelescope.h: Base class for computing the response for the LOFAR
+// lofar.h: Base class for computing the response for the LOFAR
// telescope.
//
// Copyright (C) 2020
diff --git a/cpp/telescope/mwa.cc b/cpp/telescope/mwa.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f709a5ace998b3d17eb8f4dab4c32f5678df5f69
--- /dev/null
+++ b/cpp/telescope/mwa.cc
@@ -0,0 +1,32 @@
+#include "mwa.h"
+#include "../griddedresponse/mwagrid.h"
+
+#include <casacore/tables/Tables/TableRecord.h>
+#include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
+
+using everybeam::griddedresponse::GriddedResponse;
+using everybeam::griddedresponse::MWAGrid;
+using everybeam::telescope::MWA;
+
+MWA::MWA(casacore::MeasurementSet &ms, const Options &options)
+ : Telescope(ms, options) {
+ if (GetNrStations() == 0) throw std::runtime_error("No antennae in set");
+
+ casacore::MSAntenna antenna(ms.antenna());
+ casacore::MPosition::ScalarColumn antenna_pos_col(
+ antenna, antenna.columnName(casacore::MSAntennaEnums::POSITION));
+ ms_properties_.array_position = antenna_pos_col(0);
+
+ casacore::Table mwa_tile_pointing =
+ ms.keywordSet().asTable("MWA_TILE_POINTING");
+ casacore::ArrayColumn<int> delays_col(mwa_tile_pointing, "DELAYS");
+ casacore::Array<int> delays_arr = delays_col(0);
+ casacore::Array<int>::contiter delays_arr_ptr = delays_arr.cbegin();
+ for (int i = 0; i != 16; ++i) ms_properties_.delays[i] = delays_arr_ptr[i];
+}
+
+std::unique_ptr<GriddedResponse> MWA::GetGriddedResponse(
+ const coords::CoordinateSystem &coordinate_system) {
+ std::unique_ptr<GriddedResponse> grid(new MWAGrid(this, coordinate_system));
+ return grid;
+};
\ No newline at end of file
diff --git a/cpp/telescope/mwa.h b/cpp/telescope/mwa.h
new file mode 100644
index 0000000000000000000000000000000000000000..460e893b8edfb8475163cb517faac7ad51e4d69e
--- /dev/null
+++ b/cpp/telescope/mwa.h
@@ -0,0 +1,67 @@
+// mwa.h: Class for MWA telescopes.
+// Inherits from Telescope class.
+//
+// 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_MWA_H_
+#define EVERYBEAM_TELESCOPE_MWA_H_
+
+#include "telescope.h"
+
+#include <casacore/measures/Measures/MPosition.h>
+
+namespace everybeam {
+
+namespace griddedresponse {
+class MWAGrid;
+class GriddedResponse;
+} // namespace griddedresponse
+
+namespace telescope {
+
+class MWA final : public Telescope {
+ friend class griddedresponse::MWAGrid;
+
+ public:
+ /**
+ * @brief Construct a new MWA object
+ *
+ * @param ms MeasurementSet
+ * @param model Element Response model
+ * @param options telescope options
+ */
+ MWA(casacore::MeasurementSet &ms, const Options &options);
+
+ std::unique_ptr<griddedresponse::GriddedResponse> GetGriddedResponse(
+ const coords::CoordinateSystem &coordinate_system) override;
+
+ private:
+ struct MSProperties {
+ double delays[16];
+ casacore::MPosition array_position;
+ };
+ MSProperties ms_properties_;
+};
+
+} // namespace telescope
+} // namespace everybeam
+#endif // EVERYBEAM_TELESCOPE_MWA_H_
\ No newline at end of file
diff --git a/cpp/telescope/telescope.h b/cpp/telescope/telescope.h
index f0e58671b50b31c7d96fba6589c593a98e5b61ed..2dc075bc6cb96feb43e181305be005ea58df37e4 100644
--- a/cpp/telescope/telescope.h
+++ b/cpp/telescope/telescope.h
@@ -1,4 +1,4 @@
-// Telescope.h: Base class for computing the Telescope response.
+// telescope.h: Base class for computing the Telescope response.
//
// Copyright (C) 2020
// ASTRON (Netherlands Institute for Radio Astronomy)
diff --git a/cpp/telescope/vlatelescope.h b/cpp/telescope/vlatelescope.h
deleted file mode 100644
index 854d9a143f303e9eb12d06a5a98d0e7924fb9d98..0000000000000000000000000000000000000000
--- a/cpp/telescope/vlatelescope.h
+++ /dev/null
@@ -1 +0,0 @@
-#ifndef
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index 52cc4287b5bedc147e7f67b170988329af5f0a77..4917f084ffa1eabf00d7da6c9ce53c0e98ad82d0 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -1,7 +1,3 @@
-set(LOFAR_MOCK_MS CACHE STRING "LOFAR measurement set used for testing")
-set(VLA_MOCK_MS CACHE STRING "VLA measurement set used for testing")
-
-#------------------------------------------------------------------------------
# Find boost
find_package(Boost COMPONENTS unit_test_framework)
@@ -37,6 +33,16 @@ if(Boost_FOUND)
)
endif()
+ # Test MWA beamforming (assumes that mwa h5 coefficients file can be found!)
+ if( MWA_MOCK_MS AND MWA_COEFF_PATH )
+ set(TEST_SOURCES
+ ${TEST_SOURCES}
+ ${CMAKE_CURRENT_SOURCE_DIR}/tmwa.cc
+ )
+ elseif( MWA_MOCK_MS OR MWA_COEFF_PATH )
+ message(WARNING "To compile the MWA test, both the MWA_MOCK_MS and MWA_COEFF_PATH needs to be set. Now skipping MWA test compiling.")
+ endif()
+
# Add test executable
add_executable(runtests ${TEST_SOURCES})
target_link_libraries(runtests everybeam ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
diff --git a/cpp/test/tload_lofar.cc b/cpp/test/tload_lofar.cc
index 3c3b573ffc74ffad1526b22e6ef77d5d46399767..508b0a187f76f3024e101619605aed743d378336 100644
--- a/cpp/test/tload_lofar.cc
+++ b/cpp/test/tload_lofar.cc
@@ -69,7 +69,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
{0.108123, -5.36076e-05}};
// Compare with everybeam
- std::size_t offset_22 = (2 + 2 * height) * 4;
+ std::size_t offset_22 = (2 + 2 * width) * 4;
for (std::size_t i = 0; i < 4; ++i) {
BOOST_CHECK(std::abs(antenna_buffer_single[offset_22 + i] - lofar_p22[i]) <
1e-4);
@@ -82,7 +82,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
{0.0936919, 0.000110673}};
// Compare with everybeam
- std::size_t offset_13 = (1 + 3 * height) * 4;
+ std::size_t offset_13 = (1 + 3 * width) * 4;
for (std::size_t i = 0; i < 4; ++i) {
BOOST_CHECK(std::abs(antenna_buffer_single[offset_13 + i] - lofar_p13[i]) <
1e-4);
@@ -117,4 +117,9 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
norm_jones_mat += std::norm(antenna_buffer_diff_beam[offset_22 + i]);
}
BOOST_CHECK(std::abs(norm_jones_mat - 2.) < 1e-6);
+
+ // Print to np array
+ const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
+ npy::SaveArrayAsNumpy("lofar_station_responses.npy", false, 4, leshape,
+ antenna_buffer_single);
}
\ No newline at end of file
diff --git a/cpp/test/tmwa.cc b/cpp/test/tmwa.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a26179ae85c3f6ea38d1922c1e9d1dfcce92e68e
--- /dev/null
+++ b/cpp/test/tmwa.cc
@@ -0,0 +1,89 @@
+#include <boost/test/unit_test.hpp>
+
+#include "../load.h"
+#include "../options.h"
+#include "../griddedresponse/mwagrid.h"
+#include "../../external/npy.hpp"
+
+#include "config.h"
+#include <complex>
+#include <cmath>
+
+using namespace everybeam;
+
+BOOST_AUTO_TEST_CASE(load_mwa) {
+ Options options;
+ options.frequency_interpolation = false;
+ options.coeff_path = MWA_COEFF_PATH;
+
+ casacore::MeasurementSet ms(MWA_MOCK_MS);
+
+ std::unique_ptr<telescope::Telescope> telescope = Load(ms, options);
+
+ // Assert if we indeed have a MWA pointer
+ BOOST_CHECK(nullptr != dynamic_cast<telescope::MWA*>(telescope.get()));
+
+ // Assert if correct number of stations
+ std::size_t nstations = 128;
+ BOOST_CHECK_EQUAL(telescope->GetNrStations(), nstations);
+
+ double time = 4.87541808e+09;
+ double frequency = 133794999.99999999;
+ std::size_t width(16), height(16);
+ double ra(2.18166148), dec(-0.74612826), dl(1. * M_PI / 180.),
+ dm(1. * M_PI / 180.), shift_l(0.), shift_m(0.);
+ coords::CoordinateSystem coord_system = {.width = width,
+ .height = height,
+ .ra = ra,
+ .dec = dec,
+ .dl = dl,
+ .dm = dm,
+ .phase_centre_dl = shift_l,
+ .phase_centre_dm = shift_m};
+
+ std::unique_ptr<griddedresponse::GriddedResponse> grid_response =
+ telescope->GetGriddedResponse(coord_system);
+ BOOST_CHECK(nullptr !=
+ dynamic_cast<griddedresponse::MWAGrid*>(grid_response.get()));
+
+ std::vector<std::complex<float>> antenna_buffer(
+ grid_response->GetBufferSize(telescope->GetNrStations()));
+
+ grid_response->CalculateAllStations(antenna_buffer.data(), time, frequency,
+ 0);
+
+ // Check whether pixels are reproduced correctly at certain pixels on 16x16
+ // image MWABeam output at pixel (0, 8):
+ std::vector<std::complex<float>> mwa_p08 = {{-0.2311999, 0.0294384},
+ {-0.3220813, 0.03972803},
+ {-0.303604, 0.0423377},
+ {0.24611373, -0.03317199}};
+ // MWABeam output at pixel (10, 13):
+ std::vector<std::complex<float>> mwa_p1013 = {{-0.12771748, 0.02835142},
+ {-0.0854997, 0.01720171},
+ {-0.07546211, 0.01246299},
+ {0.13993847, -0.02099818}};
+ // MWABeam output at pixel (15, 15):
+ std::vector<std::complex<float>> mwa_p1515 = {{-0.0289045, 0.01751916},
+ {-0.01501623, 0.0077554},
+ {-0.01122625, 0.00331686},
+ {0.02796424, -0.00624412}};
+
+ // Convert pixel to buffer offsets
+ std::size_t offset_08 = (0 + 8 * width) * 4;
+ std::size_t offset_1013 = (10 + 13 * width) * 4;
+ std::size_t offset_1515 = (15 + 15 * width) * 4;
+
+ for (std::size_t i = 0; i < 4; ++i) {
+ BOOST_CHECK(std::abs(antenna_buffer[offset_08 + i] - mwa_p08[i]) < 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer[offset_1013 + i] - mwa_p1013[i]) <
+ 1e-6);
+ BOOST_CHECK(std::abs(antenna_buffer[offset_1515 + i] - mwa_p1515[i]) <
+ 1e-6);
+ }
+
+ // Print to np array, note: this spits out the transposed grid...
+ const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
+ npy::SaveArrayAsNumpy("mwa_station_responses.npy", false, 4, leshape,
+ antenna_buffer);
+}
\ No newline at end of file
diff --git a/cpp/test/tvla.cc b/cpp/test/tvla.cc
index 8711a8473093d411224eadac33bbd87908aff743..3a97092bd05977fd25a13fd0afccf56facdf86ab 100644
--- a/cpp/test/tvla.cc
+++ b/cpp/test/tvla.cc
@@ -17,13 +17,14 @@ BOOST_AUTO_TEST_CASE(load_vla) {
std::unique_ptr<telescope::Telescope> telescope = Load(ms);
- // Assert if we indeed have a LOFAR pointer
+ // Assert if we indeed have a VLA pointer
BOOST_CHECK(nullptr != dynamic_cast<telescope::Dish*>(telescope.get()));
// Assert if correct number of stations
std::size_t nstations = 25;
BOOST_CHECK_EQUAL(telescope->GetNrStations(), nstations);
+ // Time is irrelevant for dish telescope
double time = 0.5 * (4.90683119e+09 + 4.90684196e+09);
double frequency = 0.5e+09; // 1991000000.0;
std::size_t width(16), height(16);
@@ -69,10 +70,10 @@ BOOST_AUTO_TEST_CASE(load_vla) {
std::vector<std::complex<float>> vla_p1012 = {
{0.8672509, 0.}, {0, 0}, {0, 0}, {0.8672509, 0.}};
- // Compute offsets with everybeam
- std::size_t offset_00 = (0 + 0 * height) * 4;
- std::size_t offset_23 = (2 + 3 * height) * 4;
- std::size_t offset_1012 = (10 + 12 * height) * 4;
+ // Convert pixel to buffer offsets
+ std::size_t offset_00 = (0 + 0 * width) * 4;
+ std::size_t offset_23 = (2 + 3 * width) * 4;
+ std::size_t offset_1012 = (10 + 12 * width) * 4;
// Check if results are reproduced
BOOST_CHECK_EQUAL_COLLECTIONS(antenna_buffer.begin() + offset_00,
@@ -85,7 +86,7 @@ BOOST_AUTO_TEST_CASE(load_vla) {
antenna_buffer.begin() + offset_1012 + 4,
vla_p1012.begin(), vla_p1012.end());
- // Print to np array
+ // Print to np array, note: this spits out the transposed grid
const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
npy::SaveArrayAsNumpy("vla_station_responses.npy", false, 4, leshape,
antenna_buffer);
diff --git a/external/README.md b/external/README.md
index 7cd736aba8f8dfd06812e7dcefbaad4d75e7d022..18b169d710a50d55c412eccf60083686ccf4f981 100644
--- a/external/README.md
+++ b/external/README.md
@@ -1,9 +1,9 @@
External dependencies are included via git submodules, see .gitmodules in the root of the source.
-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:
+Please note that the dependencies are fixed to a specific commit, to make sure tests do not break as a
+result of silently updating the submodules. Dependencies are checked out on the following commits:
-- `aocommon`: 1f77c0c5e0d70507f227892353df8c5410cd87a4
+- `aocommon`: 05917d37e32a1bcd3fc985bc775b0b7b13d12532 (https://gitlab.com/aroffringa/aocommon)
- `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 138f28f4e8b125910dd2e4afca26bf25507161a5..05917d37e32a1bcd3fc985bc775b0b7b13d12532 160000
--- a/external/aocommon
+++ b/external/aocommon
@@ -1 +1 @@
-Subproject commit 138f28f4e8b125910dd2e4afca26bf25507161a5
+Subproject commit 05917d37e32a1bcd3fc985bc775b0b7b13d12532
diff --git a/scripts/clang-format-check.sh b/scripts/clang-format-check.sh
deleted file mode 100755
index a750894f39e445da27378ccdc9a0fdeab1cd03ba..0000000000000000000000000000000000000000
--- a/scripts/clang-format-check.sh
+++ /dev/null
@@ -1,36 +0,0 @@
-#!/bin/sh
-# Do check on clang-format, script is inspired on:
-# - https://dev.to/10xlearner/formatting-cpp-c-javascript-and-other-stuff-2pof
-# - https://gitlab.freedesktop.org/monado/monado/-/blob/master/scripts/format-and-spellcheck.sh
-#
-# Author: Jakob Maljaars
-# Email: jakob.maljaars_@_stcorp.nl
-#
-# To hook this script to pre-commit include the line "./scripts/clang-format-check.sh" to .git/hooks/pre-commit
-# and make sure pre-commit is an executable shell script.
-
-set -e
-PATCH_NAME=clang-fixes.diff
-
-echo "Running clang-format..."
-echo
-
-# Run clang-format from run-clang-format.sh
-SCRIPT_PATH=$(dirname "$0")
-$SCRIPT_PATH/run-clang-format.sh
-
-# Can't use tee because it hides the exit code
-if git diff --patch --exit-code > $PATCH_NAME; then
- echo
- # print in bold-face green
- echo "\e[1m\e[32mclang-format and codespell changed nothing."
- rm -f $PATCH_NAME
- exit 0;
-else
- echo
- # Print in bold-face red
- echo "\e[1m\e[31mclang-format made at least one change. Run clang-format before pushing!"
- rm -f $PATCH_NAME
- exit 1;
-fi
-
diff --git a/scripts/download_lofar_ms.sh b/scripts/download_lofar_ms.sh
index 487f878315e6d87e5cccf707ca75dafd993e5fd3..07f89f7169f4e080e40af8aa6c4e87eab230495f 100755
--- a/scripts/download_lofar_ms.sh
+++ b/scripts/download_lofar_ms.sh
@@ -19,7 +19,7 @@ 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
+ wget -q https://www.astron.nl/citt/EveryBeam/L258627-one-timestep.tar.bz2 -O $LOFAR_MOCK_ARCHIVE
fi
if [ -d $LOFAR_MOCK_MS ]
diff --git a/scripts/download_mwa_coeff.sh b/scripts/download_mwa_coeff.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a4bd4efe13ae80f673e9ff78e88ad950629de2d9
--- /dev/null
+++ b/scripts/download_mwa_coeff.sh
@@ -0,0 +1,26 @@
+#!/bin/sh
+#
+# Author: Jakob Maljaars
+# Email: jakob.maljaars_@_stcorp.nl
+
+# Script for downloading a mock VLA 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/
+
+MWA_COEFF_ARCHIVE=MWA_COEFF.tar.bz2
+MWA_COEFF_H5=MWA_COEFF.H5
+
+if [ ! -f "$MWA_COEFF_ARCHIVE" ]; then
+ wget -q www.astron.nl/citt/EveryBeam/mwa_full_embedded_element_pattern.tar.bz2 -O $MWA_COEFF_ARCHIVE
+fi
+
+tar -xf $MWA_COEFF_ARCHIVE
+rm $MWA_COEFF_ARCHIVE
\ No newline at end of file
diff --git a/scripts/download_mwa_ms.sh b/scripts/download_mwa_ms.sh
new file mode 100755
index 0000000000000000000000000000000000000000..80f9b7636cdf4bf4c32811f6ecfdf8867898825e
--- /dev/null
+++ b/scripts/download_mwa_ms.sh
@@ -0,0 +1,33 @@
+#!/bin/sh
+#
+# Author: Jakob Maljaars
+# Email: jakob.maljaars_@_stcorp.nl
+
+# Script for downloading a mock VLA 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/
+
+MWA_MOCK_ARCHIVE=MWA_ARCHIVE.tar.bz2
+MWA_MOCK_MS=MWA_MOCK.ms
+
+if [ ! -f "$MWA_MOCK_ARCHIVE" ]; then
+ wget -q www.astron.nl/citt/EveryBeam/MWA-single-timeslot.tar.bz2 -O $MWA_MOCK_ARCHIVE
+fi
+
+if [ -d $MWA_MOCK_MS ]
+then
+ echo "Directory already exists"
+else
+ mkdir $MWA_MOCK_MS
+fi
+
+tar -xf $MWA_MOCK_ARCHIVE -C $MWA_MOCK_MS --strip-components=1
+rm $MWA_MOCK_ARCHIVE
\ No newline at end of file
diff --git a/scripts/download_vla_ms.sh b/scripts/download_vla_ms.sh
index 7c18ca8ed70e25d788240fb0b9cddd7337bcef51..e07507977f619a5bd97ca8233a6417b64bf8a834 100755
--- a/scripts/download_vla_ms.sh
+++ b/scripts/download_vla_ms.sh
@@ -19,7 +19,7 @@ VLA_MOCK_ARCHIVE=VLA_ARCHIVE.tar.bz2
VLA_MOCK_MS=VLA_MOCK.ms
if [ ! -f "$VLA_MOCK_ARCHIVE" ]; then
- wget https://www.astron.nl/citt/EveryBeam/small-vla-set.tar.bz2 -O $VLA_MOCK_ARCHIVE
+ wget -q https://www.astron.nl/citt/EveryBeam/small-vla-set.tar.bz2 -O $VLA_MOCK_ARCHIVE
fi
if [ -d $VLA_MOCK_MS ]
diff --git a/scripts/run-clang-format.sh b/scripts/run-clang-format.sh
index 4191ca4ca414d11e978337aba83a6520beaee78f..3fd0419c0f1af1fb78f1dcceae3987d09dc46c0b 100755
--- a/scripts/run-clang-format.sh
+++ b/scripts/run-clang-format.sh
@@ -1,14 +1,55 @@
-#!/bin/sh
+#!/bin/bash
#
-# Author: Jakob Maljaars
-# Email: jakob.maljaars_@_stcorp.nl
+# run-clang-format.sh: Formats source code in this repo in accordance with .clang-format file.
+# This file is part of the DP3 software package.
+# Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
+# License: GNU General Public License version 3 or any later version
+#
+# To hook this script to pre-commit include the line
+# "./scripts/run-clang-format.sh" to .git/hooks/pre-commit
+# and make sure pre-commit is an executable shell script.
+
+#Script configuration for this repo. Adjust it when copying to a different repo.
+
+#The directory that contains the source files, which clang-format should format.
+SOURCE_DIR=$(dirname "$0")/..
+
+#Directories that must be excluded from formatting. These paths are
+#relative to SOURCE_DIR.
+EXCLUDE_DIRS=(external)
+
+#The extensions of the source files, which clang-format should format.
+SOURCE_EXT=(*.cc *.h)
+
+#End script configuration.
set -e
-SCRIPT_PATH=$(dirname "$0")
-cd $SCRIPT_PATH
-# Move up to parent folder which contains the source
-cd ..
+# print in bold-face
+echo -e "\e[1mRunning clang-format...\e[0m"
+
+# Convert SOURCE_EXT into "-name ext1 -o -name ext2 -o name ext3 ..."
+FIND_NAMES="-name ${SOURCE_EXT[0]}"
+for i in `seq 1 $((${#SOURCE_EXT[*]} - 1))`; do
+ FIND_NAMES+=" -o -name ${SOURCE_EXT[$i]}"
+done
+
+# Convert EXCLUDE_DIRS into "-path ./dir1 -prune -o -path ./dir2 -prune -o ..."
+FIND_EXCLUDES=
+for e in ${EXCLUDE_DIRS[*]}; do
+ FIND_EXCLUDES+="-path ./$e -prune -o "
+done
+
+cd $SOURCE_DIR
+find . $FIND_EXCLUDES -type f \( $FIND_NAMES \) \
+ -exec clang-format -i -style=file \{\} +
-# Find all cpp headers ("*.h") and code files ("*.cc") source tree and execute clang-format. Exclude ./external director
-find . -path ./external -prune -o -type f \( -name "*.h" -o -name "*.cc" \) -exec clang-format -i -style=file \{\} +
+if git diff --exit-code --quiet; then
+ # print in bold-face green
+ echo -e "\e[1m\e[32mGreat job, git shows no changed files!\e[0m"
+ exit 0;
+else
+ # Print in bold-face red
+ echo -e "\e[1m\e[31mGit shows at least one changed file now!\e[0m"
+ exit 1;
+fi