From 1b02e3ca1713b95ed6ec03c447ad2d7da8b34eb4 Mon Sep 17 00:00:00 2001
From: Bas van der Tol <tol@astron.nl>
Date: Wed, 16 Sep 2020 08:48:43 +0200
Subject: [PATCH] Resolve "DP3 tApplyBeam test fails due to bug in
Station::ComputeElementResponse"
---
cpp/antenna.cc | 33 ++++++++++++++++++++++-
cpp/antenna.h | 31 ++++++++++++++++++++++
cpp/beamformer.cc | 19 ++++++++++++++
cpp/beamformer.h | 20 ++++++++++++++
cpp/beamformeridenticalantennas.cc | 7 +++++
cpp/beamformeridenticalantennas.h | 2 ++
cpp/element.cc | 9 +++++++
cpp/element.h | 2 ++
cpp/lofarreadutils.cc | 9 -------
cpp/msv3readutils.cc | 9 -------
cpp/station.cc | 42 ++++++++++++++++++++----------
cpp/station.h | 7 ++---
12 files changed, 152 insertions(+), 38 deletions(-)
diff --git a/cpp/antenna.cc b/cpp/antenna.cc
index 6f100299..0d99cd30 100644
--- a/cpp/antenna.cc
+++ b/cpp/antenna.cc
@@ -3,6 +3,7 @@
#include "common/mathutils.h"
namespace everybeam {
+
vector3r_t Antenna::TransformToLocalDirection(const vector3r_t &direction) {
vector3r_t local_direction{
dot(coordinate_system_.axes.p, direction),
@@ -12,4 +13,34 @@ vector3r_t Antenna::TransformToLocalDirection(const vector3r_t &direction) {
return local_direction;
}
-} // namespace everybeam
\ No newline at end of file
+
+void Antenna::Transform(const CoordinateSystem &coordinate_system) {
+ coordinate_system_.axes.p =
+ coordinate_system_.axes.p[0] * coordinate_system.axes.p +
+ coordinate_system_.axes.p[1] * coordinate_system.axes.q +
+ coordinate_system_.axes.p[2] * coordinate_system.axes.r;
+
+ coordinate_system_.axes.q =
+ coordinate_system_.axes.q[0] * coordinate_system.axes.p +
+ coordinate_system_.axes.q[1] * coordinate_system.axes.q +
+ coordinate_system_.axes.q[2] * coordinate_system.axes.r;
+
+ coordinate_system_.axes.r =
+ coordinate_system_.axes.r[0] * coordinate_system.axes.p +
+ coordinate_system_.axes.r[1] * coordinate_system.axes.q +
+ coordinate_system_.axes.r[2] * coordinate_system.axes.r;
+
+ coordinate_system_.origin =
+ coordinate_system.origin +
+ coordinate_system_.origin[0] * coordinate_system.axes.p +
+ coordinate_system_.origin[1] * coordinate_system.axes.q +
+ coordinate_system_.origin[2] * coordinate_system.axes.r;
+
+ phase_reference_position_ =
+ coordinate_system.origin +
+ phase_reference_position_[0] * coordinate_system.axes.p +
+ phase_reference_position_[1] * coordinate_system.axes.q +
+ phase_reference_position_[2] * coordinate_system.axes.r;
+}
+
+} // namespace everybeam
diff --git a/cpp/antenna.h b/cpp/antenna.h
index ebc6658e..e4a76d62 100644
--- a/cpp/antenna.h
+++ b/cpp/antenna.h
@@ -115,6 +115,37 @@ class Antenna {
phase_reference_position_(phase_reference_position),
enabled_{true, true} {}
+ /**
+ * @brief Makes a copy of this Antenna object
+ *
+ * The method is virtual, so that copies can be created from a pointer
+ * to the base (Antenna) class.
+ * The original remains unchanged, therefore the method is const.
+ * The method has no implementation in the Antenna class, because
+ * Antenna is abstract, so no copy can be instantiated.
+ *
+ * This method is used by the ExtractAntenna method of the BeamFormer
+ * class to create a copy of one of the Antennas it contains.
+ */
+ virtual Ptr Clone() const = 0;
+
+ /**
+ * @brief Transform internal coordinate systems and positions
+ *
+ * @param coordinate_system to apply in the transformation
+ *
+ * This method is used by BeamFormer::ExtractAntenna to lift
+ * an antenna out of the beamformer.
+ *
+ * The transformation is needed because the coordinate system of
+ * an antenna in a beamformer is expressed in terms of
+ * the coordinate system of the beamformer.
+ * To turn an embedded antenna into a stand-alone antenna,
+ * the coordinate system of the beamformer needs to be
+ * applied to the coordinate system of the antenna
+ */
+ void Transform(const CoordinateSystem &coordinate_system);
+
/**
* @brief Compute the %Antenna Response
*
diff --git a/cpp/beamformer.cc b/cpp/beamformer.cc
index daf3048c..8d887ffe 100644
--- a/cpp/beamformer.cc
+++ b/cpp/beamformer.cc
@@ -6,6 +6,25 @@
#include <cmath>
namespace everybeam {
+
+Antenna::Ptr BeamFormer::Clone() const {
+ auto beamformer_clone = BeamFormer::Ptr(
+ new BeamFormer(coordinate_system_, phase_reference_position_));
+
+ // antennas_ is a vector of pointers to Antennas, so
+ // this creates a shallow copy, in the sense that
+ // the antennas are not copied, only the pointers.
+ beamformer_clone->antennas_ = antennas_;
+
+ return beamformer_clone;
+}
+
+Antenna::Ptr BeamFormer::ExtractAntenna(size_t antenna_index) const {
+ Antenna::Ptr antenna = antennas_[antenna_index]->Clone();
+ antenna->Transform(coordinate_system_);
+ return antenna;
+}
+
vector3r_t BeamFormer::TransformToLocalPosition(const vector3r_t &position) {
// Get antenna position relative to coordinate system origin
vector3r_t dposition{position[0] - coordinate_system_.origin[0],
diff --git a/cpp/beamformer.h b/cpp/beamformer.h
index 6bd52d91..a9811e25 100644
--- a/cpp/beamformer.h
+++ b/cpp/beamformer.h
@@ -52,6 +52,8 @@ class BeamFormer : public Antenna {
TransformToLocalPosition(phase_reference_position_);
}
+ Antenna::Ptr Clone() const override;
+
/**
* @brief Add an antenna to the m_antenna array.
*
@@ -59,6 +61,24 @@ class BeamFormer : public Antenna {
*/
void AddAntenna(Antenna::Ptr antenna) { antennas_.push_back(antenna); }
+ /**
+ * @brief Extracts an antenna from the beamformer
+ *
+ * @param antenna_index index of antenna to extact
+ * @returns pointer to a copy of antenna with index antenna_index
+ *
+ * The antenna is extracted such that it can be used stand-alone,
+ * independent of the beamformer. The coordinate system of the extracted
+ * antenna is transformed from internal representation to external
+ * representation by application of the beamformer coordinate system to
+ * the antenna coordinate system.
+ *
+ * The returned antenna can be either an Element or a BeamFormer.
+ *
+ * The beamformer itself remains unchanged.
+ */
+ Antenna::Ptr ExtractAntenna(size_t antenna_index) const;
+
protected:
vector3r_t
local_phase_reference_position_; // in coordinate system of Antenna
diff --git a/cpp/beamformeridenticalantennas.cc b/cpp/beamformeridenticalantennas.cc
index 3aa01a7a..75b3ba99 100644
--- a/cpp/beamformeridenticalantennas.cc
+++ b/cpp/beamformeridenticalantennas.cc
@@ -7,6 +7,13 @@
namespace everybeam {
+Antenna::Ptr BeamFormerIdenticalAntennas::Clone() const {
+ auto beamformer_clone = std::make_shared<BeamFormerIdenticalAntennas>(
+ coordinate_system_, phase_reference_position_);
+ beamformer_clone->antennas_ = antennas_;
+ return beamformer_clone;
+}
+
matrix22c_t BeamFormerIdenticalAntennas::LocalResponse(
real_t time, real_t freq, const vector3r_t &direction,
const Options &options) const {
diff --git a/cpp/beamformeridenticalantennas.h b/cpp/beamformeridenticalantennas.h
index 31185fd9..26aedab2 100644
--- a/cpp/beamformeridenticalantennas.h
+++ b/cpp/beamformeridenticalantennas.h
@@ -35,6 +35,8 @@ class BeamFormerIdenticalAntennas : public BeamFormer {
BeamFormerIdenticalAntennas(vector3r_t phase_reference_position)
: BeamFormer(phase_reference_position) {}
+ Antenna::Ptr Clone() const override;
+
private:
// Compute the BeamFormer response in certain direction of arrival (ITRF, m)
// and return (Jones) matrix of response
diff --git a/cpp/element.cc b/cpp/element.cc
index f3c119ae..5803d1ce 100644
--- a/cpp/element.cc
+++ b/cpp/element.cc
@@ -2,6 +2,15 @@
#include "common/mathutils.h"
namespace everybeam {
+
+Antenna::Ptr Element::Clone() const {
+ auto element_clone =
+ Element::Ptr(new Element(coordinate_system_, element_response_, id_));
+ element_clone->enabled_[0] = enabled_[0];
+ element_clone->enabled_[1] = enabled_[1];
+ return element_clone;
+}
+
matrix22c_t Element::LocalResponse(real_t time, real_t freq,
const vector3r_t &direction, size_t id,
const Options &options) const {
diff --git a/cpp/element.h b/cpp/element.h
index 39c0ec32..6e8e0722 100644
--- a/cpp/element.h
+++ b/cpp/element.h
@@ -32,6 +32,8 @@ class Element : public Antenna {
id_(id),
element_response_(element_response) {}
+ Antenna::Ptr Clone() const override;
+
/**
* @brief Compute the local response of the element.
*
diff --git a/cpp/lofarreadutils.cc b/cpp/lofarreadutils.cc
index 2fa2b397..9344c16b 100644
--- a/cpp/lofarreadutils.cc
+++ b/cpp/lofarreadutils.cc
@@ -313,15 +313,6 @@ Station::Ptr ReadLofarStation(const MeasurementSet &ms, unsigned int id,
station->SetAntenna(beam_former);
- size_t field_id = 0;
- size_t element_id = 0;
- Antenna::CoordinateSystem coordinate_system =
- common::ReadCoordinateSystem(tab_field, field_id);
- auto model = station->GetElementResponse();
- // TODO: rotate coordinate system for antenna
- auto element =
- Element::Ptr(new Element(coordinate_system, model, element_id));
- station->SetElement(element);
} else if (telescope_name == "AARTFAAC") {
ROScalarColumn<String> ant_type_col(common::GetSubTable(ms, "OBSERVATION"),
"AARTFAAC_ANTENNA_TYPE");
diff --git a/cpp/msv3readutils.cc b/cpp/msv3readutils.cc
index ddcfe851..4e343949 100644
--- a/cpp/msv3readutils.cc
+++ b/cpp/msv3readutils.cc
@@ -137,15 +137,6 @@ Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
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;
}
diff --git a/cpp/station.cc b/cpp/station.cc
index 38cd35a6..332d9734 100644
--- a/cpp/station.cc
+++ b/cpp/station.cc
@@ -80,6 +80,21 @@ const vector3r_t &Station::GetPhaseReference() const {
return phase_reference_;
}
+void Station::SetAntenna(Antenna::Ptr antenna) {
+ antenna_ = antenna;
+
+ // The antenna can be either an Element or a BeamFormer
+ // If it is a BeamFormer we recursively extract the first antenna
+ // until we have an Element.
+ // The extraction returns copies so antenna_ remains unchanged.
+ // The element that is found is used in ComputeElementResponse to
+ // compute the element response.
+ while (auto beamformer = std::dynamic_pointer_cast<BeamFormer>(antenna)) {
+ antenna = beamformer->ExtractAntenna(0);
+ }
+ element_ = std::dynamic_pointer_cast<Element>(antenna);
+}
+
// ========================================================
matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
const vector3r_t &direction,
@@ -102,13 +117,20 @@ matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
const vector3r_t &direction,
const bool rotate) const {
- // if (rotate)
- // return itsElement->response(time, freq, direction)
- // * Rotation(time, direction);
- // else
- // return itsElement->response(time, freq, direction);
+ Antenna::Options options;
+ options.rotate = rotate;
- return element_->Response(time, freq, direction);
+ if (options.rotate) {
+ vector3r_t ncp_ = NCP(time);
+ vector3r_t east = normalize(cross(ncp_, direction));
+ vector3r_t north = cross(direction, east);
+ options.east = east;
+ options.north = north;
+ }
+
+ matrix22c_t response;
+ response = element_->Response(time, freq, direction, options);
+ return response;
}
matrix22c_t Station::Response(real_t time, real_t freq,
@@ -129,14 +151,6 @@ matrix22c_t Station::Response(real_t time, real_t freq,
matrix22c_t response = antenna_->Response(time, freq, direction, options);
- // if (rotate) {
- // std::cout << "rotate" << std::endl;
- // auto r = Rotation(time, direction);
- // std::cout << r[0][0] << ", " << r[0][1] << std::endl;
- // std::cout << r[1][0] << ", " << r[1][1] << std::endl;
- // response = response * r;
- // std::cout << response[0][0] << std::endl;
- // }
return response;
}
diff --git a/cpp/station.h b/cpp/station.h
index 85197a71..105b5842 100644
--- a/cpp/station.h
+++ b/cpp/station.h
@@ -333,12 +333,9 @@ class Station {
}
//! Set antenna attribute, usually a BeamFormer, but can also be an Element
- void SetAntenna(Antenna::Ptr antenna) { antenna_ = antenna; }
+ void SetAntenna(Antenna::Ptr antenna);
- Antenna::Ptr GetAntenna() { return antenna_; }
-
- //! Set Element attribute
- void SetElement(Element::Ptr element) { element_ = element; }
+ Antenna::Ptr GetAntenna() const { return antenna_; }
private:
vector3r_t NCP(real_t time) const;
--
GitLab