LOFAR HBA station commits

cpp/CMakeLists.txt

@@ -14,7 +14,9 @@ add_library(everybeam SHARED
   elementresponse.cc
   beamformer.cc
   beamformeridenticalantennas.cc
+  beamformerlofarhba.cc
   element.cc
+  elementhamaker.cc
   load.cc
   common/fftresampler.cc
   coords/itrfconverter.cc
@@ -56,8 +58,10 @@ install (
 install (FILES
   antenna.h
   beamformer.h
+  beamformerlofarhba.h
   beamformeridenticalantennas.h
   element.h
+  elementhamaker.h
   elementresponse.h
   lofarreadutils.h
   msv3readutils.h

cpp/antenna.cc

@@ -4,7 +4,8 @@
 
 namespace everybeam {
 
-vector3r_t Antenna::TransformToLocalDirection(const vector3r_t &direction) {
+vector3r_t Antenna::TransformToLocalDirection(
+    const vector3r_t &direction) const {
   vector3r_t local_direction{
       dot(coordinate_system_.axes.p, direction),
       dot(coordinate_system_.axes.q, direction),

cpp/antenna.h

@@ -155,8 +155,9 @@ class Antenna {
    * @param options
    * @return matrix22c_t Jones matrix
    */
-  matrix22c_t Response(real_t time, real_t freq, const vector3r_t &direction,
-                       const Options &options = {}) {
+  virtual matrix22c_t Response(real_t time, real_t freq,
+                               const vector3r_t &direction,
+                               const Options &options = {}) {
     // Transform direction and directions in options to local coordinatesystem
     vector3r_t local_direction = TransformToLocalDirection(direction);
     Options local_options = {
@@ -206,7 +207,7 @@ class Antenna {
     return {1.0, 1.0};
   }
 
-  vector3r_t TransformToLocalDirection(const vector3r_t &direction);
+  vector3r_t TransformToLocalDirection(const vector3r_t &direction) const;
 };
 
 }  // namespace everybeam

cpp/beamformer.h

@@ -55,7 +55,7 @@ class BeamFormer : public Antenna {
   Antenna::Ptr Clone() const override;
 
   /**
-   * @brief Add an antenna to the m_antenna array.
+   * @brief Add an antenna to the antennas_ array.
    *
    * @param antenna
    */

cpp/beamformerlofarhba.cc

+#include "beamformerlofarhba.h"
+
+#include "common/constants.h"
+#include "common/mathutils.h"
+
+#include <cmath>
+
+namespace everybeam {
+
+Antenna::Ptr BeamFormerLofarHBA::Clone() const {
+  auto beamformer_clone = std::make_shared<BeamFormerLofarHBA>(
+      BeamFormerLofarHBA(coordinate_system_, phase_reference_position_));
+
+  // NOTE: this is an incomplete clone, only creating a deep-copy of the
+  // element. In fact, it also hides an upcast from an ElementHamaker into
+  // an Element object.
+  // The sole and single purpose of Clone() is to be used in
+  // Station::SetAntenna!
+  Element element_copy = *element_;
+  beamformer_clone->SetElement(std::make_shared<Element>(element_copy));
+  return beamformer_clone;
+}
+
+std::vector<std::complex<double>> BeamFormerLofarHBA::ComputeGeometricResponse(
+    std::vector<vector3r_t> phase_reference_positions,
+    const vector3r_t &direction) const {
+  // Initialize and fill result vector by looping over phase_reference_positions
+  std::vector<std::complex<double>> result(phase_reference_positions.size());
+
+  for (std::size_t idx = 0; idx < phase_reference_positions.size(); ++idx) {
+    // Simplified dot product, since local_phase_reference_position_ = [0, 0, 0]
+    const double dl = direction[0] * phase_reference_positions[idx][0] +
+                      direction[1] * phase_reference_positions[idx][1] +
+                      direction[2] * phase_reference_positions[idx][2];
+
+    // Note that the frequency weighting is already implicit in "direction"
+    double phase = -2 * M_PI * dl / common::c;
+    result[idx] = {std::cos(phase), std::sin(phase)};
+  }
+  return result;
+}
+
+diag22c_t BeamFormerLofarHBA::LocalArrayFactor(real_t time, real_t freq,
+                                               const vector3r_t &direction,
+                                               const Options &options) const {
+  diag22c_t result = {0};
+
+  // Compute the array factor of the field
+  diag22c_t array_factor_field =
+      FieldArrayFactor(time, freq, direction, options);
+
+  // Compute the array factor of a tile
+  std::complex<double> array_factor_tile =
+      TileArrayFactor(time, freq, direction, options);
+
+  result[0] = array_factor_tile * array_factor_field[0];
+  result[1] = array_factor_tile * array_factor_field[1];
+
+  return result;
+}
+
+diag22c_t BeamFormerLofarHBA::FieldArrayFactor(real_t time, real_t freq,
+                                               const vector3r_t &direction,
+                                               const Options &options) const {
+  // Weighted subtraction of the directions, with weights given by corresponding
+  // freqs. Purpose is to correctly handle the case in which options.freq0 !=
+  // freq
+  vector3r_t delta_direction =
+      options.freq0 * options.station0 - freq * direction;
+
+  // Get geometric response for pointing direction
+  std::vector<std::complex<double>> geometric_response =
+      ComputeGeometricResponse(tile_positions_, delta_direction);
+
+  double weight_sum[2] = {0.0, 0.0};
+  diag22c_t result;
+
+  for (std::size_t idx = 0; idx < tile_positions_.size(); ++idx) {
+    result[0] += geometric_response[idx] * (1.0 * tile_enabled_[idx][0]);
+    result[1] += geometric_response[idx] * (1.0 * tile_enabled_[idx][1]);
+    weight_sum[0] += (1.0 * tile_enabled_[idx][0]);
+    weight_sum[1] += (1.0 * tile_enabled_[idx][1]);
+  }
+
+  // Normalize the weight by the number of enabled tiles
+  result[0] /= weight_sum[0];
+  result[1] /= weight_sum[1];
+
+  return result;
+}
+
+std::complex<double> BeamFormerLofarHBA::TileArrayFactor(
+    real_t time, real_t freq, const vector3r_t &direction,
+    const Options &options) const {
+  // Weighted subtraction of the directions, with weights given by corresponding
+  // freqs. Purpose is to correctly handle the case in which options.freq0 !=
+  // freq
+  vector3r_t delta_direction = options.freq0 * options.tile0 - freq * direction;
+
+  // Get geometric response for the difference vector stored in "pointing"
+  std::vector<std::complex<double>> geometric_response =
+      ComputeGeometricResponse(element_positions_, delta_direction);
+
+  // Initialize and fill result
+  std::complex<double> result = 0;
+  for (std::size_t idx = 0; idx < element_positions_.size(); ++idx) {
+    result += geometric_response[idx];
+  }
+
+  // Normalize the result by the number of tiles
+  double weight = element_positions_.size();
+  result /= weight;
+
+  return result;
+}
+
+matrix22c_t BeamFormerLofarHBA::LocalResponse(real_t time, real_t freq,
+                                              const vector3r_t &direction,
+                                              const Options &options) const {
+  matrix22c_t result = {0};
+
+  // Compute the combined array factor
+  diag22c_t array_factor = LocalArrayFactor(time, freq, direction, options);
+
+  // NOTE: there are maybe some redundant transformations in element-> response
+  matrix22c_t element_response =
+      element_->Response(time, freq, direction, options);
+
+  result[0][0] = (array_factor[0] * element_response[0][0]);
+  result[0][1] = (array_factor[0] * element_response[0][1]);
+  result[1][0] = (array_factor[1] * element_response[1][0]);
+  result[1][1] = (array_factor[1] * element_response[1][1]);
+
+  return result;
+}
+}  // namespace everybeam

cpp/beamformerlofarhba.h

+#ifndef EVERYBEAM_BEAMFORMERLOFARHBA_H
+#define EVERYBEAM_BEAMFORMERLOFARHBA_H
+
+#include <complex>
+#include <vector>
+#include <memory>
+
+#include "beamformer.h"
+#include "elementhamaker.h"
+#include "common/types.h"
+#include "antenna.h"
+
+namespace everybeam {
+/**
+ * @brief Optimized implementation of the BeamFormer class for the LOFAR HBA
+ * telescope in combination with Hamaker element response model.
+ *
+ */
+class BeamFormerLofarHBA : public Antenna {
+ public:
+  /**
+   * @brief Construct a new BeamFormer object
+   *
+   */
+  BeamFormerLofarHBA() : Antenna() {}
+
+  /**
+   * @brief Construct a new BeamFormerLofarHBA object given a coordinate system.
+   *
+   * @param coordinate_system
+   */
+  BeamFormerLofarHBA(const CoordinateSystem &coordinate_system)
+      : Antenna(coordinate_system) {}
+
+  /**
+   * @brief Construct a new BeamFormerLofarHBA object given a coordinate system
+   * and a phase reference position
+   *
+   * @param coordinate_system
+   * @param phase_reference_position
+   */
+  BeamFormerLofarHBA(CoordinateSystem coordinate_system,
+                     vector3r_t phase_reference_position)
+      : Antenna(coordinate_system, phase_reference_position) {}
+
+  BeamFormerLofarHBA(vector3r_t phase_reference_position)
+      : Antenna(phase_reference_position) {}
+
+  /**
+   * @brief Returns an (incomplete!) clone of the BeamFormerLofarHBA class
+   * only the element_ is copied. This method is intended to be exclusively
+   * used in Station::SetAntenna!
+   *
+   * @return Antenna::Ptr
+   */
+  Antenna::Ptr Clone() const override;
+
+  /**
+   * @brief Set the (unique) Tile for the BeamFormerLofarHBA object.
+   *
+   * @param beamformer
+   */
+  void SetTile(BeamFormer::Ptr beamformer) { tile_ = beamformer; }
+
+  /**
+   * @brief Set the (unique) Element object for the BeamFormerLofarHBA object.
+   *
+   * @param element
+   */
+  void SetElement(Element::Ptr element) { element_ = element; }
+  // void SetElement(Element::Ptr element) { element_ = element; }
+
+  /**
+   * @brief Add tile position to the tile_positions_ array
+   *
+   * @param position
+   */
+  void AddTilePosition(const vector3r_t &position) {
+    tile_positions_.push_back(position);
+  }
+
+  /**
+   * @brief Mark whether tile is enabled by pushing back boolean array to
+   * tile_enabled_ array
+   *
+   * @param enabled
+   */
+  void AddTileEnabled(const std::array<bool, 2> enabled) {
+    tile_enabled_.push_back(enabled);
+  }
+
+  /**
+   * @brief Add element position to the element_positions_ array
+   *
+   * @param position
+   */
+  void AddElementPosition(const vector3r_t &position) {
+    element_positions_.push_back(position);
+  }
+
+  /**
+   * @brief Returns the element_ object of the BeamFormerLofarHBA
+   *
+   * @return std::shared_ptr<Element>
+   */
+  std::shared_ptr<Element> GetElement() const { return element_; };
+
+ private:
+  // Compute the geometric response either for the tiles, or for the elements
+  // within the BeamFormerLofarHBA. Method assumes that the direction is
+  // specified as the (frequency weighted) difference between the pointing_dir
+  // and the probing direction
+  std::vector<std::complex<double>> ComputeGeometricResponse(
+      std::vector<vector3r_t> phase_reference_positions,
+      const vector3r_t &direction) const;
+
+  // Array factor
+  virtual diag22c_t LocalArrayFactor(real_t time, real_t freq,
+                                     const vector3r_t &direction,
+                                     const Options &options) const override;
+
+  // Array factor at Field level, weighting the tiles
+  virtual diag22c_t FieldArrayFactor(real_t time, real_t freq,
+                                     const vector3r_t &direction,
+                                     const Options &options) const;
+
+  // Array Factor at tile level, weighting the elements
+  virtual std::complex<double> TileArrayFactor(real_t time, real_t freq,
+                                               const vector3r_t &direction,
+                                               const Options &options) const;
+
+  // Override of the LocalResponse method
+  virtual matrix22c_t LocalResponse(real_t time, real_t freq,
+                                    const vector3r_t &direction,
+                                    const Options &options) const override;
+
+  // BeamFormerLofarHBA has only one unique tile, and ntiles
+  // unique tile positions
+  std::shared_ptr<BeamFormer> tile_;
+  std::vector<vector3r_t> tile_positions_;
+
+  // Is tile enabled?
+  std::vector<std::array<bool, 2>> tile_enabled_;
+
+  // Each BeamFormerLofarHBA stores one unique element, and 16 unique positions
+  // Usually, the Element will be of type ElementHamaker.
+  std::shared_ptr<Element> element_;
+  std::vector<vector3r_t> element_positions_;
+};
+}  // namespace everybeam
+#endif

cpp/common/constants.h

@@ -31,7 +31,10 @@
 namespace everybeam {
 namespace common {
 /** Speed of light (m/s) */
-const real_t c = 2.99792458e+08;
+constexpr real_t c = 2.99792458e+08;
+
+/** 0.25*pi */
+constexpr real_t pi_4 = 0.7853981633974483096156608;
 }  // namespace common
 }  // namespace everybeam
 

cpp/element.cc

@@ -36,10 +36,4 @@ matrix22c_t Element::LocalResponse(real_t time, real_t freq,
   }
   return result;
 }
-
-matrix22c_t Element::LocalResponse(real_t time, real_t freq,
-                                   const vector3r_t &direction,
-                                   const Options &options) const {
-  return LocalResponse(time, freq, direction, id_, options);
-}
 }  // namespace everybeam

cpp/element.h

@@ -44,14 +44,16 @@ class Element : public Antenna {
    * @param options
    * @return matrix22c_t
    */
-  matrix22c_t LocalResponse(real_t time, real_t freq,
-                            const vector3r_t &direction, size_t id,
-                            const Options &options) const;
-
- private:
-  virtual matrix22c_t LocalResponse(
-      real_t time, real_t freq, const vector3r_t &direction,
-      const Options &options) const final override;
+  virtual matrix22c_t LocalResponse(real_t time, real_t freq,
+                                    const vector3r_t &direction, size_t id,
+                                    const Options &options) const;
+
+ protected:
+  virtual matrix22c_t LocalResponse(real_t time, real_t freq,
+                                    const vector3r_t &direction,
+                                    const Options &options) const override {
+    return LocalResponse(time, freq, direction, id_, options);
+  };
 
   int id_;
   ElementResponse::Ptr element_response_;

cpp/elementhamaker.cc

+#include "elementhamaker.h"
+#include "common/mathutils.h"
+#include "common/constants.h"
+
+namespace everybeam {
+Antenna::Ptr ElementHamaker::Clone() const {
+  auto element_clone = std::make_shared<ElementHamaker>(
+      ElementHamaker(coordinate_system_, element_response_, id_));
+  element_clone->enabled_[0] = enabled_[0];
+  element_clone->enabled_[1] = enabled_[1];
+  return element_clone;
+}
+
+matrix22c_t ElementHamaker::LocalResponse(real_t time, real_t freq,
+                                          const vector3r_t &direction,
+                                          size_t id,
+                                          const Options &options) const {
+  vector2r_t thetaphi = cart2thetaphi(direction);
+  thetaphi[1] -= 5.0 * common::pi_4;
+
+  matrix22c_t result;
+  static_assert(sizeof(std::complex<double>[2][2]) == sizeof(matrix22c_t));
+  element_response_->Response(
+      id, freq, thetaphi[0], thetaphi[1],
+      reinterpret_cast<std::complex<double>(&)[2][2]>(result));
+
+  if (options.rotate) {
+    vector3r_t up = {0.0, 0.0, 1.0};
+    vector3r_t e_phi = normalize(cross(up, direction));
+    vector3r_t e_theta = cross(e_phi, direction);
+    matrix22r_t rotation;
+    rotation[0] = {dot(e_theta, options.north), dot(e_theta, options.east)};
+    rotation[1] = {dot(e_phi, options.north), dot(e_phi, options.east)};
+    result = result * rotation;
+  }
+  return result;
+}
+}  // namespace everybeam

cpp/elementhamaker.h

+#ifndef EVERYBEAM_ELEMENT_HAMAKER_H
+#define EVERYBEAM_ELEMENT_HAMAKER_H
+
+#include <complex>
+#include <memory>
+
+#include "antenna.h"
+#include "element.h"
+#include "elementresponse.h"
+#include "common/types.h"
+
+namespace everybeam {
+
+/**
+ * @brief Elementary antenna, optimized for LOFAR Hamaker model. Derived from
+ * the Element class.
+ *
+ */
+class ElementHamaker : public Element {
+ public:
+  /**
+   * @brief Construct a new Element object
+   *
+   * @param coordinate_system (antenna) CoordinateSystem
+   * @param element_response ElementResponseModel
+   * @param id
+   */
+  ElementHamaker(const CoordinateSystem &coordinate_system,
+                 ElementResponse::Ptr element_response, int id)
+      : Element(coordinate_system, element_response, id) {}
+
+  Antenna::Ptr Clone() const override;
+
+  virtual matrix22c_t Response(real_t time, real_t freq,
+                               const vector3r_t &direction,
+                               const Options &options) override {
+    // The only transform that is needed is hard-coded in LocalResponse
+    matrix22c_t response = LocalResponse(time, freq, direction, options);
+    return response;
+  }
+
+  virtual matrix22c_t LocalResponse(real_t time, real_t freq,
+                                    const vector3r_t &direction, size_t id,
+                                    const Options &options) const override;
+
+ private:
+  virtual matrix22c_t LocalResponse(
+      real_t time, real_t freq, const vector3r_t &direction,
+      const Options &options) const final override {
+    return LocalResponse(time, freq, direction, id_, options);
+  };
+};
+}  // namespace everybeam
+
+#endif

cpp/lofarreadutils.cc

@@ -23,9 +23,12 @@
 
 #include "lofarreadutils.h"
 #include "beamformeridenticalantennas.h"
+#include "beamformerlofarhba.h"
 #include "common/mathutils.h"
 #include "common/casautils.h"
 
+#include <memory>
+
 #include <casacore/measures/Measures/MDirection.h>
 #include <casacore/measures/Measures/MPosition.h>
 #include <casacore/measures/Measures/MCDirection.h>
@@ -144,11 +147,12 @@ vector3r_t TransformToFieldCoordinates(
 //     return system;
 // }
 
-BeamFormer::Ptr MakeTile(unsigned int id, const vector3r_t &position,
-                         const TileConfig &tile_config,
-                         ElementResponse::Ptr element_response) {
-  BeamFormer::Ptr tile =
-      BeamFormer::Ptr(new BeamFormerIdenticalAntennas(position));
+std::shared_ptr<BeamFormer> MakeTile(unsigned int id,
+                                     const vector3r_t &position,
+                                     const TileConfig &tile_config,
+                                     ElementResponse::Ptr element_response) {
+  std::shared_ptr<BeamFormer> tile =
+      std::make_shared<BeamFormer>(BeamFormerIdenticalAntennas(position));
 
   for (unsigned int id = 0; id < tile_config.size(); id++) {
     vector3r_t antenna_position = tile_config[id];
@@ -165,20 +169,35 @@ BeamFormer::Ptr MakeTile(unsigned int id, const vector3r_t &position,
   return tile;
 }
 
-BeamFormer::Ptr ReadAntennaField(const Table &table, unsigned int id,
-                                 ElementResponse::Ptr element_response) {
+// Make a dedicated HBA "Hamaker" tile, saving only one element, and 16
+// element positions
+void MakeTile(std::shared_ptr<BeamFormerLofarHBA> beamformer,
+              const vector3r_t &position, const TileConfig &tile_config,
+              ElementResponse::Ptr element_response) {
+  for (unsigned int id = 0; id < tile_config.size(); id++) {
+    vector3r_t antenna_position = tile_config[id];
+
+    Antenna::CoordinateSystem antenna_coordinate_system;
+    antenna_coordinate_system.origin = antenna_position;
+    antenna_coordinate_system.axes = lofar_antenna_orientation;
+
+    std::shared_ptr<ElementHamaker> antenna = std::make_shared<ElementHamaker>(
+        ElementHamaker(antenna_coordinate_system, element_response, id));
+
+    // Only element 1 needs to be stored as an element
+    if (id == 0) {
+      beamformer->SetElement(antenna);
+    }
+    // All positions need to be stored, however
+    beamformer->AddElementPosition(antenna_position);
+  }
+}
+
+Antenna::Ptr ReadAntennaField(const Table &table, unsigned int id,
+                              ElementResponse::Ptr element_response,
+                              ElementResponseModel element_response_model) {
   Antenna::CoordinateSystem coordinate_system =
       common::ReadCoordinateSystem(table, id);
-  //     std::cout << "coordinate_system: " << std::endl;
-  //     std::cout << "  axes.p: " << coordinate_system.axes.p[0] << ", " <<
-  //     coordinate_system.axes.p[1] << ", " << coordinate_system.axes.p[2] <<
-  //     std::endl; std::cout << "  axes.q: " << coordinate_system.axes.q[0] <<
-  //     ", " << coordinate_system.axes.q[1] << ", " <<
-  //     coordinate_system.axes.q[2] << std::endl; std::cout << "  axes.r: " <<
-  //     coordinate_system.axes.r[0] << ", " << coordinate_system.axes.r[1] <<
-  //     ", " << coordinate_system.axes.r[2] << std::endl;
-  BeamFormer::Ptr beam_former(
-      new BeamFormerIdenticalAntennas(coordinate_system));
 
   ROScalarColumn<String> c_name(table, "NAME");
   ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
@@ -195,6 +214,25 @@ BeamFormer::Ptr ReadAntennaField(const Table &table, unsigned int id,
 
   TileConfig tile_config;
   if (name != "LBA") tile_config = ReadTileConfig(table, id);
+
+  std::shared_ptr<Antenna> beam_former;
+  // Cast to the beam_former corresponding to the element response
+  // model and LBA/HBA configuration
+  if (element_response_model == ElementResponseModel::kHamaker) {
+    if (name != "LBA") {
+      // Then HBA, HBA0 or HBA1
+      beam_former = std::make_shared<BeamFormerLofarHBA>(
+          BeamFormerLofarHBA(coordinate_system));
+    } else {
+      // TODO: will become a dedicated BeamFormerLofarLBA in the near future
+      beam_former = std::make_shared<BeamFormerIdenticalAntennas>(
+          BeamFormerIdenticalAntennas(coordinate_system));
+    }
+  } else {
+    // Tiles / element should be kept unique, so work with generic BeamFormer
+    beam_former = std::make_shared<BeamFormer>(BeamFormer(coordinate_system));
+  }
+
   TransformToFieldCoordinates(tile_config, coordinate_system.axes);
 
   for (size_t i = 0; i < aips_offset.ncolumn(); ++i) {
@@ -206,17 +244,50 @@ BeamFormer::Ptr ReadAntennaField(const Table &table, unsigned int id,
     Antenna::Ptr antenna;
     Antenna::CoordinateSystem antenna_coordinate_system{
         antenna_position, lofar_antenna_orientation};
+
     if (name == "LBA") {
-      antenna = Element::Ptr(
-          new Element(antenna_coordinate_system, element_response, id));
+      antenna = std::make_shared<Element>(
+          Element(antenna_coordinate_system, element_response, id));
+      antenna->enabled_[0] = !aips_flag(0, i);
+      antenna->enabled_[1] = !aips_flag(1, i);
+
+      if (element_response_model == kHamaker) {
+        // NOTE: no cast needed as yet, but it already hints towards
+        // future implementation
+        std::shared_ptr<BeamFormerIdenticalAntennas> beam_former_lba =
+            std::static_pointer_cast<BeamFormerIdenticalAntennas>(beam_former);
+        beam_former_lba->AddAntenna(antenna);
+      } else {
+        std::shared_ptr<BeamFormer> beam_former_lba =
+            std::static_pointer_cast<BeamFormer>(beam_former);
+        beam_former_lba->AddAntenna(antenna);
+      }
     } else {
-      // name is HBA, HBA0, HBA1
-      antenna = MakeTile(id, antenna_position, tile_config, element_response);
+      // name is HBA, HBA0 or HBA1
+      if (element_response_model == kHamaker) {
+        std::shared_ptr<BeamFormerLofarHBA> beam_former_hba =
+            std::static_pointer_cast<BeamFormerLofarHBA>(beam_former);
+
+        // Tile positions are uniques
+        beam_former_hba->AddTilePosition(antenna_position);
+
+        // Store only one tile
+        if (i == 0) {
+          MakeTile(beam_former_hba, antenna_position, tile_config,
+                   element_response);
+        }
+        // Tile enabled in x/y?
+        beam_former_hba->AddTileEnabled(
+            std::array<bool, 2>{!aips_flag(0, i), !aips_flag(0, i)});
+      } else {
+        std::shared_ptr<BeamFormerIdenticalAntennas> beam_former_hba =
+            std::static_pointer_cast<BeamFormerIdenticalAntennas>(beam_former);
+        antenna = MakeTile(id, antenna_position, tile_config, element_response);
+        antenna->enabled_[0] = !aips_flag(0, i);
+        antenna->enabled_[1] = !aips_flag(1, i);
+        beam_former_hba->AddAntenna(antenna);
+      }
     }
-
-    antenna->enabled_[0] = !aips_flag(0, i);
-    antenna->enabled_[1] = !aips_flag(1, i);
-    beam_former->AddAntenna(antenna);
   }
   return beam_former;
 }
@@ -282,7 +353,8 @@ Station::Ptr ReadLofarStation(const MeasurementSet &ms, unsigned int id,
   const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
 
   // Create station.
-  Station::Ptr station(new Station(name, position, model));
+  Station::Ptr station =
+      std::make_shared<Station>(Station(name, position, model));
 
   // Read phase reference position (if available).
   station->SetPhaseReference(ReadStationPhaseReference(ms.antenna(), id));
@@ -299,20 +371,19 @@ Station::Ptr ReadLofarStation(const MeasurementSet &ms, unsigned int id,
     // The Station will consist of a BeamFormer that combines the fields
     // coordinate system is ITRF
     // phase reference is station position
-    auto beam_former = BeamFormer::Ptr(new BeamFormer(
+    auto beam_former = std::make_shared<BeamFormer>(BeamFormer(
         Antenna::IdentityCoordinateSystem, station->GetPhaseReference()));
 
     for (size_t i = 0; i < tab_field.nrow(); ++i) {
       beam_former->AddAntenna(
-          ReadAntennaField(tab_field, i, station->GetElementResponse()));
+          ReadAntennaField(tab_field, i, station->GetElementResponse(),
+                           station->GetElementResponseModel()));
     }
 
     // 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);
-
   } else if (telescope_name == "AARTFAAC") {
     ROScalarColumn<String> ant_type_col(common::GetSubTable(ms, "OBSERVATION"),
                                         "AARTFAAC_ANTENNA_TYPE");

cpp/station.cc

@@ -22,6 +22,7 @@
 
 #include "station.h"
 #include "common/mathutils.h"
+#include "beamformerlofarhba.h"
 
 #include "hamaker/hamakerelementresponse.h"
 #include "oskar/oskarelementresponse.h"
@@ -34,8 +35,9 @@ Station::Station(const std::string &name, const vector3r_t &position,
     : name_(name),
       position_(position),
       phase_reference_(position),
+      element_response_model_(model),
       element_response_(nullptr) {
-  SetResponseModel(model);
+  SetResponseModel(element_response_model_);
   vector3r_t ncp = {{0.0, 0.0, 1.0}};
   ncp_.reset(new coords::ITRFDirection(ncp));
   vector3r_t ncppol0 = {{1.0, 0.0, 0.0}};
@@ -85,13 +87,25 @@ void Station::SetAntenna(Antenna::Ptr 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.
+  // until we have a BeamFormerLofarHBA or 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);
   }
+
+  // If we can cast to BeamFormerLofarHBA, then we extract the Element - please
+  // note that the Element is upcasted from an ElementHamaker into an Element in
+  // BeamFormerLofarHBA::Clone()!- and Transform the Element with the
+  // coordinate_system of the HBA BeamFormer
+  if (auto beamformer_hba =
+          std::dynamic_pointer_cast<BeamFormerLofarHBA>(antenna)) {
+    antenna = beamformer_hba->GetElement();
+    antenna->Transform(beamformer_hba->coordinate_system_);
+  }
   element_ = std::dynamic_pointer_cast<Element>(antenna);
 }
 

cpp/station.h

@@ -42,7 +42,6 @@ class Station {
  public:
   typedef std::shared_ptr<Station> Ptr;
   typedef std::shared_ptr<const Station> ConstPtr;
-  //     typedef std::vector<AntennaField::ConstPtr>  FieldList;
 
   /*!
    *  \brief Construct a new Station instance.
@@ -53,8 +52,6 @@ class Station {
   Station(const std::string &name, const vector3r_t &position,
           const ElementResponseModel model);
 
-  void SetResponseModel(const ElementResponseModel model);
-
   void SetResponse(std::shared_ptr<ElementResponse> element_response);
 
   //! Return the name of the station.
@@ -299,6 +296,11 @@ class Station {
   //! Returns a pointer to the ElementResponse class
   const ElementResponse::Ptr GetElementResponse() { return element_response_; }
 
+  //! Returns an enum of the chosen ElementResponse model
+  const ElementResponseModel GetElementResponseModel() const {
+    return element_response_model_;
+  }
+
   /**
    * @brief Compute the Jones matrix for the element response
    *
@@ -338,6 +340,8 @@ class Station {
   Antenna::Ptr GetAntenna() const { return antenna_; }
 
  private:
+  void SetResponseModel(const ElementResponseModel model);
+
   vector3r_t NCP(real_t time) const;
   vector3r_t NCPPol0(real_t time) const;
   //! Compute the parallactic rotation.
@@ -346,7 +350,8 @@ class Station {
   std::string name_;
   vector3r_t position_;
   vector3r_t phase_reference_;
-  ElementResponseModel element_response_model_ = ElementResponseModel::kUnknown;
+  ElementResponseModel
+      element_response_model_;  // = ElementResponseModel::kUnknown;
   ElementResponse::Ptr element_response_;
   Element::Ptr element_;
 

cpp/test/tlofar_hba.cc

@@ -206,7 +206,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
   // +++ b/lofar/lbeamimagemaker.cpp
   // @@ -380,7 +380,7 @@ void LBeamImageMaker::makeBeamSnapshot(
   //                                 MC4x4::KroneckerProduct(
-  //                                     stationGains[a1].HermTranspose().Transpose(),
+  // stationGains[a1].HermTranspose().Transpose(),
   //                                     stationGains[a2]);
   // -          double w = weights.Value(a1, a2);
   // +          double w = 1.;