From 19628f8ec83d0f3e567dc7351270fd6dde28938a Mon Sep 17 00:00:00 2001
From: Sebastiaan van der Tol <tol@astron.nl>
Date: Wed, 5 Aug 2020 08:35:18 +0200
Subject: [PATCH] Add station response to comparision-oskar
---
cpp/CMakeLists.txt | 9 +
cpp/load.cc | 8 +
cpp/load.h | 6 +-
cpp/msv3readutils.cc | 199 +++++++++++++++++
cpp/msv3readutils.h | 74 ++++++
cpp/station.h | 2 +
cpp/telescope/CMakeLists.txt | 2 +
cpp/telescope/oskar.cc | 34 +++
cpp/telescope/oskar.h | 65 ++++++
demo/comparison-oskar/CMakeLists.txt | 8 +-
.../generate_basefunction_plots.sh | 1 +
demo/comparison-oskar/main.cpp | 3 +-
demo/comparison-oskar/stationresponse.cpp | 211 ++++++++++++++++++
scripts/misc/add_beaminfo.py | 136 +++++++++++
14 files changed, 753 insertions(+), 5 deletions(-)
create mode 100644 cpp/msv3readutils.cc
create mode 100644 cpp/msv3readutils.h
create mode 100644 cpp/telescope/oskar.cc
create mode 100644 cpp/telescope/oskar.h
create mode 100644 demo/comparison-oskar/stationresponse.cpp
create mode 100755 scripts/misc/add_beaminfo.py
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index d491e047..182a1ef5 100644
--- a/cpp/CMakeLists.txt
+++ b/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
diff --git a/cpp/load.cc b/cpp/load.cc
index e41686ba..f34aa3fd 100644
--- a/cpp/load.cc
+++ b/cpp/load.cc
@@ -25,6 +25,8 @@ TelescopeType GetTelescopeType(const casacore::MeasurementSet &ms) {
return kATCATelescope;
else if (telescope_name == "MWA")
return kMWATelescope;
+ else if (str.rfind("OSKAR", 0) == 0)
+ return kOskarTelescope;
else
return kUnknownTelescope;
}
@@ -58,6 +60,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, model, options));
+ return telescope;
+ }
default:
casacore::ScalarColumn<casacore::String> telescope_name_col(
ms.observation(), "TELESCOPE_NAME");
diff --git a/cpp/load.h b/cpp/load.h
index 2fdc0ae0..5b470d69 100644
--- a/cpp/load.h
+++ b/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
};
/**
@@ -65,4 +67,4 @@ std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
const Options &options);
} // namespace everybeam
-#endif // EVERYBEAM_LOAD_H_
\ No newline at end of file
+#endif // EVERYBEAM_LOAD_H_
diff --git a/cpp/msv3readutils.cc b/cpp/msv3readutils.cc
new file mode 100644
index 00000000..2c65c112
--- /dev/null
+++ b/cpp/msv3readutils.cc
@@ -0,0 +1,199 @@
+// 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);
+
+inline Antenna::CoordinateSystem ReadCoordinateSystem(
+ const casacore::Table &table, unsigned int id) {
+ casacore::ArrayQuantColumn<casacore::Double> c_position(table, "POSITION",
+ "m");
+ casacore::ArrayQuantColumn<casacore::Double> c_axes(table, "COORDINATE_SYSTEM",
+ "m");
+
+ // Read antenna field center (ITRF).
+ casacore::Vector<casacore::Quantity> aips_position = c_position(id);
+ assert(aips_position.size() == 3);
+
+ vector3r_t position = {{aips_position(0).getValue(),
+ aips_position(1).getValue(),
+ aips_position(2).getValue()}};
+
+ // Read antenna field coordinate axes (ITRF).
+ casacore::Matrix<casacore::Quantity> aips_axes = c_axes(id);
+ assert(aips_axes.shape().isEqual(casacore::IPosition(2, 3, 3)));
+
+ vector3r_t p = {{aips_axes(0, 0).getValue(), aips_axes(1, 0).getValue(),
+ aips_axes(2, 0).getValue()}};
+ vector3r_t q = {{aips_axes(0, 1).getValue(), aips_axes(1, 1).getValue(),
+ aips_axes(2, 1).getValue()}};
+ vector3r_t r = {{aips_axes(0, 2).getValue(), aips_axes(1, 2).getValue(),
+ aips_axes(2, 2).getValue()}};
+
+ Antenna::CoordinateSystem coordinate_system = {position, {p, q, r}};
+ return coordinate_system;
+}
+
+BeamFormer::Ptr ReadMSv3AntennaField(const Table &table, unsigned int id,
+ ElementResponse::Ptr element_response) {
+ Antenna::CoordinateSystem coordinate_system =
+ ReadCoordinateSystem(table, id);
+ BeamFormer::Ptr beam_former(new BeamFormerIdenticalAntennas(coordinate_system));
+// BeamFormer::Ptr beam_former(new BeamFormer(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);
+ std::cout << aips_offset.shape() << std::endl;
+ std::cout << IPosition(2, aips_offset.nrow(),3) << std::endl;
+
+ 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) {
+ std::cout << i << std::endl;
+ 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;
+}
+
+vector3r_t ReadStationPhaseReference(const Table &table, unsigned int id);
+// {
+// vector3r_t phase_reference = {0.0, 0.0, 0.0};
+// const string columnName("LOFAR_PHASE_REFERENCE");
+// if (common::HasColumn(table, columnName)) {
+// ROScalarMeasColumn<MPosition> c_reference(table, columnName);
+// MPosition mReference =
+// MPosition::Convert(c_reference(id), MPosition::ITRF)();
+// MVPosition mvReference = mReference.getValue();
+// phase_reference = {mvReference(0), mvReference(1), mvReference(2)};
+// }
+// return phase_reference;
+// }
+
+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));
+
+ // Read phase reference position (if available).
+// station->SetPhaseReference(ReadStationPhaseReference(ms.antenna(), id));
+
+ 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());
+
+ // TODO
+ // If There is only one field, the top level beamformer is not needed
+ // and the station antenna can be set the the beamformer of the field
+
+ station->SetAntenna(beam_former);
+
+ size_t field_id = 0;
+ size_t element_id = 0;
+ Antenna::CoordinateSystem coordinate_system =
+ ReadCoordinateSystem(tab_phased_array, field_id);
+ auto element_response = station->GetElementResponse();
+ // TODO: rotate coordinate system for antenna
+ auto element =
+ Element::Ptr(new Element(coordinate_system, element_response, element_id));
+ station->SetElement(element);
+
+ return station;
+}
+
+} // namespace everybeam
diff --git a/cpp/msv3readutils.h b/cpp/msv3readutils.h
new file mode 100644
index 00000000..c7738e72
--- /dev/null
+++ b/cpp/msv3readutils.h
@@ -0,0 +1,74 @@
+// lofarreadutils.h: 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$
+
+#ifndef EVERYBEAM_MSV3READUTILS_H_
+#define EVERYBEAM_MSV3READUTILS_H_
+
+// \file
+// Utility functions to read the meta data relevant for simulating the beam from
+// LOFAR observations 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 ReadLofarStation 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_
diff --git a/cpp/station.h b/cpp/station.h
index 9f47426e..51327638 100644
--- a/cpp/station.h
+++ b/cpp/station.h
@@ -333,6 +333,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; }
diff --git a/cpp/telescope/CMakeLists.txt b/cpp/telescope/CMakeLists.txt
index 81dbc809..abbe58da 100644
--- a/cpp/telescope/CMakeLists.txt
+++ b/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")
diff --git a/cpp/telescope/oskar.cc b/cpp/telescope/oskar.cc
new file mode 100644
index 00000000..80b28eda
--- /dev/null
+++ b/cpp/telescope/oskar.cc
@@ -0,0 +1,34 @@
+#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 ElementResponseModel model,
+ const Options &options)
+ : Telescope(ms, model, options) {
+ ReadAllStations(ms, model);
+
+}
+
+std::unique_ptr<griddedresponse::GriddedResponse> OSKAR::GetGriddedResponse(
+ const coords::CoordinateSystem &coordinate_system) {
+ // 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;
+}
diff --git a/cpp/telescope/oskar.h b/cpp/telescope/oskar.h
new file mode 100644
index 00000000..6794cff8
--- /dev/null
+++ b/cpp/telescope/oskar.h
@@ -0,0 +1,65 @@
+// LOFARTelescope.h: Base class for computing the response for the LOFAR
+// telescope.
+//
+// Copyright (C) 2020
+// ASTRON (Netherlands Institute for Radio Astronomy)
+// P.O.Box 2, 7990 AA Dwingeloo, The Netherlands
+//
+// This file is part of the EveryBeam software suite.
+// The EveryBeam software suite is free software: you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as published
+// by the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// The EveryBeam software suite is distributed in the hope that it will be
+// useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with the EveryBeam software suite. If not, see
+// <http://www.gnu.org/licenses/>.
+//
+// $Id$
+
+#ifndef EVERYBEAM_TELESCOPE_OSKAR_H_
+#define EVERYBEAM_TELESCOPE_OSKAR_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 {
+
+//! LOFAR 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 ElementResponseModel model,
+ const Options &options);
+
+ std::unique_ptr<griddedresponse::GriddedResponse> GetGriddedResponse(
+ const coords::CoordinateSystem &coordinate_system) override;
+
+ private:
+ Station::Ptr ReadStation(const casacore::MeasurementSet &ms,
+ const std::size_t id,
+ const ElementResponseModel model) const override;
+};
+} // namespace telescope
+} // namespace everybeam
+
+#endif // EVERYBEAM_TELESCOPE_OSKAR_H_
diff --git a/demo/comparison-oskar/CMakeLists.txt b/demo/comparison-oskar/CMakeLists.txt
index 3ac1b9f1..68e1dd6a 100644
--- a/demo/comparison-oskar/CMakeLists.txt
+++ b/demo/comparison-oskar/CMakeLists.txt
@@ -1,9 +1,15 @@
#------------------------------------------------------------------------------
# CMake file for compiling a comparison between OSKAR and EveryBeam
add_executable(comparison-oskar-generate-beampattern main.cpp)
-
target_link_libraries(comparison-oskar-generate-beampattern oskar)
+add_executable(comparison-oskar-station-response stationresponse.cpp)
+target_link_libraries(comparison-oskar-station-response everybeam)
+
+# Required to get the config.h header
+target_include_directories(comparison-oskar-station-response PRIVATE "${CMAKE_BINARY_DIR}")
+
+
file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/layout.txt"
DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/telescope.tm")
diff --git a/demo/comparison-oskar/generate_basefunction_plots.sh b/demo/comparison-oskar/generate_basefunction_plots.sh
index 9be17d4a..906fa253 100755
--- a/demo/comparison-oskar/generate_basefunction_plots.sh
+++ b/demo/comparison-oskar/generate_basefunction_plots.sh
@@ -1,5 +1,6 @@
#!/bin/sh
export PATH=$EXTRA_PATH:$PATH
+export TOLERANCE="1e-12"
python3 -B `dirname "${0}"`/generate_basefunction_plots.py
diff --git a/demo/comparison-oskar/main.cpp b/demo/comparison-oskar/main.cpp
index a8ea22d5..a10ba951 100644
--- a/demo/comparison-oskar/main.cpp
+++ b/demo/comparison-oskar/main.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() {
diff --git a/demo/comparison-oskar/stationresponse.cpp b/demo/comparison-oskar/stationresponse.cpp
new file mode 100644
index 00000000..a38322fd
--- /dev/null
+++ b/demo/comparison-oskar/stationresponse.cpp
@@ -0,0 +1,211 @@
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+
+#include <oskarelementresponse.h>
+
+#include "../../external/npy.hpp" // to save arrays in numpy format
+
+#include "load.h"
+#include "options.h"
+#include "config.h"
+
+using namespace everybeam;
+
+int main(int argc, char** argv){
+
+ ElementResponseModel response_model = ElementResponseModel::kOSKARSphericalWave;
+ Options options;
+
+ casacore::MeasurementSet ms("/home/vdtol/skalowmini/skalowmini-coef1.MS");
+
+ // Load OSKAR Telescope
+ std::unique_ptr<telescope::Telescope> telescope =
+ Load(ms, response_model, options);
+
+ Station::Ptr station = telescope->GetStation(0);
+
+ Antenna::Ptr antenna = station->GetAntenna();
+
+ auto p = antenna->coordinate_system_.axes.p;
+ auto q = antenna->coordinate_system_.axes.q;
+ auto r = antenna->coordinate_system_.axes.r;
+
+ const vector3r_t station0 = r;
+ const vector3r_t tile0 = r;
+
+
+ double freq = 50e6;
+
+ int N;
+ if (argc == 1){
+ N = 256;
+ }
+ else{
+ N = atoi(argv[1]);
+ }
+
+ std::vector<std::complex<double>> result(N*N*2*2);
+ typedef std::complex<double>result_arr_t[N][N][2][2];
+ result_arr_t &result_arr = * (result_arr_t*) result.data();
+
+
+ for(int i=0; i<N; ++i) {
+ std::cout << i << std::endl;
+ double x = (2.0*i)/(N-1) - 1.0;
+ for(int j=0; j<N; ++j) {
+ double y = (2.0*j)/(N-1) - 1.0;
+
+ double z = sqrt(1.0 - x*x - y*y);
+
+
+// double theta = asin(sqrt(x*x + y*y));
+// double phi = atan2(y,x);
+
+ real_t time = 0;
+ real_t freq = 50e6;
+ const vector3r_t direction = {
+ x*p[0] + y*q[0] + z*r[0],
+ x*p[1] + y*q[1] + z*r[1],
+ x*p[2] + y*q[2] + z*r[2]
+ };
+ real_t freq0 = 50e6;
+
+ auto result = station->Response(time, freq, direction, freq0, station0, tile0);
+ result_arr[i][j][0][0] = result[0][0];
+ result_arr[i][j][0][1] = result[0][1];
+ result_arr[i][j][1][0] = result[1][0];
+ result_arr[i][j][1][1] = result[1][1];
+
+// element_response.Response(0, freq, theta, phi, result_arr[i][j]);
+ }
+ }
+
+ const long unsigned leshape [] = {(long unsigned int) N, (long unsigned int) N, 2, 2};
+ npy::SaveArrayAsNumpy("station-response.npy", false, 4, leshape, result);
+}
+
+
+/*
+
+#include "../options.h"
+#include "../griddedresponse/lofargrid.h"
+#include "../elementresponse.h"
+#include "../../external/npy.hpp"
+
+#include <complex>
+#include <cmath>
+
+using namespace everybeam;
+
+BOOST_AUTO_TEST_CASE(load_lofar) {
+ ElementResponseModel response_model = ElementResponseModel::kHamaker;
+ Options options;
+
+ casacore::MeasurementSet ms(LOFAR_MOCK_MS);
+
+ // Load LOFAR Telescope
+ std::unique_ptr<telescope::Telescope> telescope =
+ Load(ms, response_model, options);
+
+ // Assert if we indeed have a LOFAR pointer
+ BOOST_CHECK(nullptr != dynamic_cast<telescope::LOFAR*>(telescope.get()));
+
+ // Assert if correct number of stations
+ std::size_t nstations = 70;
+ BOOST_CHECK_EQUAL(telescope->GetNrStations(), nstations);
+
+ // Assert if GetStation(stationd_id) behaves properly
+ BOOST_CHECK_EQUAL(telescope->GetStation(0)->GetName(), "CS001HBA0");
+
+ // Properties extracted from MS
+ double time = 4929192878.008341;
+ double frequency = 138476562.5;
+ std::size_t width(4), height(4);
+ double ra(2.15374123), dec(0.8415521), dl(0.5 * M_PI / 180.),
+ dm(0.5 * 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::LOFARGrid*>(grid_response.get()));
+
+ // Define buffer and get gridded responses
+ std::vector<std::complex<float>> antenna_buffer_single(
+ grid_response->GetBufferSize(1));
+ grid_response->CalculateStation(antenna_buffer_single.data(), time, frequency,
+ 23);
+ BOOST_CHECK_EQUAL(antenna_buffer_single.size(),
+ std::size_t(width * height * 2 * 2));
+
+ // LOFARBeam output at pixel (2,2):
+ std::vector<std::complex<float>> lofar_p22 = {{-0.175908, -0.000478397},
+ {-0.845988, -0.00121503},
+ {-0.89047, -0.00125383},
+ {0.108123, -5.36076e-05}};
+
+ // Compare with everybeam
+ std::size_t offset_22 = (2 + 2 * height) * 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);
+ }
+
+ // LOFARBeam output at pixel (1,3):
+ std::vector<std::complex<float>> lofar_p13 = {{-0.158755, -0.000749433},
+ {-0.816165, -0.00272568},
+ {-0.863389, -0.00283979},
+ {0.0936919, 0.000110673}};
+
+ // Compare with everybeam
+ std::size_t offset_13 = (1 + 3 * height) * 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);
+ }
+
+ // 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;
+ options_diff_beam.use_differential_beam = true;
+
+ // Load LOFAR Telescope
+ std::unique_ptr<telescope::Telescope> telescope_diff_beam =
+ Load(ms, response_model, options_diff_beam);
+
+ std::unique_ptr<griddedresponse::GriddedResponse> grid_response_diff_beam =
+ telescope_diff_beam->GetGriddedResponse(coord_system);
+
+ std::vector<std::complex<float>> antenna_buffer_diff_beam(
+ grid_response_diff_beam->GetBufferSize(1));
+ grid_response_diff_beam->CalculateStation(antenna_buffer_diff_beam.data(),
+ time, frequency, 15);
+
+ double norm_jones_mat = 0.;
+ for (std::size_t i = 0; i < 4; ++i) {
+ 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);
+}
+
+*/
diff --git a/scripts/misc/add_beaminfo.py b/scripts/misc/add_beaminfo.py
new file mode 100755
index 00000000..d1ff1ed0
--- /dev/null
+++ b/scripts/misc/add_beaminfo.py
@@ -0,0 +1,136 @@
+#!/usr/bin/env python3
+
+""" Add beam info to OSKAR generated MeasurementSet """
+
+
+import casacore.tables as pt
+import pandas as pd
+from lofarantpos import geo
+from numpy.linalg import norm
+import numpy as np
+import logging
+import sys
+
+
+def read_telescope_center(oskar_telescope_dir: str):
+ """Read telescope center from OSKAR metadata and convert to ITRF"""
+ oskar_telescope_center = pd.read_csv(f"{oskar_telescope_dir}/position.txt",
+ header=None, sep=" ")
+
+ if len(oskar_telescope_center.columns) == 2:
+ oskar_telescope_center["Height"] = [0.]
+
+ oskar_telescope_center.columns = ["Lon", "Lat", "Height"]
+ telescope_center_itrf = geo.xyz_from_geographic(np.deg2rad(oskar_telescope_center["Lon"][0]),
+ np.deg2rad(oskar_telescope_center["Lat"][0]),
+ oskar_telescope_center["Height"][0])
+ return telescope_center_itrf
+
+
+def fix_antenna(oskar_ms_name: str, telescope_center_itrf: np.array):
+ """Make POSITION in ::ANTENNA subtable absolute"""
+ anttable = pt.table(f"{oskar_ms_name}::ANTENNA", readonly=False, ack=False)
+ if np.isclose(norm(anttable[0]["POSITION"]), 0, atol=1e-3):
+ logging.info("Positions are already absolute")
+ return
+ else:
+ anttable.putcol("POSITION", anttable.getcol("POSITION") + telescope_center_itrf)
+ anttable.close()
+
+
+def add_array_center(oskar_ms_name: str, telescope_center_itrf: np.array):
+ """Add ARRAY_CENTER column to ::OBSERVATION subtable"""
+ anttable = pt.table(f"{oskar_ms_name}::ANTENNA", ack=False)
+ coldesc = anttable.getcoldesc("POSITION")
+ coldesc["name"] = "ARRAY_CENTER"
+ coldesc["comment"] = "Reference position for array"
+
+ obstable = pt.table(f"{oskar_ms_name}::OBSERVATION", readonly=False, ack=False)
+ #obstable.addcols(coldesc)
+
+ obstable.putcol("ARRAY_CENTER", np.array([telescope_center_itrf]))
+ obstable.close()
+
+
+def add_phased_array_table(oskar_ms_name: str):
+ """Add PHASED_ARRAY subtable to measurement set"""
+
+ anttable = pt.table(f"{oskar_ms_name}::ANTENNA", ack=False)
+ phasedarraytable = pt.table(f"{oskar_ms_name}/PHASED_ARRAY",
+ pt.maketabdesc([]),
+ nrow=anttable.nrows())
+
+ oskar_ms = pt.table(f"{oskar_ms_name}", readonly=False, ack=False)
+ oskar_ms.putkeyword("PHASED_ARRAY", phasedarraytable, makesubrecord=True)
+ oskar_ms.close()
+
+ position_coldesc = anttable.getcoldesc("POSITION")
+ position_coldesc['comment'] = 'Position of antenna ï¬eld'
+ position_coldesc['name'] = 'POSITION'
+ phasedarraytable.addcols(position_coldesc)
+
+ coordinate_system_coldesc = pt.makearrcoldesc("COORDINATE_SYSTEM", 0.,
+ shape=[3, 3], comment="Local coordinate system",
+ valuetype='double')
+ phasedarraytable.addcols(coordinate_system_coldesc)
+ pt.taql("UPDATE $phasedarraytable SET COORDINATE_SYSTEM=0.");
+
+ element_offset_coldesc = pt.makearrcoldesc("ELEMENT_OFFSET", 0.,
+ ndim=2, comment="Offset per element",
+ valuetype='double',
+ keywords={"MEASINFO": {"type": "position",
+ "Ref": "ITRF"}})
+ phasedarraytable.addcols(element_offset_coldesc)
+
+ element_flag_coldesc = pt.makearrcoldesc("ELEMENT_FLAG", 0.,
+ ndim=2, comment="Offset per element",
+ valuetype='bool')
+ phasedarraytable.addcols(element_flag_coldesc)
+ phasedarraytable.close()
+
+
+def fill_phased_array(oskar_ms_name: str, oskar_telescope_dir: str):
+ """Fill the ::PHASED_ARRAY subtable with info from the OSKAR directory"""
+ element_locations = pd.read_csv(f"{oskar_telescope_dir}/station/layout.txt",
+ header=None, sep=" ")
+ if len(element_locations.columns) == 2:
+ element_locations["Z"] = np.zeros_like(element_locations[0])
+
+ oskar_telescope_center = pd.read_csv(f"{oskar_telescope_dir}/position.txt",
+ header=None, sep=" ")
+ oskar_telescope_center.columns = ["Lon", "Lat"]
+ telescope_center_itrf = read_telescope_center(oskar_telescope_dir)
+
+ normal_vector_ellipsoid = geo.normal_vector_ellipsoid(np.deg2rad(oskar_telescope_center["Lon"][0]),
+ np.deg2rad(oskar_telescope_center["Lat"][0]))
+
+ local_to_itrf_projection_matrix = geo.projection_matrix(telescope_center_itrf,
+ normal_vector_ellipsoid)
+
+ element_locations_itrf = local_to_itrf_projection_matrix @ element_locations.to_numpy().T
+ nr_pol = 2
+ all_unflagged = np.zeros((nr_pol, element_locations_itrf.shape[1]))
+
+ phasedarraytable = pt.table(f"{oskar_ms_name}::PHASED_ARRAY", readonly=False, ack=False)
+
+ for stationnr in range(len(phasedarraytable)):
+ phasedarraytable.putcell("COORDINATE_SYSTEM", stationnr, local_to_itrf_projection_matrix.T)
+ phasedarraytable.putcell("ELEMENT_OFFSET", stationnr, element_locations_itrf)
+ phasedarraytable.putcell("ELEMENT_FLAG", stationnr, all_unflagged)
+
+ anttable = pt.table(f"{oskar_ms_name}::ANTENNA", ack=False)
+ phasedarraytable.putcol("POSITION", anttable.getcol("POSITION"))
+ phasedarraytable.close()
+
+
+def main(oskar_ms_name: str, oskar_telescope_dir: str):
+ """Add beam info to OSKAR generated MeasurementSet"""
+ telescope_center_itrf = read_telescope_center(oskar_telescope_dir)
+ fix_antenna(oskar_ms_name, telescope_center_itrf)
+ add_array_center(oskar_ms_name, telescope_center_itrf)
+ add_phased_array_table(oskar_ms_name)
+ fill_phased_array(oskar_ms_name, oskar_telescope_dir)
+
+
+if __name__ == "__main__":
+ main(sys.argv[1], sys.argv[2])
--
GitLab