Merge branch '46-add-stationresponse-comparison-for-oskar' into 'master'

Resolve "Add stationresponse comparison for oskar"

Closes #46

See merge request !59

.gitlab-ci.yml

@@ -78,9 +78,10 @@ build-compare-oskar:
   before_script:
     - apt-get -y install python3-pip
     # Install python requirements for the OSKAR "integration" test
-    - pip3 install numpy==1.19.0 scipy h5py astropy tqdm matplotlib
+    - pip3 install numpy==1.19.0 scipy h5py astropy tqdm matplotlib pandas lofarantpos
     - mkdir -p /opt/oskar/build
-    - cd /opt/oskar && git clone https://github.com/OxfordSKA/OSKAR.git
+    - cd /opt/oskar && git clone https://github.com/basvandertol/OSKAR.git
+    - cd OSKAR && git checkout fix-transformations
     - cd /opt/oskar/build
   script:
     # OSKAR cpp install
@@ -94,7 +95,10 @@ build-compare-oskar:
     # Run OSKAR comparison, set some env variables for this session
     - export NPIXELS=8 && export APPLY_TRANSPOSE=OFF && MAX_ORDER=3 && TOLERANCE=1e-12
     - cd /opt/everybeam/build
-    - make VERBOSE=1 comparison-oskar
+    - make VERBOSE=1 comparison-oskar-basefunctions
+    # Run OSKAR stationresponse comparison
+    - export NPIXELS=32 TOLERANCE=1e-5
+    - make VERBOSE=1 comparison-oskar-station-response
 
 deploy-doc:
   stage: deploy

cpp/CMakeLists.txt

@@ -19,10 +19,12 @@ add_library(everybeam SHARED
   coords/itrfconverter.cc
   coords/itrfdirection.cc
   lofarreadutils.cc
+  msv3readutils.cc
   station.cc
   telescope/lofar.cc
   telescope/dish.cc
   telescope/mwa.cc
+  telescope/oskar.cc
   griddedresponse/lofargrid.cc
   griddedresponse/dishgrid.cc
   griddedresponse/mwagrid.cc
@@ -33,6 +35,12 @@ add_library(everybeam SHARED
   mwabeam/beam2016implementation.cc
 )
 
+# Make sure that when other targets within this project link against the everybeam target,
+# they can find the include files.
+target_include_directories(everybeam PUBLIC
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+)
+
 target_include_directories(everybeam PUBLIC ${CASACORE_INCLUDE_DIR})
 target_link_libraries(everybeam PUBLIC hamaker lobes oskar)
 target_link_libraries(everybeam PUBLIC ${CASACORE_LIBRARIES} ${HDF5_LIBRARIES})
@@ -49,6 +57,7 @@ install (FILES
   element.h
   elementresponse.h
   lofarreadutils.h
+  msv3readutils.h
   station.h
   # Related to new API:
   load.h

cpp/load.cc

@@ -25,6 +25,9 @@ TelescopeType GetTelescopeType(const casacore::MeasurementSet &ms) {
     return kATCATelescope;
   else if (telescope_name == "MWA")
     return kMWATelescope;
+  // check if telescope_name starts with "OSKAR"
+  else if (telescope_name.rfind("OSKAR", 0) == 0)
+    return kOSKARTelescope;
   else
     return kUnknownTelescope;
 }
@@ -58,6 +61,12 @@ std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
               new telescope::MWA(ms, options));
       return telescope;
     }
+    case kOSKARTelescope: {
+      std::unique_ptr<telescope::Telescope> telescope =
+          std::unique_ptr<telescope::Telescope>(
+              new telescope::OSKAR(ms, options));
+      return telescope;
+    }
     default:
       casacore::ScalarColumn<casacore::String> telescope_name_col(
           ms.observation(), "TELESCOPE_NAME");

cpp/load.h

@@ -28,6 +28,7 @@
 #include "./telescope/lofar.h"
 #include "./telescope/dish.h"
 #include "./telescope/mwa.h"
+#include "./telescope/oskar.h"
 #include "options.h"
 
 namespace everybeam {
@@ -41,7 +42,8 @@ enum TelescopeType {
   kAARTFAAC,
   kVLATelescope,
   kATCATelescope,
-  kMWATelescope
+  kMWATelescope,
+  kOSKARTelescope
 };
 
 /**

cpp/msv3readutils.cc

+// lofarreadutils.cc: Utility functions to read the meta data relevant for
+// simulating the beam from LOFAR observations stored in MS format.
+//
+// Copyright (C) 2013
+// ASTRON (Netherlands Institute for Radio Astronomy)
+// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//
+// This file is part of the LOFAR software suite.
+// The LOFAR software suite is free software: you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// The LOFAR software suite is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//
+// $Id$
+
+#include "lofarreadutils.h"
+#include "beamformeridenticalantennas.h"
+#include "common/mathutils.h"
+#include "common/casautils.h"
+
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MCDirection.h>
+#include <casacore/measures/Measures/MCPosition.h>
+#include <casacore/measures/Measures/MeasTable.h>
+#include <casacore/measures/Measures/MeasConvert.h>
+#include <casacore/measures/TableMeasures/ScalarMeasColumn.h>
+
+#include <cassert>
+#include <stdexcept>
+
+#include <casacore/ms/MeasurementSets/MSAntenna.h>
+#include <casacore/ms/MSSel/MSSelection.h>
+#include <casacore/ms/MSSel/MSAntennaParse.h>
+#include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
+#include <casacore/ms/MeasurementSets/MSDataDescription.h>
+#include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
+#include <casacore/ms/MeasurementSets/MSField.h>
+#include <casacore/ms/MeasurementSets/MSFieldColumns.h>
+#include <casacore/ms/MeasurementSets/MSObservation.h>
+#include <casacore/ms/MeasurementSets/MSObsColumns.h>
+#include <casacore/ms/MeasurementSets/MSPolarization.h>
+#include <casacore/ms/MeasurementSets/MSPolColumns.h>
+#include <casacore/ms/MeasurementSets/MSSpectralWindow.h>
+#include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
+
+namespace everybeam {
+
+constexpr Antenna::CoordinateSystem::Axes antenna_orientation = {
+    {
+        1.0,
+        0.0,
+        0.0,
+    },
+    {
+        0.0,
+        1.0,
+        0.0,
+    },
+    {0.0, 0.0, 1.0},
+};
+
+using namespace casacore;
+
+vector3r_t TransformToFieldCoordinates(
+    const vector3r_t &position, const Antenna::CoordinateSystem::Axes &axes);
+
+BeamFormer::Ptr ReadMSv3AntennaField(const Table &table, unsigned int id,
+                                     ElementResponse::Ptr element_response) {
+  Antenna::CoordinateSystem coordinate_system =
+      common::ReadCoordinateSystem(table, id);
+  BeamFormer::Ptr beam_former(
+      new BeamFormerIdenticalAntennas(coordinate_system));
+
+  ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
+  ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG");
+
+  // Read element offsets and flags.
+  Matrix<Quantity> aips_offset = c_offset(id);
+
+  assert(aips_offset.shape().isEqual(IPosition(2, aips_offset.nrow(), 3)));
+
+  Matrix<Bool> aips_flag = c_flag(id);
+  assert(aips_flag.shape().isEqual(IPosition(2, aips_offset.nrow(), 2)));
+
+  for (size_t i = 0; i < aips_offset.nrow(); ++i) {
+    vector3r_t antenna_position = {aips_offset(i, 0).getValue(),
+                                   aips_offset(i, 1).getValue(),
+                                   aips_offset(i, 2).getValue()};
+    antenna_position =
+        TransformToFieldCoordinates(antenna_position, coordinate_system.axes);
+    Antenna::Ptr antenna;
+    Antenna::CoordinateSystem antenna_coordinate_system{antenna_position,
+                                                        antenna_orientation};
+    antenna = Element::Ptr(
+        new Element(antenna_coordinate_system, element_response, id));
+
+    antenna->enabled_[0] = !aips_flag(i, 0);
+    antenna->enabled_[1] = !aips_flag(i, 1);
+    beam_former->AddAntenna(antenna);
+  }
+  return beam_former;
+}
+
+Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
+                             const ElementResponseModel model) {
+  ROMSAntennaColumns antenna(ms.antenna());
+  assert(antenna.nrow() > id && !antenna.flagRow()(id));
+
+  // Get station name.
+  const string name(antenna.name()(id));
+
+  // Get station position (ITRF).
+  MPosition mPosition =
+      MPosition::Convert(antenna.positionMeas()(id), MPosition::ITRF)();
+  MVPosition mvPosition = mPosition.getValue();
+  const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
+
+  // Create station.
+  Station::Ptr station(new Station(name, position, model));
+
+  Table tab_phased_array = common::GetSubTable(ms, "PHASED_ARRAY");
+
+  // The Station will consist of a BeamFormer that combines the fields
+  // coordinate system is ITRF
+  // phase reference is station position
+  auto beam_former =
+      ReadMSv3AntennaField(tab_phased_array, id, station->GetElementResponse());
+
+  station->SetAntenna(beam_former);
+
+  size_t field_id = 0;
+  size_t element_id = 0;
+  Antenna::CoordinateSystem coordinate_system =
+      common::ReadCoordinateSystem(tab_phased_array, field_id);
+  auto element_response = station->GetElementResponse();
+  auto element = Element::Ptr(
+      new Element(coordinate_system, element_response, element_id));
+  station->SetElement(element);
+
+  return station;
+}
+
+}  // namespace everybeam

cpp/msv3readutils.h

+// msv3readutils.h: Utility functions to read the meta data relevant for
+// simulating the beam from OSKAR simulations stored in MS format.
+//
+// Copyright (C) 2013
+// ASTRON (Netherlands Institute for Radio Astronomy)
+// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//
+// This file is part of the LOFAR software suite.
+// The LOFAR software suite is free software: you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// The LOFAR software suite is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with the LOFAR software suite. If not, see <http://www.gnu.org/licenses/>.
+//
+// $Id$
+
+#ifndef EVERYBEAM_MSV3READUTILS_H_
+#define EVERYBEAM_MSV3READUTILS_H_
+
+// \file
+// Utility functions to read the meta data relevant for simulating the beam from
+// OSKAR simulations stored in MS format.
+
+#include "station.h"
+#include "elementresponse.h"
+
+#include <casacore/ms/MeasurementSets/MeasurementSet.h>
+#include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
+#include <casacore/measures/Measures/MDirection.h>
+
+namespace everybeam {
+const ElementResponseModel defaultElementResponseModel =
+    ElementResponseModel::kUnknown;
+
+/**
+ * @brief Read single station from MeasurementSet
+ *
+ * @param ms Measurement set
+ * @param id Station id
+ * @param model Element response model
+ * @return Station::Ptr
+ */
+Station::Ptr ReadMSv3Station(
+    const casacore::MeasurementSet &ms, unsigned int id,
+    const ElementResponseModel model = defaultElementResponseModel);
+
+/**
+ * @brief Read multiple stations from measurment set into buffer out_it
+ * Loops over ReadMSv3Station for all the antennas in MeasurementSet
+ *
+ * @tparam T Template type
+ * @param ms Measurement set
+ * @param out_it Out buffer
+ * @param model Element Response buffer
+ */
+template <typename T>
+void ReadMSv3Stations(
+    const casacore::MeasurementSet &ms, T out_it,
+    const ElementResponseModel model = defaultElementResponseModel) {
+  casacore::ROMSAntennaColumns antenna(ms.antenna());
+  for (unsigned int i = 0; i < antenna.nrow(); ++i) {
+    *out_it++ = ReadMSv3Station(ms, i, model);
+  }
+}
+
+}  // namespace everybeam
+#endif  // EVERYBEAM_MSV3READUTILS_H_

cpp/oskar/CMakeLists.txt

@@ -30,3 +30,6 @@ install(
   TARGETS oskar
   EXPORT EveryBeamTargets
   DESTINATION lib)
+
+#install oskar coefficients
+install(FILES "${CMAKE_SOURCE_DIR}/coeffs/oskar.h5" DESTINATION ${CMAKE_INSTALL_DATA_DIR})

cpp/oskar/oskar_evaluate_spherical_wave_sum.cc

@@ -62,7 +62,7 @@ void oskar_evaluate_spherical_wave_sum(int num_points, const FP* theta,
           oskar_legendre2(l, abs_m, cos_t, sin_t, p, pds, dpms);
           if (abs_m == 0) {
             sin_p = (FP)0;
-            cos_p = -sqrt(f_);
+            cos_p = sqrt(f_);
             const FP4c alpha_ = alpha[ind0];
             oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.a, alpha_.b, Xt,
                            Xp);
@@ -74,10 +74,10 @@ void oskar_evaluate_spherical_wave_sum(int num_points, const FP* theta,
             d_fact = std::tgamma(d_ + 1);
             s_fact = std::tgamma(s_ + 1);
             const FP ff = f_ * d_fact / s_fact;
-            const FP nf = sqrt(ff) * (2 * (abs_m & 1) - 1);
+            const FP nf = sqrt(ff);
             const FP4c alpha_m = alpha[ind0 + abs_m];
             const FP4c alpha_p = alpha[ind0 - abs_m];
-            p = -abs_m * phi_x_;
+            p = abs_m * phi_x_;
             oskar_sincos(p, &sin_p, &cos_p);
             sin_p *= nf;
             cos_p *= nf;
@@ -86,7 +86,7 @@ void oskar_evaluate_spherical_wave_sum(int num_points, const FP* theta,
             sin_p = -sin_p;
             oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.a, alpha_p.b,
                            Xt, Xp);
-            p = -abs_m * phi_y_;
+            p = abs_m * phi_y_;
             oskar_sincos(p, &sin_p, &cos_p);
             sin_p *= nf;
             cos_p *= nf;

cpp/oskar/oskarelementresponse.cc

@@ -51,7 +51,15 @@ void OSKARElementResponseSphericalWave::Response(
   std::complex<double>* alpha_ptr = dataset->GetAlphaPtr(element_id);
 
   double phi_x = phi;
-  double phi_y = phi + M_PI_2;
+  double phi_y = phi;
+
+  // TODO: phi_x and phi_y can have different values if there is only one set
+  // of coefficients that is is used for both dipoles.
+  // In that case it is assumed the Y dipole rotated 90deg with respect
+  // to the X dipole, so then phi_y = phi+ M_PI_2.
+  // That case needs to be detected when the coefficients are read,
+  // and here phi_y needs to be set accordingly.
+
   oskar_evaluate_spherical_wave_sum_double(1, &theta, &phi_x, &phi_y, l_max,
                                            alpha_ptr, response_ptr);
 }

cpp/station.cc

@@ -64,6 +64,10 @@ void Station::SetResponseModel(const ElementResponseModel model) {
   }
 }
 
+void Station::SetResponse(std::shared_ptr<ElementResponse> element_response) {
+  element_response_.set(element_response);
+}
+
 const std::string &Station::GetName() const { return name_; }
 
 const vector3r_t &Station::GetPosition() const { return position_; }

cpp/station.h

@@ -55,6 +55,8 @@ class Station {
 
   void SetResponseModel(const ElementResponseModel model);
 
+  void SetResponse(std::shared_ptr<ElementResponse> element_response);
+
   //! Return the name of the station.
   const std::string &GetName() const;
 
@@ -333,6 +335,8 @@ class Station {
   //! Set antenna attribute, usually a BeamFormer, but can also be an Element
   void SetAntenna(Antenna::Ptr antenna) { antenna_ = antenna; }
 
+  Antenna::Ptr GetAntenna() { return antenna_; }
+
   //! Set Element attribute
   void SetElement(Element::Ptr element) { element_ = element; }
 

cpp/telescope/CMakeLists.txt

@@ -2,6 +2,8 @@ install (FILES
   dish.h
   lofar.h
   mwa.h 
+  lofar.h 
+  oskar.h
   telescope.h
 DESTINATION "include/${CMAKE_PROJECT_NAME}/telescope")
 

cpp/telescope/oskar.cc

+#include "oskar.h"
+#include "../common/mathutils.h"
+#include "../common/casautils.h"
+#include "../msv3readutils.h"
+
+#include <aocommon/banddata.h>
+#include <cassert>
+#include <casacore/measures/TableMeasures/ArrayMeasColumn.h>
+
+using namespace everybeam;
+using namespace everybeam::telescope;
+using namespace casacore;
+
+OSKAR::OSKAR(MeasurementSet &ms, const Options &options)
+    : Telescope(ms, options) {
+  stations_.resize(nstations_);
+  ReadAllStations(ms, options_.element_response_model);
+}
+
+std::unique_ptr<griddedresponse::GriddedResponse> OSKAR::GetGriddedResponse(
+    const coords::CoordinateSystem &coordinate_system) {
+  throw std::runtime_error(
+      "GetGriddedResponse() is not implemented for OSKAR Telescope");
+
+  // TODO: return an OSKARGrid here, in a similar way to the commented out code
+  // below
+
+  // 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;
+};
+
+Station::Ptr OSKAR::ReadStation(const MeasurementSet &ms, std::size_t id,
+                                const ElementResponseModel model) const {
+  Station::Ptr station = ReadMSv3Station(ms, id, model);
+  return station;
+}

cpp/telescope/oskar.h

+// OSKARTelescope.h: Base class for computing the response for the OSKAR
+// telescope.
+//
+// Copyright (C) 2020
+// ASTRON (Netherlands Institute for Radio Astronomy)
+// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//
+// This file is part of the EveryBeam software suite.
+// The EveryBeam software suite is free software: you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// The EveryBeam software suite is distributed in the hope that it will be
+// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with the EveryBeam software suite. If not, see
+// <http://www.gnu.org/licenses/>.
+//
+// $Id$
+
+#ifndef EVERYBEAM_TELESCOPE_OSKAR_H_
+#define EVERYBEAM_TELESCOPE_OSKAR_H_
+
+#include "../station.h"
+#include "../elementresponse.h"
+#include "telescope.h"
+
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MEpoch.h>
+#include <memory>
+
+namespace everybeam {
+
+namespace telescope {
+
+//! OSKAR telescope class
+class OSKAR final : public Telescope {
+ public:
+  /**
+   * @brief Construct a new OSKAR object
+   *
+   * @param ms MeasurementSet
+   * @param model Element Response model
+   * @param options telescope options
+   */
+  OSKAR(casacore::MeasurementSet &ms, const Options &options);
+
+  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:
+  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;
+
+  std::vector<Station::Ptr> stations_;
+};
+}  // namespace telescope
+}  // namespace everybeam
+
+#endif  // EVERYBEAM_TELESCOPE_OSKAR_H_

demo/comparison-oskar/CMakeLists.txt

 #------------------------------------------------------------------------------
 # CMake file for compiling a comparison between OSKAR and EveryBeam
-add_executable(comparison-oskar-generate-beampattern main.cpp)
+add_executable(make_element_response_image make_element_response_image.cpp)
+target_link_libraries(make_element_response_image oskar)
+
+add_executable(make_station_response_image make_station_response_image.cpp)
+target_link_libraries(make_station_response_image everybeam)
+
+# Required to get the config.h header
+target_include_directories(make_station_response_image PRIVATE "${CMAKE_BINARY_DIR}")
 
-target_link_libraries(comparison-oskar-generate-beampattern oskar)
 
 file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/layout.txt"
      DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/telescope.tm")
@@ -13,15 +19,23 @@ file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/position.txt"
 file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/station000/layout.txt"
      DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/station000/telescope.tm")
 
+execute_process( COMMAND ${CMAKE_COMMAND} -E copy_directory
+    "${CMAKE_CURRENT_SOURCE_DIR}/skalowmini-coef.tm" "${CMAKE_CURRENT_BINARY_DIR}/skalowmini-coef.tm")
+
+
 #------------------------------------------------------------------------------
 # comparison-oskar knits together the cpp code and the python scripts     
-add_custom_target(comparison-oskar
-     COMMAND ${CMAKE_COMMAND} -E env EXTRA_PATH="${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_SOURCE_DIR}/scripts/coeff_scripts"
+add_custom_target(comparison-oskar-basefunctions
+     COMMAND ${CMAKE_COMMAND} -E env
+         EXTRA_PATH="${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_SOURCE_DIR}/scripts/coeff_scripts"
      "${CMAKE_CURRENT_SOURCE_DIR}/generate_basefunction_plots.sh"
-     DEPENDS comparison-oskar-generate-beampattern
+     DEPENDS make_element_response_image
 )
 
-add_test(NAME comparison-oskar-test
-     CONFIGURATIONS integration
-     COMMAND make comparison-oskar)
+add_custom_target(comparison-oskar-station-response
+     COMMAND ${CMAKE_COMMAND} -E env
+         EXTRA_PATH="${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_SOURCE_DIR}/scripts/coeff_scripts:${CMAKE_SOURCE_DIR}/scripts/misc"
+     "${CMAKE_CURRENT_SOURCE_DIR}/compare_stationresponse.sh"
+     DEPENDS make_station_response_image
+)
 

demo/comparison-oskar/compare_stationresponse.py

+import os
+import sys
+import numpy as np
+from read_oskar_beams import read_oskar_beams
+import run_oskar_simulation
+import subprocess
+
+
+tolerance = float(os.environ["TOLERANCE"]) if "TOLERANCE" in os.environ else 0.0
+npixels = int(os.environ["NPIXELS"]) if "NPIXELS" in os.environ else 256
+
+
+run_oskar_simulation.main(npixels)
+subprocess.check_call(["add_beaminfo.py", "skalowmini-coef.MS", "skalowmini-coef.tm"])
+subprocess.check_call(["oskar_csv_to_hdf5.py", "skalowmini-coef.tm", "oskar.h5"])
+subprocess.check_call(["./make_station_response_image", str(npixels)])
+
+A = read_oskar_beams()
+B = np.load('station-response.npy')
+
+if tolerance:
+    difference = np.nanmax(np.abs(A - B))
+    if difference > tolerance:
+        sys.exit(
+            "Difference between OSKAR and EveryBeam spherical wave model is {}, which is larger than the tolerance {}".format(
+                difference, tolerance
+            )
+        )
+

demo/comparison-oskar/compare_stationresponse.sh

+#!/bin/sh
+
+export PATH=$EXTRA_PATH:$PATH
+python3 -B `dirname "${0}"`/compare_stationresponse.py
+

demo/comparison-oskar/generate_basefunction_plots.py

@@ -68,7 +68,7 @@ for em_idx in range(2):
             generate_oskar_csv(l * l - 1 + l - m, em_idx)
 
         subprocess.check_call(["oskar_csv_to_hdf5.py", "telescope.tm", "oskar.h5"])
-        subprocess.check_call(["comparison-oskar-generate-beampattern", str(npixels)])
+        subprocess.check_call(["make_element_response_image", str(npixels)])
 
         B = np.load("response.npy")
 

demo/comparison-oskar/main.cpp → demo/comparison-oskar/make_element_response_image.cpp

@@ -4,8 +4,7 @@
 
 #include <oskarelementresponse.h>
 
-#include "../../external/npy.hpp"
-// #include "npy.hpp"  // to save arrays in numpy format
+#include "../../external/npy.hpp"   // to save arrays in numpy format
 
 int main(int argc, char** argv){
 // int main() {