diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0265fb93bc18d19cc0322be530ded251107fe5c0..bcc721e533894259cb89b211d209e063b84d33c8 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -31,14 +31,15 @@ test-and-coverage:
- pip3 install gcovr
script:
- mkdir build
- # Download LOFAR Mock measurement set
- - ./scripts/download_lofar_ms.sh
- - cd test_data/LOFAR_MOCK.ms/ && export LOFAR_MOCK_MS=$PWD
- - cd ../../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
# 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" ..
+ - 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" ..
- make install -j8
- ctest -T test
+ # Explicit test that vla is checked
+ - ./cpp/test/runtests --log_level=test_suite --run_test=load_vla
- 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 7074966eb17c8e2639d9a5b1effdcf887ef385fc..29fe68a0093c1fd592cba6628021bd23fe03b328 100644
--- a/CMake/config.h.in
+++ b/CMake/config.h.in
@@ -4,5 +4,6 @@
#define EVERYBEAM_DATA_DIR "@CMAKE_INSTALL_DATA_DIR@"
#define TEST_MEASUREMENTSET "@TEST_MEASUREMENTSET@"
#define LOFAR_MOCK_MS "@LOFAR_MOCK_MS@"
+#define VLA_MOCK_MS "@VLA_MOCK_MS@"
#endif
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index a614621860a9e19419c3980fbdfc95bda370ccea..c1746ca68264a01002eb35586498c5c3adc1057f 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -6,6 +6,7 @@ add_subdirectory(hamaker)
add_subdirectory(lobes)
add_subdirectory(oskar)
add_subdirectory(telescope)
+add_subdirectory(circularsymmetric)
#------------------------------------------------------------------------------
add_library(everybeam SHARED
@@ -19,7 +20,11 @@ add_library(everybeam SHARED
lofarreadutils.cc
station.cc
telescope/lofar.cc
+ telescope/dish.cc
griddedresponse/lofargrid.cc
+ griddedresponse/dishgrid.cc
+ circularsymmetric/voltagepattern.cc
+ circularsymmetric/vlabeam.cc
)
target_include_directories(everybeam PUBLIC ${CASACORE_INCLUDE_DIR})
diff --git a/cpp/circularsymmetric/CMakeLists.txt b/cpp/circularsymmetric/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..1dbf27541ad0d67e7ac7abcc2d7f06af305a285a
--- /dev/null
+++ b/cpp/circularsymmetric/CMakeLists.txt
@@ -0,0 +1,4 @@
+install (FILES
+ vlabeam.h
+ voltagepattern.h
+DESTINATION "include/${CMAKE_PROJECT_NAME}/circularsymmetric")
\ No newline at end of file
diff --git a/cpp/circularsymmetric/vlabeam.cc b/cpp/circularsymmetric/vlabeam.cc
new file mode 100644
index 0000000000000000000000000000000000000000..e5d9145b958caa420ed987fe32ffe905e58100f1
--- /dev/null
+++ b/cpp/circularsymmetric/vlabeam.cc
@@ -0,0 +1,383 @@
+#include "vlabeam.h"
+
+#include <cmath>
+#include <vector>
+
+using namespace everybeam::circularsymmetric;
+
+std::array<double, 5> VLABeam::GetCoefficients(const std::string& bandName,
+ double freq) {
+ char band = '?';
+
+ size_t sharp = bandName.find('#');
+ if (sharp != std::string::npos) {
+ if (sharp > 5 && bandName.substr(0, 5) == "EVLA_") band = bandName[5];
+ }
+ if (band == '?') {
+ band = DetermineFeed(freq);
+ }
+ LimitFreqForBand(band, freq);
+
+ auto coeffmap = GetCoefficients();
+
+ const std::array<double, 5>* coeff;
+ double freqMHz = freq * 1e-6;
+ double freqUsed = freq * 1e-6;
+ std::map<int, std::array<double, 5>>::iterator low, prev;
+ low = coeffmap.lower_bound(freqMHz);
+ if (low == coeffmap.end()) {
+ --low;
+ coeff = &low->second;
+ freqUsed = low->first;
+ } else if (low == coeffmap.begin()) {
+ coeff = &low->second;
+ freqUsed = low->first;
+ } else {
+ prev = low;
+ --prev;
+ if (std::fabs(freqMHz - prev->first) < std::fabs(low->first - freqMHz)) {
+ coeff = &prev->second;
+ freqUsed = prev->first;
+ } else {
+ coeff = &low->second;
+ freqUsed = low->first;
+ }
+ }
+ return *coeff;
+}
+
+char VLABeam::DetermineFeed(double freq, double freqCenter) {
+ if ((freqCenter > 224e6 && freqCenter < 480e6) ||
+ (freq > 224e6 && freq < 480e6))
+ return 'P';
+ if ((freqCenter > 900e6 && freqCenter < 2003.0e6) ||
+ (freq > 900e6 && freq < 2003e6))
+ return 'L';
+ if ((freqCenter > 1990e6 && freqCenter < 4001.0e6) ||
+ (freq > 1990e6 && freq < 4001e6))
+ return 'S';
+ if ((freqCenter > 3990e6 && freqCenter < 8001.0e6) ||
+ (freq > 3990e6 && freq < 8001e6))
+ return 'C';
+ if ((freqCenter > 7990e6 && freqCenter < 12001.0e6) ||
+ (freq > 7990e6 && freq < 12001e6))
+ return 'X';
+ if ((freqCenter > 12000e6 && freqCenter < 18000.0e6) ||
+ (freq > 12000e6 && freq < 18000e6))
+ return 'U';
+ if ((freqCenter > 19000e6 && freqCenter < 26000.0e6) ||
+ (freq > 19000e6 && freq < 26000e6))
+ return 'K';
+ if ((freqCenter > 28000e6 && freqCenter < 38000.0e6) ||
+ (freq > 28000e6 && freq < 38000e6))
+ return 'A';
+ if ((freqCenter > 41000e6 && freqCenter < 50000.0e6) ||
+ (freq > 41000e6 && freq < 50000e6))
+ return 'Q';
+
+ return '?';
+}
+
+// From PBMath1DEVLA::limitFreqForBand
+void VLABeam::LimitFreqForBand(char band, double& freq) {
+ if (band == 'P') {
+ if (freq <= 224e6) freq = 232e6;
+ if (freq >= 480e6) freq = 470e6;
+ } else if (band == 'L') {
+ if (freq <= 900e6) freq = 1040e6;
+ if (freq >= 2000e6) freq = 2000e6;
+ } else if (band == 'S') {
+ if (freq < 2052e6) freq = 2052e6;
+ if (freq >= 3948e6) freq = 3948e6;
+ } else if (band == 'C') {
+ if (freq < 4052e6) freq = 4052e6;
+ if (freq >= 7948e6) freq = 7948e6;
+ } else if (band == 'X') {
+ if (freq < 8052e6) freq = 8052e6;
+ if (freq >= 11948e6) freq = 11948e6;
+ } else if (band == 'U') {
+ if (freq < 12052e6) freq = 12052e6;
+ if (freq >= 17948e6) freq = 17948e6;
+ } else if (band == 'K') {
+ if (freq < 19052e6) freq = 19052e6;
+ if (freq >= 25948e6) freq = 25948e6;
+ } else if (band == 'A') {
+ if (freq < 28052e6) freq = 28052e6;
+ if (freq >= 38048e6) freq = 38048e6;
+ } else if (band == 'Q') {
+ if (freq < 41052e6) freq = 41052e6;
+ if (freq >= 43948e6) freq = 43948e6;
+ }
+}
+
+std::map<char, double> VLABeam::GetFeedConf() {
+ std::map<char, double> feed_conf;
+ feed_conf['L'] =
+ (-185.9) * M_PI / 180.0; // squint orientation, rads, North of +AZ axis
+ feed_conf['S'] = (-11.61) * M_PI / 180.0;
+ feed_conf['C'] = (-104.8) * M_PI / 180.0;
+ feed_conf['X'] = (-113.7) * M_PI / 180.0;
+ feed_conf['U'] = (42.4) * M_PI / 180.0;
+ feed_conf['K'] = (64.4) * M_PI / 180.0;
+ feed_conf['A'] = (106.9) * M_PI / 180.0;
+ feed_conf['Q'] = (85.5) * M_PI / 180.0;
+ return feed_conf;
+}
+
+std::map<int, std::array<double, 5>> VLABeam::GetCoefficients() {
+ // This comes from PBMath1DEVLA::init()
+ // Also see https://library.nrao.edu/public/memos/evla/EVLAM_195.pdf
+ std::vector<double> wFreqs_p{232., 246., 281., 296., 312., 328., 344., 357.,
+ 382., 392., 403., 421., 458., 470., 1040, 1104,
+ 1168, 1232, 1296, 1360, 1424, 1488, 1552, 1680,
+ 1744, 1808, 1872, 1936, 2000};
+ for (double& f : wFreqs_p) f *= 1e6;
+
+ std::map<int, std::array<double, 5>> coeffmap;
+ ////P
+ coeffmap[232] = {1.0, -1.137e-3, 5.19e-7, -1.04e-10, 0.71e-14};
+ coeffmap[246] = {1.0, -1.130e-3, 5.04e-7, -1.02e-10, 0.77e-14};
+ coeffmap[281] = {1.0, -1.106e-3, 5.11e-7, -1.10e-10, 0.91e-14};
+ coeffmap[296] = {1.0, -1.125e-3, 5.27e-7, -1.14e-10, 0.96e-14};
+ coeffmap[312] = {1.0, -1.030e-3, 4.44e-7, -0.89e-10, 0.68e-14};
+ coeffmap[328] = {1.0, -0.980e-3, 4.25e-7, -0.87e-10, 0.69e-14};
+ coeffmap[344] = {1.0, -0.974e-3, 4.09e-7, -0.76e-10, 0.53e-14};
+ coeffmap[357] = {1.0, -0.996e-3, 4.23e-7, -0.79e-10, 0.51e-14};
+ coeffmap[382] = {1.0, -1.002e-3, 4.39e-7, -0.88e-10, 0.64e-14};
+ coeffmap[392] = {1.0, -1.067e-3, 5.13e-7, -1.12e-10, 0.90e-14};
+ coeffmap[403] = {1.0, -1.057e-3, 4.90e-7, -1.06e-10, 0.87e-14};
+ coeffmap[421] = {1.0, -1.154e-3, 5.85e-7, -1.33e-10, 1.08e-14};
+ coeffmap[458] = {1.0, -0.993e-3, 4.67e-7, -1.04e-10, 0.88e-14};
+ coeffmap[470] = {1.0, -1.010e-3, 4.85e-7, -1.07e-10, 0.86e-14};
+ /////////L
+ coeffmap[1040] = {1.000, -1.529e-3, 8.69e-7, -1.88e-10};
+ coeffmap[1104] = {1.000, -1.486e-3, 8.15e-7, -1.68e-10};
+ coeffmap[1168] = {1.000, -1.439e-3, 7.53e-7, -1.45e-10};
+ coeffmap[1232] = {1.000, -1.450e-3, 7.87e-7, -1.63e-10};
+ coeffmap[1296] = {1.000, -1.428e-3, 7.62e-7, -1.54e-10};
+ coeffmap[1360] = {1.000, -1.449e-3, 8.02e-7, -1.74e-10};
+ coeffmap[1424] = {1.000, -1.462e-3, 8.23e-7, -1.83e-10};
+ coeffmap[1488] = {1.000, -1.455e-3, 7.92e-7, -1.63e-10};
+ coeffmap[1552] = {1.000, -1.435e-3, 7.54e-7, -1.49e-10};
+ coeffmap[1680] = {1.000, -1.443e-3, 7.74e-7, -1.57e-10};
+ coeffmap[1744] = {1.000, -1.462e-3, 8.02e-7, -1.69e-10};
+ coeffmap[1808] = {1.000, -1.488e-3, 8.38e-7, -1.83e-10};
+ coeffmap[1872] = {1.000, -1.486e-3, 8.26e-7, -1.75e-10};
+ coeffmap[1936] = {1.000, -1.459e-3, 7.93e-7, -1.62e-10};
+ coeffmap[2000] = {1.000, -1.508e-3, 8.31e-7, -1.68e-10};
+ ////////S
+ coeffmap[2052] = {1.000, -1.429e-3, 7.52e-7, -1.47e-10};
+ coeffmap[2180] = {1.000, -1.389e-3, 7.06e-7, -1.33e-10};
+ coeffmap[2436] = {1.000, -1.377e-3, 6.90e-7, -1.27e-10};
+ coeffmap[2564] = {1.000, -1.381e-3, 6.92e-7, -1.26e-10};
+ coeffmap[2692] = {1.000, -1.402e-3, 7.23e-7, -1.40e-10};
+ coeffmap[2820] = {1.000, -1.433e-3, 7.62e-7, -1.54e-10};
+ coeffmap[2948] = {1.000, -1.433e-3, 7.46e-7, -1.42e-10};
+ coeffmap[3052] = {1.000, -1.467e-3, 8.05e-7, -1.70e-10};
+ coeffmap[3180] = {1.000, -1.497e-3, 8.38e-7, -1.80e-10};
+ coeffmap[3308] = {1.000, -1.504e-3, 8.37e-7, -1.77e-10};
+ coeffmap[3436] = {1.000, -1.521e-3, 8.63e-7, -1.88e-10};
+ coeffmap[3564] = {1.000, -1.505e-3, 8.37e-7, -1.75e-10};
+ coeffmap[3692] = {1.000, -1.521e-3, 8.51e-7, -1.79e-10};
+ coeffmap[3820] = {1.000, -1.534e-3, 8.57e-7, -1.77e-10};
+ coeffmap[3948] = {1.000, -1.516e-3, 8.30e-7, -1.66e-10};
+ /// C
+ coeffmap[4052] = {1.000, -1.406e-3, 7.41e-7, -1.48e-10};
+ coeffmap[4180] = {1.000, -1.385e-3, 7.09e-7, -1.36e-10};
+ coeffmap[4308] = {1.000, -1.380e-3, 7.08e-7, -1.37e-10};
+ coeffmap[4436] = {1.000, -1.362e-3, 6.95e-7, -1.35e-10};
+ coeffmap[4564] = {1.000, -1.365e-3, 6.92e-7, -1.31e-10};
+ coeffmap[4692] = {1.000, -1.339e-3, 6.56e-7, -1.17e-10};
+ coeffmap[4820] = {1.000, -1.371e-3, 7.06e-7, -1.40e-10};
+ coeffmap[4948] = {1.000, -1.358e-3, 6.91e-7, -1.34e-10};
+ coeffmap[5052] = {1.000, -1.360e-3, 6.91e-7, -1.33e-10};
+ coeffmap[5180] = {1.000, -1.353e-3, 6.74e-7, -1.25e-10};
+ coeffmap[5308] = {1.000, -1.359e-3, 6.82e-7, -1.27e-10};
+ coeffmap[5436] = {1.000, -1.380e-3, 7.05e-7, -1.37e-10};
+ coeffmap[5564] = {1.000, -1.376e-3, 6.99e-7, -1.31e-10};
+ coeffmap[5692] = {1.000, -1.405e-3, 7.39e-7, -1.47e-10};
+ coeffmap[5820] = {1.000, -1.394e-3, 7.29e-7, -1.45e-10};
+ coeffmap[5948] = {1.000, -1.428e-3, 7.57e-7, -1.57e-10};
+ coeffmap[6052] = {1.000, -1.445e-3, 7.68e-7, -1.50e-10};
+ coeffmap[6148] = {1.000, -1.422e-3, 7.38e-7, -1.38e-10};
+ coeffmap[6308] = {1.000, -1.463e-3, 7.94e-7, -1.62e-10};
+ coeffmap[6436] = {1.000, -1.478e-3, 8.22e-7, -1.74e-10};
+ coeffmap[6564] = {1.000, -1.473e-3, 8.00e-7, -1.62e-10};
+ coeffmap[6692] = {1.000, -1.455e-3, 7.76e-7, -1.53e-10};
+ coeffmap[6820] = {1.000, -1.487e-3, 8.22e-7, -1.72e-10};
+ coeffmap[6948] = {1.000, -1.472e-3, 8.05e-7, -1.67e-10};
+ coeffmap[7052] = {1.000, -1.470e-3, 8.01e-7, -1.64e-10};
+ coeffmap[7180] = {1.000, -1.503e-3, 8.50e-7, -1.84e-10};
+ coeffmap[7308] = {1.000, -1.482e-3, 8.19e-7, -1.72e-10};
+ coeffmap[7436] = {1.000, -1.498e-3, 8.22e-7, -1.66e-10};
+ coeffmap[7564] = {1.000, -1.490e-3, 8.18e-7, -1.66e-10};
+ coeffmap[7692] = {1.000, -1.481e-3, 7.98e-7, -1.56e-10};
+ coeffmap[7820] = {1.000, -1.474e-3, 7.94e-7, -1.57e-10};
+ coeffmap[7948] = {1.000, -1.448e-3, 7.69e-7, -1.51e-10};
+ //////X
+ coeffmap[8052] = {1.000, -1.403e-3, 7.21e-7, -1.37e-10};
+ coeffmap[8180] = {1.000, -1.398e-3, 7.10e-7, -1.32e-10};
+ coeffmap[8308] = {1.000, -1.402e-3, 7.16e-7, -1.35e-10};
+ coeffmap[8436] = {1.000, -1.400e-3, 7.12e-7, -1.32e-10};
+ coeffmap[8564] = {1.000, -1.391e-3, 6.95e-7, -1.25e-10};
+ coeffmap[8692] = {1.000, -1.409e-3, 7.34e-7, -1.49e-10};
+ coeffmap[8820] = {1.000, -1.410e-3, 7.36e-7, -1.45e-10};
+ coeffmap[8948] = {1.000, -1.410e-3, 7.34e-7, -1.43e-10};
+ coeffmap[9052] = {1.000, -1.403e-3, 7.20e-7, -1.36e-10};
+ coeffmap[9180] = {1.000, -1.396e-3, 7.09e-7, -1.31e-10};
+ coeffmap[9308] = {1.000, -1.432e-3, 7.68e-7, -1.55e-10};
+ coeffmap[9436] = {1.000, -1.414e-3, 7.43e-7, -1.47e-10};
+ coeffmap[9564] = {1.000, -1.416e-3, 7.45e-7, -1.47e-10};
+ coeffmap[9692] = {1.000, -1.406e-3, 7.26e-7, -1.39e-10};
+ coeffmap[9820] = {1.000, -1.412e-3, 7.36e-7, -1.43e-10};
+ coeffmap[9948] = {1.000, -1.409e-3, 7.29e-7, -1.39e-10};
+ coeffmap[10052] = {1.000, -1.421e-3, 7.46e-7, -1.45e-10};
+ coeffmap[10180] = {1.000, -1.409e-3, 7.25e-7, -1.36e-10};
+ coeffmap[10308] = {1.000, -1.402e-3, 7.13e-7, -1.31e-10};
+ coeffmap[10436] = {1.000, -1.399e-3, 7.09e-7, -1.29e-10};
+ coeffmap[10564] = {1.000, -1.413e-3, 7.37e-7, -1.43e-10};
+ coeffmap[10692] = {1.000, -1.412e-3, 7.34e-7, -1.41e-10};
+ coeffmap[10820] = {1.000, -1.401e-3, 7.12e-7, -1.31e-10};
+ coeffmap[10948] = {1.000, -1.401e-3, 7.12e-7, -1.31e-10};
+ coeffmap[10052] = {1.000, -1.401e-3, 7.12e-7, -1.31e-10};
+ coeffmap[11180] = {1.000, -1.394e-3, 6.99e-7, -1.24e-10};
+ coeffmap[11308] = {1.000, -1.394e-3, 7.01e-7, -1.26e-10};
+ coeffmap[11436] = {1.000, -1.391e-3, 6.94e-7, -1.22e-10};
+ coeffmap[11564] = {1.000, -1.389e-3, 6.92e-7, -1.22e-10};
+ coeffmap[11692] = {1.000, -1.386e-3, 6.80e-7, -1.15e-10};
+ coeffmap[11820] = {1.000, -1.391e-3, 6.88e-7, -1.19e-10};
+ coeffmap[11948] = {1.000, -1.399e-3, 6.97e-7, -1.22e-10};
+ /// U
+ coeffmap[12052] = {1.000, -1.399e-3, 7.17e-7, -1.34e-10};
+ coeffmap[12180] = {1.000, -1.392e-3, 7.07e-7, -1.31e-10};
+ coeffmap[12308] = {1.000, -1.393e-3, 7.19e-7, -1.38e-10};
+ coeffmap[12436] = {1.000, -1.393e-3, 7.20e-7, -1.40e-10};
+ coeffmap[12564] = {1.000, -1.395e-3, 7.19e-7, -1.38e-10};
+ coeffmap[12692] = {1.000, -1.397e-3, 7.20e-7, -1.37e-10};
+ coeffmap[12820] = {1.000, -1.388e-3, 7.06e-7, -1.32e-10};
+ coeffmap[12948] = {1.000, -1.397e-3, 7.18e-7, -1.36e-10};
+ coeffmap[13052] = {1.000, -1.400e-3, 7.27e-7, -1.40e-10};
+ coeffmap[13180] = {1.000, -1.406e-3, 7.44e-7, -1.50e-10};
+ coeffmap[13308] = {1.000, -1.403e-3, 7.37e-7, -1.47e-10};
+ coeffmap[13436] = {1.000, -1.392e-3, 7.08e-7, -1.31e-10};
+ coeffmap[13564] = {1.000, -1.384e-3, 6.94e-7, -1.24e-10};
+ coeffmap[13692] = {1.000, -1.382e-3, 6.95e-7, -1.25e-10};
+ coeffmap[13820] = {1.000, -1.376e-3, 6.88e-7, -1.24e-10};
+ coeffmap[13948] = {1.000, -1.384e-3, 6.98e-7, -1.28e-10};
+ coeffmap[14052] = {1.000, -1.400e-3, 7.36e-7, -1.48e-10};
+ coeffmap[14180] = {1.000, -1.397e-3, 7.29e-7, -1.45e-10};
+ coeffmap[14308] = {1.000, -1.399e-3, 7.32e-7, -1.45e-10};
+ coeffmap[14436] = {1.000, -1.396e-3, 7.25e-7, -1.42e-10};
+ coeffmap[14564] = {1.000, -1.393e-3, 7.20e-7, -1.39e-10};
+ coeffmap[14692] = {1.000, -1.384e-3, 7.03e-7, -1.31e-10};
+ coeffmap[14820] = {1.000, -1.388e-3, 7.06e-7, -1.32e-10};
+ coeffmap[14948] = {1.000, -1.393e-3, 7.16e-7, -1.37e-10};
+ coeffmap[15052] = {1.000, -1.402e-3, 7.38e-7, -1.48e-10};
+ coeffmap[15180] = {1.000, -1.407e-3, 7.47e-7, -1.53e-10};
+ coeffmap[15308] = {1.000, -1.406e-3, 7.41e-7, -1.48e-10};
+ coeffmap[15436] = {1.000, -1.399e-3, 7.31e-7, -1.44e-10};
+ coeffmap[15564] = {1.000, -1.397e-3, 7.28e-7, -1.43e-10};
+ coeffmap[15692] = {1.000, -1.401e-3, 7.35e-7, -1.46e-10};
+ coeffmap[15820] = {1.000, -1.402e-3, 7.34e-7, -1.45e-10};
+ coeffmap[15948] = {1.000, -1.399e-3, 7.30e-7, -1.44e-10};
+ coeffmap[16052] = {1.000, -1.419e-3, 7.59e-7, -1.54e-10};
+ coeffmap[16180] = {1.000, -1.419e-3, 7.59e-7, -1.52e-10};
+ coeffmap[16308] = {1.000, -1.412e-3, 7.40e-7, -1.44e-10};
+ coeffmap[16436] = {1.000, -1.407e-3, 7.32e-7, -1.40e-10};
+ coeffmap[16564] = {1.000, -1.408e-3, 7.32e-7, -1.41e-10};
+ coeffmap[16692] = {1.000, -1.410e-3, 7.34e-7, -1.40e-10};
+ coeffmap[16820] = {1.000, -1.407e-3, 7.27e-7, -1.38e-10};
+ coeffmap[16948] = {1.000, -1.423e-3, 7.63e-7, -1.55e-10};
+ coeffmap[17052] = {1.000, -1.437e-3, 7.87e-7, -1.66e-10};
+ coeffmap[17180] = {1.000, -1.438e-3, 7.84e-7, -1.64e-10};
+ coeffmap[17308] = {1.000, -1.445e-3, 7.98e-7, -1.71e-10};
+ coeffmap[17436] = {1.000, -1.452e-3, 8.10e-7, -1.77e-10};
+ coeffmap[17564] = {1.000, -1.458e-3, 8.13e-7, -1.70e-10};
+ coeffmap[17692] = {1.000, -1.456e-3, 8.06e-7, -1.72e-10};
+ coeffmap[17820] = {1.000, -1.453e-3, 8.00e-7, -1.68e-10};
+ coeffmap[17948] = {1.000, -1.452e-3, 7.99e-7, -1.69e-10};
+ /////K
+ coeffmap[19052] = {1.000, -1.419e-3, 7.56e-7, -1.53e-10};
+ coeffmap[19180] = {1.000, -1.426e-3, 7.70e-7, -1.59e-10};
+ coeffmap[19308] = {1.000, -1.433e-3, 7.82e-7, -1.64e-10};
+ coeffmap[19436] = {1.000, -1.429e-3, 7.73e-7, -1.60e-10};
+ coeffmap[19564] = {1.000, -1.427e-3, 7.70e-7, -1.59e-10};
+ coeffmap[19692] = {1.000, -1.425e-3, 7.65e-7, -1.56e-10};
+ coeffmap[19820] = {1.000, -1.430e-3, 7.76e-7, -1.62e-10};
+ coeffmap[19948] = {1.000, -1.434e-3, 7.81e-7, -1.63e-10};
+ coeffmap[21052] = {1.000, -1.448e-3, 8.05e-7, -1.73e-10};
+ coeffmap[21180] = {1.000, -1.436e-3, 7.84e-7, -1.63e-10};
+ coeffmap[21308] = {1.000, -1.441e-3, 7.94e-7, -1.68e-10};
+ coeffmap[21436] = {1.000, -1.439e-3, 7.89e-7, -1.66e-10};
+ coeffmap[21564] = {1.000, -1.442e-3, 7.96e-7, -1.69e-10};
+ coeffmap[21692] = {1.000, -1.435e-3, 7.81e-7, -1.61e-10};
+ coeffmap[21820] = {1.000, -1.442e-3, 7.92e-7, -1.66e-10};
+ coeffmap[21948] = {1.000, -1.439e-3, 7.82e-7, -1.61e-10};
+ coeffmap[23052] = {1.000, -1.401e-3, 7.21e-7, -1.37e-10};
+ coeffmap[23180] = {1.000, -1.408e-3, 7.31e-7, -1.41e-10};
+ coeffmap[23308] = {1.000, -1.407e-3, 7.28e-7, -1.39e-10};
+ coeffmap[23436] = {1.000, -1.407e-3, 7.31e-7, -1.41e-10};
+ coeffmap[23564] = {1.000, -1.419e-3, 7.47e-7, -1.47e-10};
+ coeffmap[23692] = {1.000, -1.395e-3, 7.10e-7, -1.33e-10};
+ coeffmap[23820] = {1.000, -1.413e-3, 7.36e-7, -1.42e-10};
+ coeffmap[23948] = {1.000, -1.402e-3, 7.21e-7, -1.36e-10};
+ coeffmap[25052] = {1.000, -1.402e-3, 7.17e-7, -1.31e-10};
+ coeffmap[25180] = {1.000, -1.432e-3, 7.73e-7, -1.58e-10};
+ coeffmap[25308] = {1.000, -1.407e-3, 7.22e-7, -1.33e-10};
+ coeffmap[25436] = {1.000, -1.417e-3, 7.43e-7, -1.45e-10};
+ coeffmap[25564] = {1.000, -1.422e-3, 7.52e-7, -1.48e-10};
+ coeffmap[25692] = {1.000, -1.427e-3, 7.59e-7, -1.52e-10};
+ coeffmap[25820] = {1.000, -1.416e-3, 7.42e-7, -1.44e-10};
+ coeffmap[25948] = {1.000, -1.422e-3, 7.46e-7, -1.45e-10};
+ /// A
+ coeffmap[28052] = {1.000, -1.444e-3, 7.61e-7, -1.44e-10};
+ coeffmap[28180] = {1.000, -1.439e-3, 7.54e-7, -1.42e-10};
+ coeffmap[28308] = {1.000, -1.457e-3, 7.87e-7, -1.58e-10};
+ coeffmap[28436] = {1.000, -1.457e-3, 7.90e-7, -1.60e-10};
+ coeffmap[28564] = {1.000, -1.455e-3, 7.87e-7, -1.59e-10};
+ coeffmap[28692] = {1.000, -1.458e-3, 7.88e-7, -1.58e-10};
+ coeffmap[28820] = {1.000, -1.453e-3, 7.81e-7, -1.56e-10};
+ coeffmap[28948] = {1.000, -1.460e-3, 7.98e-7, -1.64e-10};
+ coeffmap[31052] = {1.000, -1.415e-3, 7.44e-7, -1.44e-10};
+ coeffmap[31180] = {1.000, -1.408e-3, 7.26e-7, -1.37e-10};
+ coeffmap[31308] = {1.000, -1.413e-3, 7.28e-7, -1.36e-10};
+ coeffmap[31436] = {1.000, -1.394e-3, 7.07e-7, -1.30e-10};
+ coeffmap[31564] = {1.000, -1.404e-3, 7.23e-7, -1.37e-10};
+ coeffmap[31692] = {1.000, -1.427e-3, 7.48e-7, -1.44e-10};
+ coeffmap[31820] = {1.000, -1.418e-3, 7.48e-7, -1.48e-10};
+ coeffmap[31948] = {1.000, -1.413e-3, 7.37e-7, -1.42e-10};
+ coeffmap[34052] = {1.000, -1.42e-3, 7.28e-7, -1.34e-10};
+ coeffmap[34180] = {1.000, -1.46e-3, 7.77e-7, -1.53e-10};
+ coeffmap[34308] = {1.000, -1.42e-3, 7.41e-7, -1.42e-10};
+ coeffmap[34436] = {1.000, -1.42e-3, 7.36e-7, -1.39e-10};
+ coeffmap[34564] = {1.000, -1.46e-3, 7.76e-7, -1.52e-10};
+ coeffmap[34692] = {1.000, -1.42e-3, 7.34e-7, -1.38e-10};
+ coeffmap[34820] = {1.000, -1.42e-3, 7.34e-7, -1.39e-10};
+ coeffmap[34948] = {1.000, -1.45e-3, 7.68e-7, -1.49e-10};
+ coeffmap[37152] = {1.000, -1.42e-3, 7.47e-7, -1.44e-10};
+ coeffmap[37280] = {1.000, -1.41e-3, 7.35e-7, -1.40e-10};
+ coeffmap[37408] = {1.000, -1.45e-3, 7.65e-7, -1.46e-10};
+ coeffmap[37536] = {1.000, -1.41e-3, 7.13e-7, -1.29e-10};
+ coeffmap[37664] = {1.000, -1.41e-3, 7.30e-7, -1.38e-10};
+ coeffmap[37792] = {1.000, -1.45e-3, 7.75e-7, -1.50e-10};
+ coeffmap[37820] = {1.000, -1.45e-3, 7.68e-7, -1.49e-10};
+ coeffmap[38048] = {1.000, -1.41e-3, 7.38e-7, -1.43e-10};
+ // Q
+ coeffmap[41052] = {1.000, -1.453e-3, 7.69e-7, -1.47e-10};
+ coeffmap[41180] = {1.000, -1.479e-3, 8.03e-7, -1.61e-10};
+ coeffmap[41308] = {1.000, -1.475e-3, 7.97e-7, -1.58e-10};
+ coeffmap[41436] = {1.000, -1.451e-3, 7.73e-7, -1.51e-10};
+ coeffmap[41564] = {1.000, -1.450e-3, 7.71e-7, -1.51e-10};
+ coeffmap[41692] = {1.000, -1.465e-3, 7.79e-7, -1.49e-10};
+ coeffmap[41820] = {1.000, -1.460e-3, 7.73e-7, -1.47e-10};
+ coeffmap[41948] = {1.000, -1.434e-3, 7.47e-7, -1.40e-10};
+ coeffmap[43052] = {1.000, -1.428e-3, 7.40e-7, -1.38e-10};
+ coeffmap[43180] = {1.000, -1.418e-3, 7.29e-7, -1.34e-10};
+ coeffmap[43308] = {1.000, -1.433e-3, 7.49e-7, -1.43e-10};
+ coeffmap[43436] = {1.000, -1.438e-3, 7.55e-7, -1.45e-10};
+ coeffmap[43564] = {1.000, -1.419e-3, 7.36e-7, -1.40e-10};
+ coeffmap[43692] = {1.000, -1.397e-3, 7.13e-7, -1.33e-10};
+ coeffmap[43820] = {1.000, -1.423e-3, 7.39e-7, -1.40e-10};
+ coeffmap[43948] = {1.000, -1.452e-3, 7.68e-7, -1.47e-10};
+ return coeffmap;
+}
diff --git a/cpp/circularsymmetric/vlabeam.h b/cpp/circularsymmetric/vlabeam.h
new file mode 100644
index 0000000000000000000000000000000000000000..734c6484344a6fd062566c4a9b18cf78abc94b99
--- /dev/null
+++ b/cpp/circularsymmetric/vlabeam.h
@@ -0,0 +1,23 @@
+#ifndef EVERYBEAM_CIRCULARSYMMETRIC_VLABEAM_H_
+#define EVERYBEAM_CIRCULARSYMMETRIC_VLABEAM_H_
+
+#include <array>
+#include <map>
+#include <string>
+
+namespace everybeam {
+namespace circularsymmetric {
+class VLABeam {
+ public:
+ static std::array<double, 5> GetCoefficients(const std::string& band_name,
+ double freq);
+
+ private:
+ static std::map<int, std::array<double, 5>> GetCoefficients();
+ static std::map<char, double> GetFeedConf();
+ static char DetermineFeed(double freq, double freqCenter = 0.0);
+ static void LimitFreqForBand(char band, double& freq);
+};
+} // namespace circularsymmetric
+} // namespace everybeam
+#endif // EVERYBEAM_CIRCULARSYMMETRIC_VLABEAM_H_
\ No newline at end of file
diff --git a/cpp/circularsymmetric/voltagepattern.cc b/cpp/circularsymmetric/voltagepattern.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5d461f702a84a7a00ab5130e08499d9bd6112905
--- /dev/null
+++ b/cpp/circularsymmetric/voltagepattern.cc
@@ -0,0 +1,195 @@
+#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;
+
+void VoltagePattern::EvaluatePolynomial(const UVector<double>& coefficients,
+ bool invert) {
+ // This comes from casa's: void PBMath1DIPoly::fillPBArray(), wideband case
+ size_t nsamples = 10000;
+ size_t nfreq = frequencies_.size();
+ size_t ncoef = coefficients.size() / nfreq;
+ values_.resize(nsamples * nfreq);
+ inverse_increment_radius_ = double(nsamples - 1) / maximum_radius_arc_min_;
+ double* output = values_.data();
+ for (size_t n = 0; n != nfreq; n++) {
+ const double* freqcoefficients = &coefficients[n * ncoef];
+ for (size_t i = 0; i < nsamples; i++) {
+ double taper = 0.0;
+ double x2 = double(i) / inverse_increment_radius_;
+ x2 = x2 * x2;
+ double y = 1.0;
+
+ for (size_t j = 0; j < ncoef; j++) {
+ taper += y * freqcoefficients[j];
+ y *= x2;
+ }
+ if (taper >= 0.0) {
+ if (invert && taper != 0.0) {
+ taper = 1.0 / std::sqrt(taper);
+ } else {
+ taper = std::sqrt(taper);
+ }
+ } else {
+ taper = 0.0;
+ }
+ *output = taper;
+ ++output;
+ }
+ }
+};
+
+UVector<double> VoltagePattern::InterpolateValues(double freq) const {
+ UVector<double> result;
+ size_t ifit = 0;
+ size_t nfreq = frequencies_.size();
+ for (ifit = 0; ifit != nfreq; ifit++) {
+ if (freq <= frequencies_[ifit]) break;
+ }
+ if (ifit == 0) {
+ result.assign(values_.begin(), values_.begin() + NSamples());
+ } else if (ifit == nfreq) {
+ result.assign(values_.begin() + (nfreq - 1) * NSamples(), values_.end());
+ } else {
+ size_t n = NSamples();
+ double l = (freq - frequencies_[ifit - 1]) /
+ (frequencies_[ifit] - frequencies_[ifit - 1]);
+ const double* vpA = FreqIndexValues(ifit - 1);
+ const double* vpB = FreqIndexValues(ifit);
+ result.resize(n);
+ for (size_t i = 0; i != n; ++i) {
+ result[i] = vpA[i] * (1.0 - l) + vpB[i] * l;
+ }
+ }
+ return result;
+}
+
+const double* VoltagePattern::InterpolateValues(
+ double frequency_hz, UVector<double>& interpolated_values) const {
+ if (frequencies_.size() > 1) {
+ interpolated_values = InterpolateValues(frequency_hz);
+ return interpolated_values.data();
+ } else {
+ return FreqIndexValues(0);
+ }
+}
+
+double VoltagePattern::LmMaxSquared(double frequency_hz) const {
+ double factor =
+ (180.0 / M_PI) * 60.0 * frequency_hz * 1.0e-9; // arcminutes * GHz
+ double rmax = maximum_radius_arc_min_ / factor;
+ return rmax * rmax;
+}
+
+// void VoltagePattern::Render(PrimaryBeamImageSet& beamImages, double
+// pixel_scale_x,
+// double pixel_scale_y, double phase_centre_ra,
+// double phase_centre_dec, double pointing_ra,
+// double pointing_dec, double phase_centre_dl,
+// double phase_centre_dm, double frequency_hz)
+// const {
+// size_t width = beamImages.Width(), height = beamImages.Height();
+// double lmMaxSq = LmMaxSquared(frequency_hz);
+
+// UVector<double> interpolated_values;
+// const double* vp = InterpolateValues(frequency_hz, interpolated_values);
+
+// double factor =
+// (180.0 / M_PI) * 60.0 * frequency_hz * 1.0e-9; // arcminutes * GHz
+// double l0, m0;
+// ImageCoordinates::RaDecToLM(pointing_ra, pointing_dec, phase_centre_ra,
+// phase_centre_dec, l0, m0);
+// l0 += phase_centre_dl;
+// m0 += phase_centre_dm;
+// size_t imgIndex = 0;
+// // Logger::Debug << "Interpolating 1D voltage pattern to output
+// image...\n";
+// for (size_t iy = 0; iy != height; ++iy) {
+// for (size_t ix = 0; ix != width; ++ix) {
+// double l, m, ra, dec;
+// ImageCoordinates::XYToLM(ix, iy, pixel_scale_x, pixel_scale_y, width,
+// height,
+// l, m);
+// l += phase_centre_dl;
+// m += m0;
+// ImageCoordinates::LMToRaDec(l, m, phase_centre_ra, phase_centre_dec,
+// ra, dec); ImageCoordinates::RaDecToLM(ra, dec, pointing_ra,
+// pointing_dec, l, m); l -= l0; m -= m0; double r2 = l * l + m * m;
+// double out; if (r2 > lmMaxSq) {
+// out = 0.0;
+// } else {
+// double r = std::sqrt(r2) * factor;
+// int indx = int(r * inverse_increment_radius_);
+// out = vp[indx];
+// }
+
+// beamImages[0][imgIndex] = out;
+// beamImages[1][imgIndex] = 0.0;
+// beamImages[2][imgIndex] = 0.0;
+// beamImages[3][imgIndex] = 0.0;
+// beamImages[4][imgIndex] = 0.0;
+// beamImages[5][imgIndex] = 0.0;
+// beamImages[6][imgIndex] = out;
+// beamImages[7][imgIndex] = 0.0;
+// ++imgIndex;
+// }
+// }
+// }
+
+void VoltagePattern::Render(std::complex<float>* aterm, size_t width,
+ size_t height, double pixel_scale_x,
+ double pixel_scale_y, double phase_centre_ra,
+ double phase_centre_dec, double pointing_ra,
+ double pointing_dec, double phase_centre_dl,
+ double phase_centre_dm, double frequency_hz) const {
+ double lmMaxSq = LmMaxSquared(frequency_hz);
+
+ UVector<double> interpolated_values;
+ const double* vp = InterpolateValues(frequency_hz, interpolated_values);
+
+ double factor =
+ (180.0 / M_PI) * 60.0 * frequency_hz * 1.0e-9; // arcminutes * GHz
+ double l0, m0;
+ ImageCoordinates::RaDecToLM(pointing_ra, pointing_dec, phase_centre_ra,
+ phase_centre_dec, l0, m0);
+ l0 += phase_centre_dl;
+ m0 += phase_centre_dm;
+ for (size_t iy = 0; iy != height; ++iy) {
+ std::complex<float>* row = aterm + iy * width * 4;
+ for (size_t ix = 0; ix != width; ++ix) {
+ double l, m, ra, dec;
+ ImageCoordinates::XYToLM(ix, iy, pixel_scale_x, pixel_scale_y, width,
+ height, l, m);
+ l += phase_centre_dl;
+ m += m0;
+ ImageCoordinates::LMToRaDec(l, m, phase_centre_ra, phase_centre_dec, ra,
+ dec);
+ ImageCoordinates::RaDecToLM(ra, dec, pointing_ra, pointing_dec, l, m);
+ l -= l0;
+ m -= m0;
+ double r2 = l * l + m * m;
+ double out;
+ if (r2 > lmMaxSq) {
+ out = 1e-4;
+ } else {
+ double r = std::sqrt(r2) * factor;
+ int indx = int(r * inverse_increment_radius_);
+ out = vp[indx] * (1.0 - 1e-4) + 1e-4;
+ }
+
+ std::complex<float>* ptr = row + ix * 4;
+ ptr[0] = out;
+ ptr[1] = 0.0;
+ ptr[2] = 0.0;
+ ptr[3] = out;
+ }
+ }
+}
diff --git a/cpp/circularsymmetric/voltagepattern.h b/cpp/circularsymmetric/voltagepattern.h
new file mode 100644
index 0000000000000000000000000000000000000000..756f98a3c1e59f5e4808799c3375fa5d60237559
--- /dev/null
+++ b/cpp/circularsymmetric/voltagepattern.h
@@ -0,0 +1,61 @@
+#ifndef EVERYBEAM_CIRCULARSYMMETRIC_VOLTAGE_PATTERN_H_
+#define EVERYBEAM_CIRCULARSYMMETRIC_VOLTAGE_PATTERN_H_
+
+#include <aocommon/uvector.h>
+
+#include <complex>
+
+namespace everybeam {
+namespace circularsymmetric {
+//! Holds the information for a symmetric voltage pattern
+class VoltagePattern {
+ public:
+ VoltagePattern(double frequency, double maximum_radius_arc_min) {
+ frequencies_.assign(1, frequency);
+ maximum_radius_arc_min_ = maximum_radius_arc_min;
+ };
+
+ size_t NSamples() const { return values_.size() / frequencies_.size(); }
+
+ const double* FreqIndexValues(size_t freq_index) const {
+ return &values_[freq_index * NSamples()];
+ }
+
+ void EvaluatePolynomial(const aocommon::UVector<double>& coefficients,
+ bool invert);
+
+ // void Render(class PrimaryBeamImageSet& beamImages, double pixel_scale_x,
+ // double pixel_scale_y, double phase_centre_ra, double
+ // phase_centre_dec, double pointing_ra, double pointing_dec,
+ // double phase_centre_dl, double phase_centre_dm, double
+ // frequency_hz) const;
+
+ void Render(std::complex<float>* aterm, size_t width, size_t height,
+ double pixel_scale_x, double pixel_scale_y,
+ double phase_centre_ra, double phase_centre_dec,
+ double pointing_ra, double pointing_dec, double phase_centre_dl,
+ double phase_centre_dm, double frequency_hz) const;
+
+ private:
+ // Only works when frequencies_.size() > 1
+ aocommon::UVector<double> InterpolateValues(double freq) const;
+ // Works for any frequencies_.size(), including when 1
+ const double* InterpolateValues(
+ double frequency_hz,
+ aocommon::UVector<double>& interpolated_values) const;
+
+ double LmMaxSquared(double frequency_hz) const;
+
+ double inverse_increment_radius_, maximum_radius_arc_min_;
+
+ // These are the radial (one-dimensional) values of the beam
+ // It is an array of size nsamples x nfrequencies, where the sample index is
+ // least significant (fastest changing)
+ aocommon::UVector<double> values_;
+
+ // Array of since nfrequencies
+ aocommon::UVector<double> frequencies_;
+};
+} // namespace circularsymmetric
+} // namespace everybeam
+#endif // EVERYBEAM_CIRCULARSYMMETRIC_VOLTAGE_PATTERN_H_
diff --git a/cpp/griddedresponse/CMakeLists.txt b/cpp/griddedresponse/CMakeLists.txt
index f2ec3de77b039950dfdec503945b130629bc156b..300981f2e38b7342f5b910fd834b843d95c4f028 100644
--- a/cpp/griddedresponse/CMakeLists.txt
+++ b/cpp/griddedresponse/CMakeLists.txt
@@ -1,4 +1,5 @@
install (FILES
griddedresponse.h
lofargrid.h
+ dishgrid.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/griddedresponse")
diff --git a/cpp/griddedresponse/dishgrid.cc b/cpp/griddedresponse/dishgrid.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d340730d46ef2669c993048468acd8bb5cce3c27
--- /dev/null
+++ b/cpp/griddedresponse/dishgrid.cc
@@ -0,0 +1,40 @@
+#include "dishgrid.h"
+#include "../telescope/dish.h"
+#include "../circularsymmetric/voltagepattern.h"
+#include "../circularsymmetric/vlabeam.h"
+
+#include <aocommon/uvector.h>
+
+#include <algorithm>
+
+using namespace everybeam;
+using namespace everybeam::griddedresponse;
+
+void DishGrid::CalculateStation(std::complex<float>* buffer, double,
+ double frequency, const size_t,
+ const size_t field_id) {
+ const telescope::Dish& dishtelescope =
+ static_cast<const telescope::Dish&>(*telescope_);
+
+ double pdir_ra = dishtelescope.ms_properties_.field_pointing[field_id].first,
+ pdir_dec =
+ dishtelescope.ms_properties_.field_pointing[field_id].second;
+ std::array<double, 5> coefs =
+ circularsymmetric::VLABeam::GetCoefficients("", frequency);
+ circularsymmetric::VoltagePattern vp(frequency, 53.0);
+ aocommon::UVector<double> coefs_vec(coefs.begin(), coefs.end());
+ vp.EvaluatePolynomial(coefs_vec, false);
+ vp.Render(buffer, width_, height_, dl_, dm_, ra_, dec_, pdir_ra, pdir_dec,
+ phase_centre_dl_, phase_centre_dm_, frequency);
+};
+
+void DishGrid::CalculateAllStations(std::complex<float>* buffer, double,
+ double frequency, const size_t field_id) {
+ CalculateStation(buffer, 0., frequency, 0, field_id);
+
+ // Just repeat nstations times
+ for (size_t i = 1; i != telescope_->GetNrStations(); ++i) {
+ std::copy_n(buffer, width_ * height_ * 4,
+ buffer + i * width_ * height_ * 4);
+ }
+}
diff --git a/cpp/griddedresponse/dishgrid.h b/cpp/griddedresponse/dishgrid.h
new file mode 100644
index 0000000000000000000000000000000000000000..8df4f45c87cacec07eae2d1a4e6baab28aa049cd
--- /dev/null
+++ b/cpp/griddedresponse/dishgrid.h
@@ -0,0 +1,47 @@
+// vlagrid.h: Class for computing the VLA 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_DISHGRID_H_
+#define EVERYBEAM_GRIDDEDRESPONSE_DISHGRID_H_
+
+#include "griddedresponse.h"
+
+namespace everybeam {
+namespace griddedresponse {
+
+class DishGrid final : public GriddedResponse {
+ public:
+ DishGrid(telescope::Telescope* telescope_ptr,
+ const coords::CoordinateSystem coordinate_system)
+ : 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;
+
+ void CalculateAllStations(std::complex<float>* buffer, double time,
+ double frequency, const size_t field_id) override;
+};
+} // namespace griddedresponse
+} // namespace everybeam
+#endif // EVERYBEAM_GRIDDEDRESPONSE_DISHGRID_H_
\ No newline at end of file
diff --git a/cpp/griddedresponse/griddedresponse.h b/cpp/griddedresponse/griddedresponse.h
index 67145981b86cd1486cfbc0d81daaa4ec583f47c5..0f97cd3766a1cbd87ae2a1075388da57f4ea5abe 100644
--- a/cpp/griddedresponse/griddedresponse.h
+++ b/cpp/griddedresponse/griddedresponse.h
@@ -59,7 +59,8 @@ class GriddedResponse {
* @param frequency Frequency (Hz)
*/
virtual void CalculateStation(std::complex<float>* buffer, double time,
- double freq, const size_t station_idx) = 0;
+ double freq, const size_t station_idx,
+ const size_t field_id = 0) = 0;
/**
* @brief Compute the Beam response for all stations in a Telescope
@@ -70,7 +71,8 @@ class GriddedResponse {
* @param frequency Frequency (Hz)
*/
virtual void CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency) = 0;
+ double frequency,
+ const size_t field_id = 0) = 0;
std::size_t GetBufferSize(std::size_t nstations) {
return std::size_t(nstations * width_ * height_ * 2 * 2);
diff --git a/cpp/griddedresponse/lofargrid.cc b/cpp/griddedresponse/lofargrid.cc
index 466c18eef8178e60eee9ce425faf40c0e9f65185..03abdc5c25cc4ce9d02fabd56b5481d4c824656c 100644
--- a/cpp/griddedresponse/lofargrid.cc
+++ b/cpp/griddedresponse/lofargrid.cc
@@ -26,7 +26,10 @@ LOFARGrid::LOFARGrid(telescope::Telescope* telescope_ptr,
};
void LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
- double frequency, const size_t station_idx) {
+ double frequency, const size_t station_idx,
+ const size_t) {
+ const telescope::LOFAR& lofartelescope =
+ static_cast<const telescope::LOFAR&>(*telescope_);
aocommon::Lane<Job> lane(nthreads_);
lane_ = &lane;
@@ -34,7 +37,7 @@ void LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
double sb_freq = use_channel_frequency_ ? frequency : subband_frequency_;
// Dummy calculation of gain matrix, needed for multi-threading
inverse_central_gain_.resize(1);
- matrix22c_t gain_matrix = telescope_->GetStation(station_idx)
+ matrix22c_t gain_matrix = lofartelescope.GetStation(station_idx)
->Response(time, frequency, diff_beam_centre_,
sb_freq, station0_, tile0_);
@@ -60,7 +63,9 @@ void LOFARGrid::CalculateStation(std::complex<float>* buffer, double time,
}
void LOFARGrid::CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency) {
+ double frequency, const size_t) {
+ const telescope::LOFAR& lofartelescope =
+ static_cast<const telescope::LOFAR&>(*telescope_);
aocommon::Lane<Job> lane(nthreads_);
lane_ = &lane;
@@ -69,11 +74,11 @@ void LOFARGrid::CalculateAllStations(std::complex<float>* buffer, double time,
double sb_freq = use_channel_frequency_ ? frequency : subband_frequency_;
// Dummy loop, needed for multi-threading
- inverse_central_gain_.resize(telescope_->GetNrStations());
- for (size_t i = 0; i != telescope_->GetNrStations(); ++i) {
- matrix22c_t gain_matrix = telescope_->GetStation(i)->Response(
+ inverse_central_gain_.resize(lofartelescope.GetNrStations());
+ for (size_t i = 0; i != lofartelescope.GetNrStations(); ++i) {
+ matrix22c_t gain_matrix = lofartelescope.GetStation(i)->Response(
time, frequency, diff_beam_centre_, sb_freq, station0_, tile0_);
- if (telescope_->GetOptions().use_differential_beam) {
+ if (lofartelescope.GetOptions().use_differential_beam) {
inverse_central_gain_[i][0] = gain_matrix[0][0];
inverse_central_gain_[i][1] = gain_matrix[0][1];
inverse_central_gain_[i][2] = gain_matrix[1][0];
@@ -91,7 +96,7 @@ void LOFARGrid::CalculateAllStations(std::complex<float>* buffer, double time,
}
for (size_t y = 0; y != height_; ++y) {
- for (size_t antenna_idx = 0; antenna_idx != telescope_->GetNrStations();
+ for (size_t antenna_idx = 0; antenna_idx != lofartelescope.GetNrStations();
++antenna_idx) {
lane.write(Job{
.y = y, .antenna_idx = antenna_idx, .buffer_offset = antenna_idx});
@@ -133,6 +138,8 @@ void LOFARGrid::SetITRFVectors(double time) {
void LOFARGrid::CalcThread(std::complex<float>* buffer, double time,
double frequency) {
+ const telescope::LOFAR& lofartelescope =
+ static_cast<const telescope::LOFAR&>(*telescope_);
const size_t values_per_ant = width_ * height_ * 4;
double sb_freq = use_channel_frequency_ ? frequency : subband_frequency_;
@@ -160,11 +167,11 @@ void LOFARGrid::CalcThread(std::complex<float>* buffer, double time,
std::complex<float>* ant_buffer_ptr =
base_buffer + job.buffer_offset * values_per_ant;
- matrix22c_t gain_matrix = telescope_->GetStation(job.antenna_idx)
+ matrix22c_t gain_matrix = lofartelescope.GetStation(job.antenna_idx)
->Response(time, frequency, itrf_direction,
sb_freq, station0_, tile0_);
- if (telescope_->GetOptions().use_differential_beam) {
+ if (lofartelescope.GetOptions().use_differential_beam) {
aocommon::MC2x2F station_gains;
station_gains[0] = gain_matrix[0][0];
station_gains[1] = gain_matrix[0][1];
diff --git a/cpp/griddedresponse/lofargrid.h b/cpp/griddedresponse/lofargrid.h
index 1ad8795bfc86f9338aed48449edbb5467f1f470e..e8d8f94c6d419a8986668a9e6d23c4c5c9ec19b8 100644
--- a/cpp/griddedresponse/lofargrid.h
+++ b/cpp/griddedresponse/lofargrid.h
@@ -61,7 +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) override;
+ double frequency, const size_t station_idx,
+ const size_t) override;
/**
* @brief Compute the Beam response for all stations in a Telescope
@@ -72,7 +73,7 @@ class LOFARGrid final : public GriddedResponse {
* @param freq Frequency (Hz)
*/
void CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency) override;
+ double frequency, const size_t) override;
private:
casacore::MDirection delay_dir_, tile_beam_dir_, preapplied_beam_dir_;
diff --git a/cpp/load.cc b/cpp/load.cc
index 5f8a5a0d5ba7a40658cef269057ad8c1106c8d6a..a91810d40063bf90370396652bc4e5853bff9a06 100644
--- a/cpp/load.cc
+++ b/cpp/load.cc
@@ -4,22 +4,33 @@
#include <casacore/ms/MeasurementSets/MeasurementSet.h>
#include <casacore/tables/Tables/ScalarColumn.h>
+#include <stdexcept>
+
using namespace everybeam;
namespace everybeam {
namespace {
-enum TelescopeType { kUnknownTelescope, kLofarTelescope };
+enum TelescopeType {
+ kUnknownTelescope,
+ kLofarTelescope,
+ kVLATelescope,
+ kATCATelescope
+};
-TelescopeType Convert(const std::string &str) {
- if (str == "LOFAR")
- return kLofarTelescope;
+TelescopeType Convert(const std::string &telescope_name) {
+ if (telescope_name == "LOFAR") return kLofarTelescope;
+ // Maybe more elegant with boost::to_upper_copy()?
+ else if (telescope_name.compare(0, 4, "EVLA") == 0)
+ return kVLATelescope;
+ else if (telescope_name.compare(0, 4, "ATCA") == 0)
+ return kATCATelescope;
else
return kUnknownTelescope;
}
} // namespace
std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
- const ElementResponseModel model,
- const Options &options) {
+ const Options &options,
+ const ElementResponseModel model) {
// Read Telescope name and convert to enum
casacore::ScalarColumn<casacore::String> telescope_name_col(ms.observation(),
"TELESCOPE_NAME");
@@ -32,6 +43,18 @@ std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
new telescope::LOFAR(ms, model, options));
return telescope;
}
+ case kATCATelescope: {
+ std::unique_ptr<telescope::Telescope> telescope =
+ std::unique_ptr<telescope::Telescope>(
+ new telescope::Dish(ms, options));
+ return telescope;
+ }
+ case kVLATelescope: {
+ std::unique_ptr<telescope::Telescope> telescope =
+ std::unique_ptr<telescope::Telescope>(
+ new telescope::Dish(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 19102da3d9e32b5427015f463b5f66b5cf595e6d..ff9ed40467e3d51eb07de2c766abc4f18777b8d9 100644
--- a/cpp/load.h
+++ b/cpp/load.h
@@ -26,6 +26,7 @@
#include "./telescope/telescope.h"
#include "./telescope/lofar.h"
+#include "./telescope/dish.h"
#include "options.h"
namespace everybeam {
@@ -40,7 +41,8 @@ namespace everybeam {
* @return telescope::Telescope::Ptr
*/
std::unique_ptr<telescope::Telescope> Load(
- casacore::MeasurementSet &ms, const ElementResponseModel model,
- const Options &options = Options::GetDefault());
+ casacore::MeasurementSet &ms,
+ const Options &options = Options::GetDefault(),
+ const ElementResponseModel model = ElementResponseModel::kHamaker);
} // namespace everybeam
#endif // EVERYBEAM_LOAD_H_
\ No newline at end of file
diff --git a/cpp/telescope/CMakeLists.txt b/cpp/telescope/CMakeLists.txt
index 27539ca1783f0e1849892e1634e7a351ee513ff5..23934e1b6a917fae0ae4f2d63a9d2dbf43d6644d 100644
--- a/cpp/telescope/CMakeLists.txt
+++ b/cpp/telescope/CMakeLists.txt
@@ -1,4 +1,5 @@
install (FILES
+ dish.h
lofar.h
telescope.h
DESTINATION "include/${CMAKE_PROJECT_NAME}/telescope")
diff --git a/cpp/telescope/atcatelescope.h b/cpp/telescope/atcatelescope.h
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/cpp/telescope/dish.cc b/cpp/telescope/dish.cc
new file mode 100644
index 0000000000000000000000000000000000000000..0aeba4e28b53c1cdc9eaf1c685ce68573662ed58
--- /dev/null
+++ b/cpp/telescope/dish.cc
@@ -0,0 +1,29 @@
+#include "dish.h"
+#include "../griddedresponse/dishgrid.h"
+
+#include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
+
+using namespace everybeam;
+using namespace everybeam::telescope;
+
+Dish::Dish(casacore::MeasurementSet &ms, const Options &options)
+ : Telescope(ms, options) {
+ casacore::MSField field_table = ms.field();
+ casacore::ArrayColumn<double> pointing_dir_col(
+ field_table, casacore::MSField::columnName(casacore::MSField::DELAY_DIR));
+
+ for (std::size_t field_id = 0; field_id != field_table.nrow(); ++field_id) {
+ casacore::Array<double> pdir = pointing_dir_col(field_id);
+ double pdir_ra = *pdir.cbegin();
+ double pdir_dec = *(pdir.cbegin() + 1);
+ ms_properties_.field_pointing.emplace_back(pdir_ra, pdir_dec);
+ }
+}
+
+std::unique_ptr<griddedresponse::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));
+ return grid;
+};
\ No newline at end of file
diff --git a/cpp/telescope/dish.h b/cpp/telescope/dish.h
new file mode 100644
index 0000000000000000000000000000000000000000..efa2a0b9ffcf03e7af7b4a2df9d144cabe1059c4
--- /dev/null
+++ b/cpp/telescope/dish.h
@@ -0,0 +1,63 @@
+// dishtelescope.h: Base class for dish telescopes (VLA, ATCA, ...).
+// 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_DISH_H_
+#define EVERYBEAM_TELESCOPE_DISH_H_
+
+#include "telescope.h"
+
+#include <casacore/measures/Measures/MDirection.h>
+
+namespace everybeam {
+
+namespace griddedresponse {
+class DishGrid;
+class GriddedResponse;
+} // namespace griddedresponse
+
+namespace telescope {
+
+/**
+ * This class calculates the a-terms for dishes with a circularly symmetric
+ * response.
+ */
+class Dish : public Telescope {
+ friend class griddedresponse::DishGrid;
+
+ public:
+ Dish(casacore::MeasurementSet &ms, const Options &options);
+
+ std::unique_ptr<griddedresponse::GriddedResponse> GetGriddedResponse(
+ const coords::CoordinateSystem &coordinate_system) override;
+
+ private:
+ struct MSProperties {
+ std::vector<std::pair<double, double>> field_pointing;
+ };
+ MSProperties ms_properties_;
+};
+} // namespace telescope
+} // namespace everybeam
+
+#endif // EVERYBEAM_TELESCOPE_DISH_H_
\ No newline at end of file
diff --git a/cpp/telescope/lofar.cc b/cpp/telescope/lofar.cc
index 49c800a355ee22edf669d4a48c7a22b2d69a8426..263b5e01cd9b7318edaaf8437f43729f010f7561 100644
--- a/cpp/telescope/lofar.cc
+++ b/cpp/telescope/lofar.cc
@@ -58,7 +58,8 @@ bool GetPreappliedBeamDirection(casacore::MeasurementSet &ms,
LOFAR::LOFAR(MeasurementSet &ms, const ElementResponseModel model,
const Options &options)
- : Telescope(ms, model, options) {
+ : Telescope(ms, options) {
+ stations_.resize(nstations_);
ReadAllStations(ms, model);
// Populate MeasurementSet properties struct
@@ -102,7 +103,6 @@ std::unique_ptr<griddedresponse::GriddedResponse> LOFAR::GetGriddedResponse(
// Get and return GriddedResponse ptr
std::unique_ptr<griddedresponse::GriddedResponse> grid(
new griddedresponse::LOFARGrid(this, coordinate_system));
- // griddedresponse::GriddedResponse grid(LOFARGrid(this, coordinate_system));
return grid;
};
diff --git a/cpp/telescope/lofar.h b/cpp/telescope/lofar.h
index 553116e078993b5d962416594e814c27b9c7959a..6e52d958eddd8661487630a606ac60253fd03c2b 100644
--- a/cpp/telescope/lofar.h
+++ b/cpp/telescope/lofar.h
@@ -25,6 +25,8 @@
#ifndef EVERYBEAM_TELESCOPE_LOFAR_H_
#define EVERYBEAM_TELESCOPE_LOFAR_H_
+#include "../station.h"
+#include "../elementresponse.h"
#include "telescope.h"
#include <casacore/measures/Measures/MPosition.h>
@@ -59,14 +61,44 @@ class LOFAR final : public Telescope {
std::unique_ptr<griddedresponse::GriddedResponse> GetGriddedResponse(
const coords::CoordinateSystem &coordinate_system) override;
+ /**
+ * @brief Get station by index
+ *
+ * @param station_id Station index to retrieve
+ * @return Station::Ptr
+ */
+ Station::Ptr GetStation(std::size_t station_idx) const {
+ // Assert only in DEBUG mode
+ assert(station_idx < nstations_);
+ return stations_[station_idx];
+ }
+
private:
+ /**
+ * @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) {
+ stations_[i] = ReadStation(ms, i, model);
+ }
+ };
+
Station::Ptr ReadStation(const casacore::MeasurementSet &ms,
const std::size_t id,
- const ElementResponseModel model) const override;
+ const ElementResponseModel model) const;
struct MSProperties {
double subband_freq;
casacore::MDirection delay_dir, tile_beam_dir, preapplied_beam_dir;
};
+
+ std::vector<Station::Ptr> stations_;
MSProperties ms_properties_;
};
} // namespace telescope
diff --git a/cpp/telescope/telescope.h b/cpp/telescope/telescope.h
index 979017313c5b753f31c9f7c03d198be71c891713..f0e58671b50b31c7d96fba6589c593a98e5b61ed 100644
--- a/cpp/telescope/telescope.h
+++ b/cpp/telescope/telescope.h
@@ -24,9 +24,7 @@
#ifndef EVERYBEAM_TELESCOPE_TELESCOPE_H_
#define EVERYBEAM_TELESCOPE_TELESCOPE_H_
-#include "../station.h"
#include "../options.h"
-#include "../elementresponse.h"
#include <vector>
#include <memory>
@@ -61,18 +59,6 @@ class Telescope {
virtual std::unique_ptr<griddedresponse::GriddedResponse> GetGriddedResponse(
const coords::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_idx) const {
- // Assert only in DEBUG mode
- assert(station_idx < nstations_);
- return stations_[station_idx];
- }
-
std::size_t GetNrStations() const { return nstations_; };
Options GetOptions() const { return options_; };
@@ -82,40 +68,13 @@ class Telescope {
* @brief Construct a new Telescope object
*
* @param ms MeasurementSet
- * @param model ElementResponse model
* @param options telescope options
*/
- Telescope(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) {
- stations_[i] = ReadStation(ms, i, model);
- }
- };
+ Telescope(casacore::MeasurementSet &ms, const Options &options)
+ : nstations_(ms.antenna().nrow()), options_(options){};
std::size_t nstations_;
Options options_;
- std::vector<Station::Ptr> stations_;
-
- 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
diff --git a/cpp/telescope/vlatelescope.h b/cpp/telescope/vlatelescope.h
new file mode 100644
index 0000000000000000000000000000000000000000..854d9a143f303e9eb12d06a5a98d0e7924fb9d98
--- /dev/null
+++ b/cpp/telescope/vlatelescope.h
@@ -0,0 +1 @@
+#ifndef
diff --git a/cpp/test/CMakeLists.txt b/cpp/test/CMakeLists.txt
index bc662495c34c726ee1cecf36c84be7b0443857b6..52cc4287b5bedc147e7f67b170988329af5f0a77 100644
--- a/cpp/test/CMakeLists.txt
+++ b/cpp/test/CMakeLists.txt
@@ -1,4 +1,5 @@
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
@@ -20,13 +21,21 @@ if(Boost_FOUND)
include_directories(${Boost_INCLUDE_DIR})
set(TEST_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/runtests.cc)
- # Test loading of LOFAR
+ # Test LOFAR beamforming
if( LOFAR_MOCK_MS )
set(TEST_SOURCES
${TEST_SOURCES}
${CMAKE_CURRENT_SOURCE_DIR}/tload_lofar.cc
)
- endif()
+ endif()
+
+ # Test VLA beamforming
+ if( VLA_MOCK_MS )
+ set(TEST_SOURCES
+ ${TEST_SOURCES}
+ ${CMAKE_CURRENT_SOURCE_DIR}/tvla.cc
+ )
+ endif()
# Add test executable
add_executable(runtests ${TEST_SOURCES})
diff --git a/cpp/test/tload_lofar.cc b/cpp/test/tload_lofar.cc
index c061eaa13c5a3cbcb552c0b3fd19ba060a218280..4958df8defcfaa453879ed2fe3f47356104fb7ab 100644
--- a/cpp/test/tload_lofar.cc
+++ b/cpp/test/tload_lofar.cc
@@ -20,7 +20,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
// Load LOFAR Telescope
std::unique_ptr<telescope::Telescope> telescope =
- Load(ms, response_model, options);
+ Load(ms, options, response_model);
// Assert if we indeed have a LOFAR pointer
BOOST_CHECK(nullptr != dynamic_cast<telescope::LOFAR*>(telescope.get()));
@@ -30,7 +30,9 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
BOOST_CHECK_EQUAL(telescope->GetNrStations(), nstations);
// Assert if GetStation(stationd_id) behaves properly
- BOOST_CHECK_EQUAL(telescope->GetStation(0)->GetName(), "CS001HBA0");
+ const telescope::LOFAR& lofartelescope =
+ static_cast<const telescope::LOFAR&>(*telescope.get());
+ BOOST_CHECK_EQUAL(lofartelescope.GetStation(0)->GetName(), "CS001HBA0");
// Properties extracted from MS
double time = 4929192878.008341;
@@ -86,13 +88,13 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
1e-4);
}
- // std::vector<std::complex<float>> antenna_buffer_all(
- // grid_response->GetBufferSize(telescope->GetNrStations()));
- // grid_response->CalculateAllStations(antenna_buffer_all.data(), time,
- // frequency);
- // BOOST_CHECK_EQUAL(
- // antenna_buffer_all.size(),
- // std::size_t(telescope->GetNrStations() * width * height * 2 * 2));
+ std::vector<std::complex<float>> antenna_buffer_all(
+ grid_response->GetBufferSize(telescope->GetNrStations()));
+ grid_response->CalculateAllStations(antenna_buffer_all.data(), time,
+ frequency);
+ BOOST_CHECK_EQUAL(
+ antenna_buffer_all.size(),
+ std::size_t(telescope->GetNrStations() * width * height * 2 * 2));
// Test with differential beam, single
Options options_diff_beam;
@@ -100,7 +102,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
// Load LOFAR Telescope
std::unique_ptr<telescope::Telescope> telescope_diff_beam =
- Load(ms, response_model, options_diff_beam);
+ Load(ms, options_diff_beam, response_model);
std::unique_ptr<griddedresponse::GriddedResponse> grid_response_diff_beam =
telescope_diff_beam->GetGriddedResponse(coord_system);
@@ -115,9 +117,4 @@ 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("station_responses.npy", false, 4, leshape,
- // antenna_buffer_single);
}
\ No newline at end of file
diff --git a/cpp/test/tvla.cc b/cpp/test/tvla.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4a4da0bc8de7e778ee8d7598099f6916c0e7be9d
--- /dev/null
+++ b/cpp/test/tvla.cc
@@ -0,0 +1,91 @@
+#include <boost/test/unit_test.hpp>
+
+#include "../load.h"
+#include "../options.h"
+#include "../griddedresponse/dishgrid.h"
+#include "../elementresponse.h"
+#include "../../external/npy.hpp"
+
+#include "config.h"
+#include <complex>
+#include <cmath>
+
+using namespace everybeam;
+
+BOOST_AUTO_TEST_CASE(load_vla) {
+ casacore::MeasurementSet ms(VLA_MOCK_MS);
+
+ std::unique_ptr<telescope::Telescope> telescope = Load(ms);
+
+ // Assert if we indeed have a LOFAR 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);
+
+ double time = 0.5 * (4.90683119e+09 + 4.90684196e+09);
+ double frequency = 0.5e+09; // 1991000000.0;
+ std::size_t width(16), height(16);
+ double ra(2.62880729), dec(0.02831797), dl(0.125 * M_PI / 180.),
+ dm(0.125 * 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::DishGrid*>(grid_response.get()));
+
+ std::vector<std::complex<float>> antenna_buffer(
+ grid_response->GetBufferSize(telescope->GetNrStations()));
+ grid_response->CalculateAllStations(antenna_buffer.data(), time, frequency);
+
+ // Check that XX/YY are equal
+ // Loop over pixels to check that off-diagonals are 0
+ // and diagonal entries are equal
+ for (std::size_t pxl = 0; pxl < width * height; ++pxl) {
+ BOOST_CHECK(std::abs(antenna_buffer[pxl * 4 + 1]) < 1e-8);
+ BOOST_CHECK(std::abs(antenna_buffer[pxl * 4 + 2]) < 1e-8);
+ BOOST_CHECK(
+ std::abs(antenna_buffer[pxl * 4] - antenna_buffer[pxl * 4 + 3]) < 1e-8);
+ }
+
+ // Check whether pixels are reproduced correctly at certain pixels on 16x16
+ // image VLABeam output at pixel (0, 0):
+ std::vector<std::complex<float>> vla_p00 = {
+ {0.4687362, 0.}, {0, 0}, {0, 0}, {0.4687362, 0.}};
+ // VLABeam output at pixel (2, 3):
+ std::vector<std::complex<float>> vla_p23 = {
+ {0.65400606, 0.}, {0, 0}, {0, 0}, {0.65400606, 0.}};
+ // VLABeam output at pixel (10, 12):
+ 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;
+
+ // Check if results are reproduced
+ BOOST_CHECK_EQUAL_COLLECTIONS(antenna_buffer.begin() + offset_00,
+ antenna_buffer.begin() + offset_00 + 4,
+ vla_p00.begin(), vla_p00.end());
+ BOOST_CHECK_EQUAL_COLLECTIONS(antenna_buffer.begin() + offset_23,
+ antenna_buffer.begin() + offset_23 + 4,
+ vla_p23.begin(), vla_p23.end());
+ BOOST_CHECK_EQUAL_COLLECTIONS(antenna_buffer.begin() + offset_1012,
+ antenna_buffer.begin() + offset_1012 + 4,
+ vla_p1012.begin(), vla_p1012.end());
+
+ // Print to np array
+ const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
+ npy::SaveArrayAsNumpy("vla_station_responses.npy", false, 4, leshape,
+ antenna_buffer);
+}
\ No newline at end of file
diff --git a/scripts/download_vla_ms.sh b/scripts/download_vla_ms.sh
new file mode 100755
index 0000000000000000000000000000000000000000..7c18ca8ed70e25d788240fb0b9cddd7337bcef51
--- /dev/null
+++ b/scripts/download_vla_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/
+
+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
+fi
+
+if [ -d $VLA_MOCK_MS ]
+then
+ echo "Directory already exists"
+else
+ mkdir $VLA_MOCK_MS
+fi
+
+tar -xf $VLA_MOCK_ARCHIVE -C $VLA_MOCK_MS --strip-components=1
+rm $VLA_MOCK_ARCHIVE
\ No newline at end of file