Lobes fixes

cpp/antenna.h

@@ -198,6 +198,9 @@ class Antenna {
   vector3r_t phase_reference_position_;
   bool enabled_[2];
 
+ protected:
+  vector3r_t TransformToLocalDirection(const vector3r_t &direction) const;
+
  private:
   virtual matrix22c_t LocalResponse(real_t time, real_t freq,
                                     const vector3r_t &direction,
@@ -208,8 +211,6 @@ class Antenna {
                                      const Options &options) const {
     return {1.0, 1.0};
   }
-
-  vector3r_t TransformToLocalDirection(const vector3r_t &direction) const;
 };
 
 }  // namespace everybeam

cpp/beamformer.cc

@@ -145,6 +145,10 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
     result = result * rotation;
   }
 
+  // Wipe out basefunctions cache
+  if (nullptr != field_response_.get()) {
+    field_response_->ClearFieldQuantities();
+  }
   return result;
 }
 

cpp/element.h

@@ -35,11 +35,43 @@ class Element : public Antenna {
   Antenna::Ptr Clone() const override;
 
   /**
-   * @brief Compute the local response of the element.
+   * @brief Get the Element ID object
+   *
+   * @return size_t
+   */
+  size_t GetElementID() const { return id_; }
+
+  /**
+   * @brief Convenience function to compute the %Element Response a for given
+   * element index
    *
    * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
    * @param freq Frequency of the plane wave (Hz).
    * @param direction Direction of arrival (ITRF, m).
+   * @param id Element index.
+   * @param options
+   * @return matrix22c_t Jones matrix
+   */
+  matrix22c_t ResponseID(real_t time, real_t freq, const vector3r_t &direction,
+                         size_t id, const Options &options = {}) {
+    // Transform direction and directions in options to local coordinatesystem
+    vector3r_t local_direction = TransformToLocalDirection(direction);
+    Options local_options;
+    local_options.freq0 = options.freq0;
+    local_options.station0 = TransformToLocalDirection(options.station0);
+    local_options.tile0 = TransformToLocalDirection(options.tile0);
+    local_options.rotate = options.rotate;
+    local_options.east = TransformToLocalDirection(options.east);
+    local_options.north = TransformToLocalDirection(options.north);
+    return LocalResponse(time, freq, local_direction, id, local_options);
+  }
+
+  /**
+   * @brief Compute the local response of the element.
+   *
+   * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+   * @param freq Frequency of the plane wave (Hz).
+   * @param direction Direction of arrival (East-North-Up, m).
    * @param id ID of element
    * @param options
    * @return matrix22c_t

cpp/fieldresponse.h

@@ -21,6 +21,12 @@ class FieldResponse : public ElementResponse {
    * @param phi Angle in the xy-plane wrt. x-axis  (rad)
    */
   virtual void SetFieldQuantities(double theta, double phi) = 0;
+
+  /**
+   * @brief Clear the cached field quantities
+   *
+   */
+  virtual void ClearFieldQuantities() = 0;
 };
 }  // namespace everybeam
 #endif

cpp/lobes/lobeselementresponse.cc

@@ -142,41 +142,6 @@ LOBESElementResponse::LOBESElementResponse(std::string name) {
 
   nms_.resize(dims_nms[0]);
   dataset.read(nms_.data(), H5::PredType::NATIVE_INT);
-
-  // Normalise coefficients_ with response for an eps from zenith
-  // rather than 0. This is required since F4far_new is singular for theta=0
-  BaseFunctions basefunctions_zenith = ComputeBaseFunctions(1e-6, 0.);
-  Eigen::Index nr_rows = basefunctions_zenith.rows();
-  Eigen::Index nr_freqs = frequencies_.size();
-
-  for (Eigen::Index freq_idx = 0; freq_idx < nr_freqs; ++freq_idx) {
-    std::complex<double> xx = {0}, xy = {0}, yx = {0}, yy = {0};
-    for (Eigen::Index element_id = 0; element_id < coefficients_shape_[2];
-         ++element_id) {
-      for (Eigen::Index i = 0; i < nr_rows; ++i) {
-        std::complex<double> q2 = basefunctions_zenith(i, 0);
-        std::complex<double> q3 = basefunctions_zenith(i, 1);
-        xx += q2 * coefficients_(0, freq_idx, element_id, i);
-        xy += q3 * coefficients_(0, freq_idx, element_id, i);
-        yx += q2 * coefficients_(1, freq_idx, element_id, i);
-        yy += q3 * coefficients_(1, freq_idx, element_id, i);
-      }
-    }
-    // Compute normalisation factor: sqrt(2) / sqrt(1/nr_elements * sum of
-    // response matrix elements ^ 2) double norm_factor =
-    double norm_factor =
-        std::sqrt(2.) /
-        std::sqrt(1. / double(coefficients_shape_[2]) *
-                  (std::pow(std::abs(xx), 2) + std::pow(std::abs(xy), 2) +
-                   std::pow(std::abs(yx), 2) + std::pow(std::abs(yy), 2)));
-    for (Eigen::Index element_id = 0; element_id < coefficients_shape_[2];
-         ++element_id) {
-      for (Eigen::Index i = 0; i < nr_rows; ++i) {
-        coefficients_(0, freq_idx, element_id, i) *= norm_factor;
-        coefficients_(1, freq_idx, element_id, i) *= norm_factor;
-      }
-    }
-  }
 }
 
 LOBESElementResponse::BaseFunctions LOBESElementResponse::ComputeBaseFunctions(
@@ -205,6 +170,25 @@ LOBESElementResponse::BaseFunctions LOBESElementResponse::ComputeBaseFunctions(
 void LOBESElementResponse::Response(
     int element_id, double freq, double theta, double phi,
     std::complex<double> (&response)[2][2]) const {
+  // Initialize the response to zero.
+  response[0][0] = 0.0;
+  response[0][1] = 0.0;
+  response[1][0] = 0.0;
+  response[1][1] = 0.0;
+
+  // Clip directions below the horizon.
+  if (theta >= M_PI_2) {
+    return;
+  }
+
+  // Fill basefunctions if not yet set. Set clear_basefunctions to false
+  // to disable the caching of the basefunctions
+  bool clear_basefunctions = false;
+  if (basefunctions_.rows() == 0) {
+    clear_basefunctions = true;
+    basefunctions_ = ComputeBaseFunctions(theta, phi);
+  }
+
   int freq_idx = FindFrequencyIdx(freq);
   std::complex<double> xx = {0}, xy = {0}, yx = {0}, yy = {0};
 
@@ -228,6 +212,11 @@ void LOBESElementResponse::Response(
   response[0][1] = xy;
   response[1][0] = yx;
   response[1][1] = yy;
+
+  if (clear_basefunctions) {
+    // Do a destructive resize
+    basefunctions_.resize(0, 2);
+  }
 }
 
 std::shared_ptr<LOBESElementResponse> LOBESElementResponse::GetInstance(

cpp/lobes/lobeselementresponse.h

@@ -67,10 +67,19 @@ class LOBESElementResponse : public FieldResponse {
     basefunctions_ = ComputeBaseFunctions(theta, phi);
   };
 
+  /**
+   * @brief Clear the cached basefunctions
+   *
+   */
+  virtual void ClearFieldQuantities() final override {
+    // Destructively resize the basefunctions_ to 0 rows
+    basefunctions_.resize(0, 2);
+  };
+
  private:
   // Typdef of BaseFunctions as Eigen::Array type
   typedef Eigen::Array<std::complex<double>, Eigen::Dynamic, 2> BaseFunctions;
-  BaseFunctions basefunctions_;
+  mutable BaseFunctions basefunctions_;
 
   // Compose the path to a LOBES coefficient file
   std::string GetPath(const char *) const;
@@ -87,7 +96,7 @@ class LOBESElementResponse : public FieldResponse {
   BaseFunctions ComputeBaseFunctions(double theta, double phi) const;
 
   // Store h5 coefficients in coefficients_
-  Eigen::Tensor<std::complex<double>, 4> coefficients_;
+  Eigen::Tensor<std::complex<double>, 4, Eigen::RowMajor> coefficients_;
   std::vector<unsigned int> coefficients_shape_;
 
   // Store h5 frequencies in frequencies_

cpp/station.cc

@@ -121,8 +121,8 @@ void Station::SetAntenna(Antenna::Ptr antenna) {
 // ========================================================
 matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
                                             const vector3r_t &direction,
-                                            size_t id,
-                                            const bool rotate) const {
+                                            size_t id, bool is_local,
+                                            bool rotate) const {
   Antenna::Options options;
   options.rotate = rotate;
 
@@ -134,12 +134,13 @@ matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
     options.north = north;
   }
 
-  return element_->LocalResponse(time, freq, direction, id, options);
+  return is_local ? element_->LocalResponse(time, freq, direction, id, options)
+                  : element_->ResponseID(time, freq, direction, id, options);
 }
 
 matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
                                             const vector3r_t &direction,
-                                            const bool rotate) const {
+                                            bool is_local, bool rotate) const {
   Antenna::Options options;
   options.rotate = rotate;
 
@@ -151,9 +152,9 @@ matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
     options.north = north;
   }
 
-  matrix22c_t response;
-  response = element_->Response(time, freq, direction, options);
-  return response;
+  return is_local ? element_->LocalResponse(time, freq, direction,
+                                            element_->GetElementID(), options)
+                  : element_->Response(time, freq, direction, options);
 }
 
 matrix22c_t Station::Response(real_t time, real_t freq,

cpp/station.h

@@ -306,27 +306,35 @@ class Station {
    *
    * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
    * @param freq Frequency of the plane wave (Hz).
-   * @param direction Direction of arrival (ITRF, m).
+   * @param direction Direction of arrival. If is_local is true: (ENU, m) else
+   * direction vector in global coord system is assumed.
+   * @param is_local Use local east-north-up system (true) or global coordinate
+   * system (false, default).
    * @param id Element id
    * @param rotate Boolean deciding if paralactic rotation should be applied.
    * @return matrix22c_t Jones matrix of element response
    */
   matrix22c_t ComputeElementResponse(real_t time, real_t freq,
                                      const vector3r_t &direction, size_t id,
-                                     const bool rotate) const;
+                                     bool is_local = false,
+                                     bool rotate = true) const;
 
   /**
    * @brief Compute the Jones matrix for the element response
    *
    * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
    * @param freq Frequency of the plane wave (Hz).
-   * @param direction Direction of arrival (ITRF, m).
+   * @param direction Direction of arrival. If is_local is true: (ENU, m) else
+   * direction vector in global coord system is assumed.
+   * @param is_local Use local east-north-up system (true) or global coordinate
+   * system (false, default).
    * @param rotate Boolean deciding if paralactic rotation should be applied.
    * @return matrix22c_t Jones matrix of element response
    */
   matrix22c_t ComputeElementResponse(real_t time, real_t freq,
                                      const vector3r_t &direction,
-                                     const bool rotate = true) const;
+                                     bool is_local = false,
+                                     bool rotate = true) const;
 
   //! Specialized implementation of response function.
   matrix22c_t Response(real_t time, real_t freq,

cpp/test/tlofar_hba.cc

@@ -77,7 +77,7 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
   const Station& station =
       static_cast<const Station&>(*(lofartelescope.GetStation(11).get()));
   matrix22c_t element_response =
-      station.ComputeElementResponse(time, frequency, direction, true);
+      station.ComputeElementResponse(time, frequency, direction, false);
 
   // Check whether element_response and target_element_response are "equal"
   for (size_t i = 0; i != 2; ++i) {

cpp/test/tlofar_lba.cc

@@ -73,7 +73,7 @@ BOOST_AUTO_TEST_CASE(test_hamaker) {
   const Station& station =
       static_cast<const Station&>(*(lofartelescope.GetStation(19).get()));
   matrix22c_t element_response =
-      station.ComputeElementResponse(time, frequency, direction, true);
+      station.ComputeElementResponse(time, frequency, direction, false);
 
   for (size_t i = 0; i != 2; ++i) {
     for (size_t j = 0; j != 2; ++j) {
@@ -163,15 +163,16 @@ BOOST_AUTO_TEST_CASE(test_lobes) {
   // Compare with everybeam at pixel (1, 3). This solution only is a "reference"
   // certainly not a "ground-truth"
   std::vector<std::complex<float>> everybeam_ref_p13 = {
-      {0.022321859, 0.032528818},
-      {-0.0022747002, 0.096434057},
-      {0.037618704, 0.015968619},
-      {-0.041608963, 0.021379814}};
-  std::size_t offset_13 = (3 + 1 * width) * 4;
+      {-2.1952772, -2.7803752},
+      {-1.4990979, -1.9227437},
+      {0.0086051673, 0.034685627},
+      {1.5801408, 2.0045171}};
 
+  std::size_t offset_13 = (3 + 1 * width) * 4;
   for (std::size_t i = 0; i < 4; ++i) {
+    // relative tolerance is 2e-1%, so 2e-3
     BOOST_CHECK_CLOSE(antenna_buffer_single[offset_13 + i],
-                      everybeam_ref_p13[i], 1e-4);
+                      everybeam_ref_p13[i], 2e-1);
   }
   // const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
   // npy::SaveArrayAsNumpy("lobes_station_response.npy", false, 4, leshape,

demo/lobes/lobes_image.py

+from everybeam import load_telescope, LOFAR, Options, thetaphi2cart
+import matplotlib.pyplot as plt
+import numpy as np
+import os
+
+# MS
+try:
+    datadir = os.environ["DATA_DIR"]
+except:
+    ValueError("DATA_DIR variable not found in environment variables")
+filename = "LOFAR_LBA_MOCK.ms"
+ms_path = os.path.join(datadir, filename)
+
+# Response settings
+mode = "element"
+time = 4.92183348e09
+frequency = 57812500.0
+station_id = 20
+element_id = 0
+
+# Same station0 direction as in python/test
+station0 = np.array([0.655743, -0.0670973, 0.751996])
+is_local = True  # use local coords
+rotate = True
+response_model = "lobes"
+
+# Load telescope
+telescope = load_telescope(ms_path, element_response_model=response_model)
+station_name = telescope.station_name(20)
+assert station_name == "CS302LBA"
+
+# Make coordinates for mesh
+x_v = np.linspace(-1.0, 1.0, 128)
+y_v = np.linspace(-1.0, 1.0, 128)
+
+response = np.empty((y_v.size, x_v.size, 2, 2), dtype=np.cdouble)
+response.fill(np.nan)
+for i, x in enumerate(x_v):
+    for j, y in enumerate(y_v):
+        if (x ** 2 + y ** 2) <= 1.0:
+            # Compute theta/phi and resulting direction vector
+            theta = np.arcsin(np.sqrt(x * x + y * y))
+            phi = np.arctan2(y, x)
+            direction = thetaphi2cart(theta, phi)
+            if mode == "element":
+                response[j, i, :, :] = telescope.element_response(
+                    time,
+                    station_id,
+                    element_id,
+                    frequency,
+                    direction,
+                    is_local,
+                    rotate=rotate,
+                )
+            elif mode == "station":
+                response[j, i, :, :] = telescope.station_response(
+                    time, station_id, frequency, direction, station0, rotate=rotate
+                )
+            else:
+                raise Exception(
+                    "Unrecognized response mode. Must be either station or element"
+                )
+
+# Plotting
+X, Y = np.meshgrid(x_v, y_v)
+label = ("X", "Y")
+fig, axs = plt.subplots(2, 4, figsize=(10, 5), sharex=True, sharey=True)
+fig.suptitle(f"{mode} response for station {station_name}: {response_model}")
+for i, row in enumerate(range(2)):
+    for j, col in enumerate(range(2)):
+        im1 = axs[row, col].pcolor(X, Y, np.abs(response[:, :, row, col]))
+        axs[row, col].set_aspect("equal", "box")
+        axs[row, col].set_title(f"abs({label[j]}{label[i]})")
+        fig.colorbar(im1, ax=axs[row, col])
+
+        im2 = axs[row, col + 2].pcolor(X, Y, np.angle(response[:, :, row, col]))
+        axs[row, col + 2].set_aspect("equal", "box")
+        axs[row, col + 2].set_title(f"angle({label[j]}{label[i]})")
+        fig.colorbar(im2, ax=axs[row, col + 2])
+
+plt.show()

python/wrappers/pytelescope.cc

@@ -139,6 +139,12 @@ void init_telescope(py::module &m) {
         -------
         int
        )pbdoc")
+      .def("station_name",
+           [](PhasedArray &self, size_t idx) {
+             const Station &station =
+                 static_cast<const Station &>(*(self.GetStation(idx).get()));
+             return station.GetName();
+           })
       .def("channel_frequency", &PhasedArray::GetChannelFrequency,
            R"pbdoc(
         Retrieve channel frequency for a given (zero-based) channel index.
@@ -460,7 +466,7 @@ void init_telescope(py::module &m) {
             Evaluation response at time.
             Time in modified Julian date, UTC, in seconds (MJD(UTC), s)
         station_idx: int
-            Get response for station index
+            station index
         freq: float
             Frequency of the plane wave (Hz)
         direction: np.1darray
@@ -477,6 +483,54 @@ void init_telescope(py::module &m) {
        )pbdoc",
           py::arg("time"), py::arg("station_idx"), py::arg("freq"),
           py::arg("direction"), py::arg("station0_direction"),
+          py::arg("rotate") = true)
+      .def(
+          "element_response",
+          [](PhasedArray &self, double time, size_t station_idx,
+             size_t element_idx, double freq,
+             const py::array_t<double> pydirection, bool is_local,
+             bool rotate) -> py::array_t<std::complex<double>> {
+            if (station_idx >= self.GetNrStations()) {
+              throw std::runtime_error(
+                  "Requested station index exceeds number of stations.");
+            }
+            // TODO: need such a check for element index too
+            vector3r_t direction = np2vector3r_t(pydirection);
+            const Station &station = static_cast<const Station &>(
+                *(self.GetStation(station_idx).get()));
+            matrix22c_t response = station.ComputeElementResponse(
+                time, freq, direction, element_idx, is_local, rotate);
+            return py::array_t<std::complex<double>>{cast_matrix(response)};
+          },
+          R"pbdoc(
+        Get element response given a station and an element in prescribed direction.
+
+        Parameters
+        ----------
+        time: double
+            Evaluation response at time. 
+            Time in modified Julian date, UTC, in seconds (MJD(UTC), s) 
+        station_idx: int
+            station index
+        element_idx: int
+            element index
+        freq: float
+            Frequency of the plane wave (Hz)
+        direction: np.1darray
+            Direction of arrival either in ITRF (m) or local East-North-Up (m)
+        is_local: bool, optional
+            Is the specified direction in local East-North-Up? If not, global coordinate 
+            system is assumed. [True/False] Defaults to False. 
+        rotate: bool, optional
+            Apply paralactic angle rotation? [True/False] Defaults to True
+
+        Returns
+        -------
+        np.ndarray
+            Response (Jones) matrix
+       )pbdoc",
+          py::arg("time"), py::arg("station_idx"), py::arg("element_idx"),
+          py::arg("freq"), py::arg("direction"), py::arg("is_local") = false,
           py::arg("rotate") = true);
 
   py::class_<LOFAR, PhasedArray>(m, "LOFAR").def(py::init(&create_lofar));

python/wrappers/pyutils.cc

 #include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <pybind11/stl.h>
 
 #include "elementresponse.h"
 #include "options.h"
+#include "common/mathutils.h"
 #include "coords/coordutils.h"
 
 namespace py = pybind11;
+using everybeam::cart2thetaphi;
 using everybeam::ElementResponseModel;
 using everybeam::Options;
+using everybeam::thetaphi2cart;
+using everybeam::vector2r_t;
+using everybeam::vector3r_t;
 using everybeam::coords::CoordinateSystem;
 
+namespace {
+// Convert pyarray of size 3 to vector3r_t
+vector3r_t np2vector3r_t(const py::array_t<double> pyarray) {
+  auto r = pyarray.unchecked<1>();
+  if (r.size() != 3) {
+    throw std::runtime_error("Pyarry is of incorrect size, must be 3.");
+  }
+  return vector3r_t{r[0], r[1], r[2]};
+}
+}  // namespace
+
 void init_utils(py::module &m) {
+  m.def(
+      "cart2thetaphi",
+      [](py::array_t<double> pydirection) {
+        vector3r_t direction = np2vector3r_t(pydirection);
+        vector2r_t thetaphi = cart2thetaphi(direction);
+        return thetaphi;
+      },
+      R"pbdoc(
+        Convert direction vector (ITRF or local East-North-Up)
+        to theta, phi.
+        
+        Parameters
+        ---------- 
+        direction: np.1darray
+            Direction vector
+        
+        Returns
+        -------
+        list:
+            [theta, phi]   
+       )pbdoc",
+      py::arg("direction"));
+  m.def(
+      "thetaphi2cart",
+      [](double theta, double phi) -> py::array_t<float> {
+        vector2r_t thetaphi = {theta, phi};
+        vector3r_t direction = everybeam::thetaphi2cart(thetaphi);
+        return py::array_t<float>(py::cast(direction));
+      },
+      R"pbdoc(
+        Convert theta, phi angles to direction vector 
+        
+        Parameters
+        ---------- 
+        theta: float
+            theta angle [rad]
+        phi: float
+            phi angle [rad]
+        
+        Returns
+        -------
+        np.1darray:
+            
+       )pbdoc",
+      py::arg("theta"), py::arg("phi"));
+
   // Bindings for CoordinateSystem struct
   py::class_<CoordinateSystem>(m, "CoordinateSystem",
                                R"pbdoc(