From 458afa312728ffe051f4ce106d2e24aa28c08b22 Mon Sep 17 00:00:00 2001
From: Sebastiaan van der Tol <tol@astron.nl>
Date: Fri, 7 Aug 2020 06:16:02 +0200
Subject: [PATCH] Fixes after rebase
---
cpp/load.cc | 9 ++--
cpp/msv3readutils.cc | 65 +++++++++++------------
cpp/telescope/oskar.cc | 7 +--
cpp/telescope/oskar.h | 31 ++++++++++-
demo/comparison-oskar/stationresponse.cpp | 16 ++++--
scripts/misc/add_beaminfo.py | 4 +-
6 files changed, 83 insertions(+), 49 deletions(-)
diff --git a/cpp/load.cc b/cpp/load.cc
index f34aa3fd..ad8a9977 100644
--- a/cpp/load.cc
+++ b/cpp/load.cc
@@ -25,8 +25,9 @@ TelescopeType GetTelescopeType(const casacore::MeasurementSet &ms) {
return kATCATelescope;
else if (telescope_name == "MWA")
return kMWATelescope;
- else if (str.rfind("OSKAR", 0) == 0)
- return kOskarTelescope;
+ // check if telescope_name starts with "OSKAR"
+ else if (telescope_name.rfind("OSKAR", 0) == 0)
+ return kOSKARTelescope;
else
return kUnknownTelescope;
}
@@ -60,10 +61,10 @@ std::unique_ptr<telescope::Telescope> Load(casacore::MeasurementSet &ms,
new telescope::MWA(ms, options));
return telescope;
}
- case kOskarTelescope: {
+ case kOSKARTelescope: {
std::unique_ptr<telescope::Telescope> telescope =
std::unique_ptr<telescope::Telescope>(
- new telescope::OSKAR(ms, model, options));
+ new telescope::OSKAR(ms, options));
return telescope;
}
default:
diff --git a/cpp/msv3readutils.cc b/cpp/msv3readutils.cc
index 5ee00ad5..3ad0da8b 100644
--- a/cpp/msv3readutils.cc
+++ b/cpp/msv3readutils.cc
@@ -73,39 +73,39 @@ 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;
-}
+// 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);
+ Antenna::CoordinateSystem coordinate_system = common::ReadCoordinateSystem(table, id);
BeamFormer::Ptr beam_former(
new BeamFormerIdenticalAntennas(coordinate_system));
// BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system));
@@ -115,8 +115,6 @@ BeamFormer::Ptr ReadMSv3AntennaField(const Table &table, unsigned int id,
// 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)));
@@ -124,7 +122,6 @@ BeamFormer::Ptr ReadMSv3AntennaField(const Table &table, unsigned int 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()};
@@ -194,7 +191,7 @@ Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
size_t field_id = 0;
size_t element_id = 0;
Antenna::CoordinateSystem coordinate_system =
- ReadCoordinateSystem(tab_phased_array, field_id);
+ common::ReadCoordinateSystem(tab_phased_array, field_id);
auto element_response = station->GetElementResponse();
// TODO: rotate coordinate system for antenna
auto element = Element::Ptr(
diff --git a/cpp/telescope/oskar.cc b/cpp/telescope/oskar.cc
index f280ec03..d498aa2a 100644
--- a/cpp/telescope/oskar.cc
+++ b/cpp/telescope/oskar.cc
@@ -11,10 +11,11 @@ using namespace everybeam;
using namespace everybeam::telescope;
using namespace casacore;
-OSKAR::OSKAR(MeasurementSet &ms, const ElementResponseModel model,
+OSKAR::OSKAR(MeasurementSet &ms,
const Options &options)
- : Telescope(ms, model, options) {
- ReadAllStations(ms, model);
+ : Telescope(ms, options) {
+ stations_.resize(nstations_);
+ ReadAllStations(ms, options_.element_response_model);
}
std::unique_ptr<griddedresponse::GriddedResponse> OSKAR::GetGriddedResponse(
diff --git a/cpp/telescope/oskar.h b/cpp/telescope/oskar.h
index 62b49fb2..9b85e71e 100644
--- a/cpp/telescope/oskar.h
+++ b/cpp/telescope/oskar.h
@@ -25,6 +25,8 @@
#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>
@@ -46,16 +48,41 @@ class OSKAR final : public Telescope {
* @param model Element Response model
* @param options telescope options
*/
- OSKAR(casacore::MeasurementSet &ms, const ElementResponseModel model,
+ 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 override;
+ const ElementResponseModel model) const;
+
+ std::vector<Station::Ptr> stations_;
+
};
} // namespace telescope
} // namespace everybeam
diff --git a/demo/comparison-oskar/stationresponse.cpp b/demo/comparison-oskar/stationresponse.cpp
index a38322fd..5542e779 100644
--- a/demo/comparison-oskar/stationresponse.cpp
+++ b/demo/comparison-oskar/stationresponse.cpp
@@ -16,14 +16,23 @@ int main(int argc, char** argv){
ElementResponseModel response_model = ElementResponseModel::kOSKARSphericalWave;
Options options;
+ options.element_response_model = response_model;
casacore::MeasurementSet ms("/home/vdtol/skalowmini/skalowmini-coef1.MS");
+ casacore::ScalarMeasColumn<casacore::MDirection> referenceDirColumn(
+ ms.field(),
+ casacore::MSField::columnName(casacore::MSFieldEnums::REFERENCE_DIR));
+
+ auto preapplied_beam_dir = referenceDirColumn(0);
+
+ std::cout << preapplied_beam_dir << std::endl;
+
// Load OSKAR Telescope
- std::unique_ptr<telescope::Telescope> telescope =
- Load(ms, response_model, options);
+ auto telescope = Load(ms, options);
+ auto &oskar_telescope = *dynamic_cast<telescope::OSKAR*>(telescope.get());
- Station::Ptr station = telescope->GetStation(0);
+ Station::Ptr station = oskar_telescope.GetStation(0);
Antenna::Ptr antenna = station->GetAntenna();
@@ -77,7 +86,6 @@ int main(int argc, char** argv){
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]);
}
}
diff --git a/scripts/misc/add_beaminfo.py b/scripts/misc/add_beaminfo.py
index 633144d3..7ecd7cd1 100755
--- a/scripts/misc/add_beaminfo.py
+++ b/scripts/misc/add_beaminfo.py
@@ -72,7 +72,7 @@ def add_phased_array_table(oskar_ms_name: str):
shape=[3, 3], comment="Local coordinate system",
valuetype='double')
phasedarraytable.addcols(coordinate_system_coldesc)
- pt.taql("UPDATE $phasedarraytable SET COORDINATE_SYSTEM=0.");
+ pt.taql("UPDATE $phasedarraytable SET COORDINATE_AXES=0.");
element_offset_coldesc = pt.makearrcoldesc("ELEMENT_OFFSET", 0.,
ndim=2, comment="Offset per element",
@@ -113,7 +113,7 @@ def fill_phased_array(oskar_ms_name: str, oskar_telescope_dir: str):
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("COORDINATE_AXES", stationnr, local_to_itrf_projection_matrix.T)
phasedarraytable.putcell("ELEMENT_OFFSET", stationnr, element_locations_itrf)
phasedarraytable.putcell("ELEMENT_FLAG", stationnr, all_unflagged)
--
GitLab