diff --git a/.gitignore b/.gitignore
index 2e5a2d7b5826bbbf2e186b27300b52f5d84f0655..a5dd1f20299365a163d2dc6cabf38e997f23395a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,7 +2,12 @@
 #
 # Directories
 .vscode/
+coeffs/lobes
 build/
 builds/
 test_data/
 __pycache__/
+cpp/test/thamaker_tmp.cc
+cpp/test/tlobes_tmp.cc
+download_lofar_lobes_ms.sh
+
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 3c58e752d2389a6525479d33787c9ea30e6017a9..f0e8fe282f1d4f6e3a9033b70ca2d5629afa26b1 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -22,11 +22,11 @@ clang-format:
 
 # Build a debug version of EveryBeam from the base image
 test-and-coverage:
-  except:
-    - schedules
+  # stage is always triggered
   stage: build
   image: everybeam_base:${CI_COMMIT_SHORT_SHA}
   before_script:
+    - apt-get update
     - apt-get -y install python3-pip
     - python3 -m pip install gcovr
   script:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2cd585365062893574c54aad112104fac5e4454b..0152888ead344fa6d663ceac58ebf0a30f6dc6c3 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -116,6 +116,12 @@ endif()
 # Generate config.h file
 configure_file(${CMAKE_SOURCE_DIR}/CMake/config.h.in ${CMAKE_BINARY_DIR}/config.h)
 
+#------------------------------------------------------------------------------
+# pybind11 
+# NOTE: temporary fix to build pybindings for lobes
+add_subdirectory(python/lobes)
+
+
 #------------------------------------------------------------------------------
 # Documentation
 find_package(Doxygen)
diff --git a/cpp/beamformer.cc b/cpp/beamformer.cc
index 43634b4fda0230f4b02f652540b9e7dffbab69ec..e1d98fd441759be1c33b3f45136bd18e117c8960 100644
--- a/cpp/beamformer.cc
+++ b/cpp/beamformer.cc
@@ -98,6 +98,20 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
   std::vector<std::complex<double>> geometric_response =
       ComputeGeometricResponse(freq, direction);
 
+  // Copy options into local_options. Needed to propagate
+  // the potential change in the rotate boolean downstream
+  Options local_options = options;
+  // If field_response_ not nullptr, set/precompute quantities
+  // related to the field
+  if (nullptr != field_response_.get()) {
+    // Lock the mutex, to avoid that the LOBESElementResponse basefunctions_
+    // are overwritten before response is computed
+    mtx_.lock();
+    vector2r_t thetaphi = cart2thetaphi(direction);
+    field_response_->SetFieldQuantities(thetaphi[0], thetaphi[1]);
+    local_options.rotate = false;
+  }
+
   matrix22c_t result = {0};
   for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size();
        ++antenna_idx) {
@@ -108,7 +122,7 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
         geometric_response[antenna_idx];
 
     matrix22c_t antenna_response =
-        antenna->Response(time, freq, direction, options);
+        antenna->Response(time, freq, direction, local_options);
     result[0][0] += antenna_weight.first * antenna_geometric_reponse *
                     antenna_response[0][0];
     result[0][1] += antenna_weight.first * antenna_geometric_reponse *
@@ -118,6 +132,23 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq,
     result[1][1] += antenna_weight.second * antenna_geometric_reponse *
                     antenna_response[1][1];
   }
+
+  // If the Jones matrix needs to be rotated from theta, phi directions
+  // to north, east directions, but this has not been done yet, do it here
+  if (options.rotate && !local_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;
+  }
+
+  // Unlock mutex in case it was locked
+  if (nullptr != field_response_.get()) {
+    mtx_.unlock();
+  }
   return result;
 }
 
diff --git a/cpp/beamformer.h b/cpp/beamformer.h
index b073d68ab852e3895d2a141a33ad51bf3fd13222..245bdac2d0b07ddc4ac450a57362e880cd7b0644 100644
--- a/cpp/beamformer.h
+++ b/cpp/beamformer.h
@@ -6,6 +6,9 @@
 
 #include "element.h"
 #include "common/types.h"
+#include "fieldresponse.h"
+
+#include <mutex>
 
 namespace everybeam {
 class BeamFormer : public Antenna {
@@ -16,7 +19,8 @@ class BeamFormer : public Antenna {
    * @brief Construct a new BeamFormer object
    *
    */
-  BeamFormer() : Antenna() {
+  BeamFormer(std::shared_ptr<FieldResponse> field_response = nullptr)
+      : Antenna(), field_response_(field_response) {
     local_phase_reference_position_ =
         TransformToLocalPosition(phase_reference_position_);
   }
@@ -26,8 +30,9 @@ class BeamFormer : public Antenna {
    *
    * @param coordinate_system
    */
-  BeamFormer(const CoordinateSystem &coordinate_system)
-      : Antenna(coordinate_system) {
+  BeamFormer(const CoordinateSystem &coordinate_system,
+             std::shared_ptr<FieldResponse> field_response = nullptr)
+      : Antenna(coordinate_system), field_response_(field_response) {
     local_phase_reference_position_ =
         TransformToLocalPosition(phase_reference_position_);
   }
@@ -40,14 +45,17 @@ class BeamFormer : public Antenna {
    * @param phase_reference_position
    */
   BeamFormer(CoordinateSystem coordinate_system,
-             vector3r_t phase_reference_position)
-      : Antenna(coordinate_system, phase_reference_position) {
+             vector3r_t phase_reference_position,
+             std::shared_ptr<FieldResponse> field_response = nullptr)
+      : Antenna(coordinate_system, phase_reference_position),
+        field_response_(field_response) {
     local_phase_reference_position_ =
         TransformToLocalPosition(phase_reference_position_);
   }
 
-  BeamFormer(vector3r_t phase_reference_position)
-      : Antenna(phase_reference_position) {
+  BeamFormer(vector3r_t phase_reference_position,
+             std::shared_ptr<FieldResponse> field_response = nullptr)
+      : Antenna(phase_reference_position), field_response_(field_response) {
     local_phase_reference_position_ =
         TransformToLocalPosition(phase_reference_position_);
   }
@@ -111,6 +119,12 @@ class BeamFormer : public Antenna {
 
   // List of antennas in BeamFormer
   std::vector<Antenna::Ptr> antennas_;
+
+  // (Optional) shared pointer to field response model (e.g. LOBES)
+  std::shared_ptr<FieldResponse> field_response_;
+
+ private:
+  mutable std::mutex mtx_;
 };
 }  // namespace everybeam
 #endif
diff --git a/cpp/fieldresponse.h b/cpp/fieldresponse.h
new file mode 100644
index 0000000000000000000000000000000000000000..c1642997ff08e576d829608714b186244b3606a5
--- /dev/null
+++ b/cpp/fieldresponse.h
@@ -0,0 +1,26 @@
+#ifndef EVERYBEAM_FIELDRESPONSE_H
+#define EVERYBEAM_FIELDRESPONSE_H
+
+#include "elementresponse.h"
+
+namespace everybeam {
+
+/**
+ * @brief Virtual class that inherits from the ElementResponse class but can be
+ * used to precompute quantities at "field" level. Currently, the
+ * LOBESElementResponse class inherits from this class
+ *
+ */
+class FieldResponse : public ElementResponse {
+ public:
+  /**
+   * @brief Virtual method that can be used to set field quantities based
+   * on the direction of interest, specified in (theta, phi) angles
+   *
+   * @param theta Angle wrt. z-axis (rad)
+   * @param phi Angle in the xy-plane wrt. x-axis  (rad)
+   */
+  virtual void SetFieldQuantities(double theta, double phi) = 0;
+};
+}  // namespace everybeam
+#endif
diff --git a/cpp/griddedresponse/phasedarraygrid.cc b/cpp/griddedresponse/phasedarraygrid.cc
index f33fc8592ea3c95f22aa5d2258e98f0598481323..3d8384d7565fedf39df9c7753741abe379392dce 100644
--- a/cpp/griddedresponse/phasedarraygrid.cc
+++ b/cpp/griddedresponse/phasedarraygrid.cc
@@ -17,6 +17,7 @@ PhasedArrayGrid::PhasedArrayGrid(
   // Set private members
   use_differential_beam_ = telescope_->GetOptions().use_differential_beam;
   size_t ncpus = aocommon::ThreadPool::NCPUs();
+  // TODO: toskar fails on DAS5
   nthreads_ = std::min(ncpus, telescope_->GetNrStations());
   threads_.resize(nthreads_);
 };
diff --git a/cpp/lobes/CMakeLists.txt b/cpp/lobes/CMakeLists.txt
index ab07c8bb1e3794ad52521f5b137f9823b9dc178a..c1596f25a5c0f567e40fac81de0ca1859c604e39 100644
--- a/cpp/lobes/CMakeLists.txt
+++ b/cpp/lobes/CMakeLists.txt
@@ -1,13 +1,16 @@
 #------------------------------------------------------------------------------
-# Required packages for lobes
-# pybind and eigen not required here, will be moved to python dir with pybindings
-# find_package(pybind11 REQUIRED)
-# find_package (Eigen3 REQUIRED NO_MODULE)
+# directory for config.h
+include_directories(${CMAKE_BINARY_DIR})
 
-# pybind11_add_module(pylobes SHARED lobes.cc ElementResponse.cc LOBESElementResponse.cc)
-# target_link_libraries (pylobes PUBLIC Eigen3::Eigen ${HDF5_LIBRARIES})
+#------------------------------------------------------------------------------
+# Download the LOBES coefficient files, leads to a relatively small amount of overhead
+# in make command, even if the coefficient files have been dowloaded
+add_custom_target(download_lobes_coefficients COMMAND ${CMAKE_SOURCE_DIR}/scripts/download_lobes_coeffs.sh)
 
+#------------------------------------------------------------------------------
 add_library(lobes SHARED lobeselementresponse.cc)
+add_dependencies(lobes download_lobes_coefficients)
+target_link_libraries(lobes PUBLIC ${HDF5_LIBRARIES})
 string(TOLOWER ${CMAKE_PROJECT_NAME} projectname )
 set_target_properties(lobes PROPERTIES LIBRARY_OUTPUT_NAME "${projectname}-lobes")
 
@@ -18,6 +21,10 @@ install(
   EXPORT EveryBeamTargets
   DESTINATION lib)
 
+# install coefficients
+message("install lobes coefficients in: " ${CMAKE_INSTALL_DATA_DIR}/lobes)
+install(DIRECTORY "${CMAKE_SOURCE_DIR}/coeffs/lobes" DESTINATION ${CMAKE_INSTALL_DATA_DIR} FILES_MATCHING PATTERN "LOBES_*")
+
 #------------------------------------------------------------------------------
 # TODO: can't we remove all this?
 set(MISC_DIR ${CMAKE_SOURCE_DIR}/scripts/misc)
@@ -29,6 +36,3 @@ configure_file(${MISC_DIR}/hamaker_vs_lobes.py ${CMAKE_CURRENT_BINARY_DIR}/hamak
 configure_file(${MISC_DIR}/test_beam_model.py ${CMAKE_CURRENT_BINARY_DIR}/test_beam_model.py COPYONLY)
 
 configure_file(${CMAKE_SOURCE_DIR}/coeffs/CS302_coords.mat ${CMAKE_CURRENT_BINARY_DIR}/CS302_coords.mat COPYONLY)
-
-# TODO Too large for git repo, get file from somewhere else
-# configure_file(LBA_CS302_fine.mat LBA_CS302_fine.mat COPYONLY)
diff --git a/cpp/lobes/ElementResponse.cc b/cpp/lobes/ElementResponse.cc
deleted file mode 100644
index 5c68d667122785111de9398db0019196a7d928a2..0000000000000000000000000000000000000000
--- a/cpp/lobes/ElementResponse.cc
+++ /dev/null
@@ -1,146 +0,0 @@
-// ElementResponse.cc: Functions to compute the (idealized) response of a LOFAR
-// LBA or HBA dual dipole antenna.
-//
-// Copyright (C) 2011
-// 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 "ElementResponse.h"
-
-#include <cmath>
-
-// The coefficients are kept in an unnamed namespace which effectively makes
-// them invisible outside this translation unit.
-namespace {
-#include "defaultcoefflba.cc"
-#include "defaultcoeffhba.cc"
-}  // namespace
-
-namespace everybeam {
-
-void ElementResponseLBA(double freq, double theta, double phi,
-                        std::complex<double> (&response)[2][2]) {
-  ElementResponse(freq, theta, phi, response, default_lba_freq_center,
-                  default_lba_freq_range, default_lba_coeff_shape,
-                  default_lba_coeff);
-}
-
-void ElementResponseHBA(double freq, double theta, double phi,
-                        std::complex<double> (&response)[2][2]) {
-  ElementResponse(freq, theta, phi, response, default_hba_freq_center,
-                  default_hba_freq_range, default_hba_coeff_shape,
-                  default_hba_coeff);
-}
-
-void ElementResponse(double freq, double theta, double phi,
-                     std::complex<double> (&response)[2][2], double freq_center,
-                     double freq_range, const unsigned int (&coeff_shape)[3],
-                     const std::complex<double> coeff[]) {
-  // 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;
-  }
-
-  const unsigned int nHarmonics = coeff_shape[0];
-  const unsigned int nPowerTheta = coeff_shape[1];
-  const unsigned int nPowerFreq = coeff_shape[2];
-
-  // The model is parameterized in terms of a normalized frequency in the
-  // range [-1, 1]. The appropriate conversion is taken care of below.
-  freq = (freq - freq_center) / freq_range;
-
-  // The variables sign and kappa are used to compute the value of kappa
-  // mentioned in the description of the beam model [kappa = (-1)^k * (2 * k
-  //+ 1)] incrementally.
-  int sign = 1, kappa = 1;
-
-  std::complex<double> P[2], Pj[2];
-  for (unsigned int k = 0; k < nHarmonics; ++k) {
-    // Compute the (diagonal) projection matrix P for the current harmonic.
-    // This requires the evaluation of two polynomials in theta and freq (of
-    // degree nPowerTheta in theta and nPowerFreq in freq), one for each
-    // element of P. The polynomials are evaluated using Horner's rule.
-
-    // Horner's rule requires backward iteration of the coefficients, so
-    // move the iterator to the first location past the end of the block of
-    // coefficients for the current harmonic (k).
-    coeff += nPowerTheta * nPowerFreq * 2;
-
-    // Evaluate the polynomial. Note that the iterator is always decremented
-    // before it is dereferenced, using the prefix operator--. After
-    // evaluation of the polynomial, the iterator points exactly to the
-    // beginning of the block of coefficients for the current harmonic (k),
-    // that is, all the decrements together exactly cancel the increment
-    // aplied above.
-
-    // Evaluate the highest order term. Note that the order of the
-    // assigments is important because of the iterator decrement, i.e. P[1]
-    // should be assigned first.
-    P[1] = *--coeff;
-    P[0] = *--coeff;
-
-    for (unsigned int i = 0; i < nPowerFreq - 1; ++i) {
-      P[1] = P[1] * freq + *--coeff;
-      P[0] = P[0] * freq + *--coeff;
-    }
-
-    // Evaluate the remaining terms.
-    for (unsigned int j = 0; j < nPowerTheta - 1; ++j) {
-      Pj[1] = *--coeff;
-      Pj[0] = *--coeff;
-
-      for (unsigned int i = 0; i < nPowerFreq - 1; ++i) {
-        Pj[1] = Pj[1] * freq + *--coeff;
-        Pj[0] = Pj[0] * freq + *--coeff;
-      }
-
-      P[1] = P[1] * theta + Pj[1];
-      P[0] = P[0] * theta + Pj[0];
-    }
-
-    // Because the increment and decrements cancel, the iterator points to
-    // the same location as at the beginning of this iteration of the outer
-    // loop. The next iteration should use the coefficients for the next
-    // harmonic (k), so we move the iterator to the start of that block.
-    coeff += nPowerTheta * nPowerFreq * 2;
-
-    // Compute the Jones matrix for the current harmonic, by rotating P over
-    // kappa * az, and add it to the result.
-    const double angle = sign * kappa * phi;
-    const double caz = std::cos(angle);
-    const double saz = std::sin(angle);
-
-    response[0][0] += caz * P[0];
-    response[0][1] += -saz * P[1];
-    response[1][0] += saz * P[0];
-    response[1][1] += caz * P[1];
-
-    // Update sign and kappa.
-    sign = -sign;
-    kappa += 2;
-  }
-}
-
-}  // namespace everybeam
diff --git a/cpp/lobes/ElementResponse.h b/cpp/lobes/ElementResponse.h
deleted file mode 100644
index 5fa7b98af8d6d480bfce29731dffeb3eba91b69f..0000000000000000000000000000000000000000
--- a/cpp/lobes/ElementResponse.h
+++ /dev/null
@@ -1,100 +0,0 @@
-// ElementResponse.h: Functions to compute the (idealized) response of a LOFAR
-// LBA or HBA dual dipole antenna.
-//
-// Copyright (C) 2011
-// 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 LOFAR_ELEMENTRESPONSE_H
-#define LOFAR_ELEMENTRESPONSE_H
-
-// \file
-// Functions to compute the (idealized) response of a LOFAR LBA or HBA dual
-// dipole antenna.
-
-#include <complex>
-
-namespace everybeam {
-
-// \addtogroup ElementResponse
-// @{
-
-// Compute the response of an idealized LOFAR LBA dual dipole antenna to
-// radiation at frequency freq (Hz) arriving from the direction given by theta,
-// phi (rad). The antenna model is described in a spherical coordinate system
-// with coordinates theta (zenith angle) and phi (azimith). The +X dipole is at
-// azimuth zero, the +Y dipole is at azimuth PI / 2.0.
-//
-// Preconditions:
-// --------------
-// freq: Frequency in Hz in the range [10 MHz, 100 MHz].
-// theta: Zenith angle in rad in the range [0.0, PI / 2.0].
-// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
-//
-void ElementResponseLBA(double freq, double theta, double phi,
-                        std::complex<double> (&response)[2][2]);
-
-// Compute the response of an idealized LOFAR HBA dual dipole antenna to
-// radiation at frequency freq (Hz) arriving from the direction given by theta,
-// phi (rad). The antenna model is described in a spherical coordinate system
-// with coordinates theta (zenith angle) and phi (azimith). The +X dipole is at
-// azimuth zero, the +Y dipole is at azimuth PI / 2.0.
-//
-// Preconditions:
-// --------------
-// freq: Frequency in Hz in the range [120 MHz, 240 MHz].
-// theta: Zenith angle in rad in the range [0.0, PI / 2.0].
-// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
-//
-void ElementResponseHBA(double freq, double theta, double phi,
-                        std::complex<double> (&response)[2][2]);
-
-// Compute the response of an idealized LOFAR dual dipole antenna to radiation
-// at frequency freq (Hz) arriving from the direction given by theta, phi (rad).
-// The antenna model is described in a spherical coordinate system with
-// coordinates theta (zenith angle) and phi (azimith). The +X dipole is at
-// azimuth zero, the +Y dipole is at azimuth PI / 2.0.
-//
-// This function uses a set of user defined coefficients to evaluate the beam
-// model. The coeff_shape parameter defines the shape of the coefficient array
-// as no. of harmonics x degree in theta x degree in frequency x 2. The last
-// dimension is implicit and always equal to 2. The coeff parameter points to an
-// array of coefficients of the proper size, stored in row-major order
-// ("C"-order). The freq_center and freq_range parameters define the frequency
-// range covered by the model described by the set of coefficients.
-//
-// Preconditions:
-// --------------
-// freq: Frequency in Hz in the range [freq_center - freq_range,
-//     freq_center + freq_range].
-// theta: Zenith angle in rad in the range [0.0, PI / 2.0].
-// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
-// freq_range, freq_center: Frequency center and range in Hz, should be > 0.
-// coeff_shape: Shape of the coefficient array, all dimensions should be > 0.
-//
-void ElementResponse(double freq, double theta, double phi,
-                     std::complex<double> (&response)[2][2], double freq_center,
-                     double freq_range, const unsigned int (&coeff_shape)[3],
-                     const std::complex<double> coeff[]);
-
-// @}
-
-}  // namespace everybeam
-
-#endif
diff --git a/cpp/lobes/defaultcoeffhba.cc b/cpp/lobes/defaultcoeffhba.cc
deleted file mode 100644
index 7d660b92338dd0e4153ce9f98432edf7f7a62210..0000000000000000000000000000000000000000
--- a/cpp/lobes/defaultcoeffhba.cc
+++ /dev/null
@@ -1,123 +0,0 @@
-// Beam model coefficients converted by convert_coeff.py.
-// Conversion performed on 2011/09/14/08:23:16 UTC using:
-//     convert_coeff.py element_beam_HBA.coeff defaultcoeffhba.cc default_hba
-
-#include <complex>
-
-// Center frequency, and frequency range for which the beam model coefficients
-// are valid. The beam model is parameterized in terms of a normalized
-// frequency f' in the range [-1.0, 1.0]. The appropriate conversion is:
-//
-//     f' = (f - center) / range
-//
-const double default_hba_freq_center = 180e6;
-const double default_hba_freq_range = 60e6;
-
-// Shape of the coefficient array: 2x5x5x2 (the size of the last dimension is
-// implied, and always equal to 2).
-//
-const unsigned int default_hba_coeff_shape[3] = {2, 5, 5};
-
-// The array of coefficients in row-major order ("C"-order).
-//
-const std::complex<double> default_hba_coeff[100] = {
-    std::complex<double>(0.9989322499459223, 0.0003305895124867),
-    std::complex<double>(1.0030546028872600, 0.0002157249025076),
-    std::complex<double>(0.0003002209532403, 0.0007909077657054),
-    std::complex<double>(0.0022051270911392, 0.0003834815341981),
-    std::complex<double>(-0.0003856663268042, 0.0008435910525861),
-    std::complex<double>(0.0004887765294093, 0.0002777796480946),
-    std::complex<double>(-0.0000699366665322, 0.0005136144371953),
-    std::complex<double>(0.0001520602842105, 0.0001303481681886),
-    std::complex<double>(0.0000512381993616, 0.0001550785137302),
-    std::complex<double>(0.0000819244737818, 0.0000466470412396),
-    std::complex<double>(0.0249658150445263, -0.0122024663463393),
-    std::complex<double>(-0.0917825091832822, -0.0062606338208358),
-    std::complex<double>(-0.0083709499453879, -0.0289759752488368),
-    std::complex<double>(-0.0689260153643395, -0.0111348626546314),
-    std::complex<double>(0.0116296166994115, -0.0307342946951178),
-    std::complex<double>(-0.0171249717275797, -0.0080642275561593),
-    std::complex<double>(0.0012408055399100, -0.0191295543986957),
-    std::complex<double>(-0.0051031652662961, -0.0037143632875100),
-    std::complex<double>(-0.0022414352263751, -0.0060474723525871),
-    std::complex<double>(-0.0024377933436567, -0.0012852163337395),
-    std::complex<double>(-0.6730977722052307, 0.0940030437973656),
-    std::complex<double>(0.3711597596859299, 0.0557089394867947),
-    std::complex<double>(0.2119250520015808, 0.2155514942677135),
-    std::complex<double>(0.6727380529527980, 0.0989550572104158),
-    std::complex<double>(-0.0419944347289523, 0.2355624543349744),
-    std::complex<double>(0.1917656461134636, 0.0732470381581913),
-    std::complex<double>(0.0048918921441903, 0.1588912409502319),
-    std::complex<double>(0.0575369727210951, 0.0344677222786687),
-    std::complex<double>(0.0241014578366618, 0.0547046570516960),
-    std::complex<double>(0.0219986510834463, 0.0112189146988984),
-    std::complex<double>(0.0665319393516388, -0.1418009730472832),
-    std::complex<double>(-0.7576728614553603, -0.0472040122949963),
-    std::complex<double>(-0.1017024786435272, -0.3302620837788515),
-    std::complex<double>(-0.5600906156274197, -0.0797555201430585),
-    std::complex<double>(0.0889729243872774, -0.3406964719938829),
-    std::complex<double>(-0.1342560801672904, -0.0515926960946038),
-    std::complex<double>(-0.0149335262655201, -0.2084962323582034),
-    std::complex<double>(-0.0327252678958813, -0.0172371907472848),
-    std::complex<double>(-0.0362395089905272, -0.0661322227928722),
-    std::complex<double>(-0.0141568558526096, -0.0042676979206835),
-    std::complex<double>(0.1121669548152054, 0.0504713119323919),
-    std::complex<double>(0.1882531376700409, 0.0088411256350159),
-    std::complex<double>(0.0066968933526899, 0.1181452711088882),
-    std::complex<double>(0.0981630367567397, 0.0129921405004959),
-    std::complex<double>(-0.0347327225501659, 0.1186585563636635),
-    std::complex<double>(0.0102831315790362, 0.0046275244914932),
-    std::complex<double>(0.0070209144233666, 0.0689639468490938),
-    std::complex<double>(-0.0020239346031291, -0.0025499069613344),
-    std::complex<double>(0.0132702874173192, 0.0207916487187541),
-    std::complex<double>(0.0004387107229914, -0.0017223838914815),
-    std::complex<double>(-0.0004916757488397, 0.0000266213616248),
-    std::complex<double>(0.0006516553273188, -0.0000433166563288),
-    std::complex<double>(-0.0004357897643121, 0.0000320567996700),
-    std::complex<double>(0.0005818285824826, -0.0001021069650381),
-    std::complex<double>(-0.0001047488648808, -0.0000302146563592),
-    std::complex<double>(0.0001593350153828, -0.0000879125663990),
-    std::complex<double>(-0.0000141882506567, -0.0000941521783975),
-    std::complex<double>(-0.0000004226298134, -0.0000245060763932),
-    std::complex<double>(-0.0000177429496833, -0.0000561890408003),
-    std::complex<double>(-0.0000018388829279, 0.0000032387726477),
-    std::complex<double>(0.0162495046881796, -0.0010736997976255),
-    std::complex<double>(-0.0175635905033026, 0.0012997068962173),
-    std::complex<double>(0.0138897851110661, -0.0014876219938565),
-    std::complex<double>(-0.0150211436594772, 0.0029712291209158),
-    std::complex<double>(0.0031705620225488, 0.0004838463688512),
-    std::complex<double>(-0.0034418973689263, 0.0024603729467258),
-    std::complex<double>(0.0003028387544878, 0.0026905629457281),
-    std::complex<double>(0.0006768121359769, 0.0005901486396051),
-    std::complex<double>(0.0004634797107989, 0.0016976603895716),
-    std::complex<double>(0.0003344773954073, -0.0001499932789294),
-    std::complex<double>(-0.1492097398080444, 0.0123735410547393),
-    std::complex<double>(0.1393121453502456, -0.0121117146246749),
-    std::complex<double>(-0.1217628319418324, 0.0222643129255504),
-    std::complex<double>(0.1108579917761457, -0.0262986164183475),
-    std::complex<double>(-0.0273147374272124, 0.0098595182007132),
-    std::complex<double>(0.0208992817013466, -0.0205929453727953),
-    std::complex<double>(-0.0002152227668601, -0.0089220757225133),
-    std::complex<double>(-0.0074792188817697, -0.0043562231368076),
-    std::complex<double>(-0.0012019994038721, -0.0079939660050373),
-    std::complex<double>(-0.0035807498769946, 0.0014801422733613),
-    std::complex<double>(0.1567990061437258, -0.0143275575385193),
-    std::complex<double>(-0.1043118778001582, 0.0106756004832779),
-    std::complex<double>(0.1151024257152241, -0.0225518489392044),
-    std::complex<double>(-0.0593437249231851, 0.0216080058910987),
-    std::complex<double>(0.0142781186223020, -0.0057037138045721),
-    std::complex<double>(0.0151043140114779, 0.0141435752121475),
-    std::complex<double>(-0.0057143555179676, 0.0141142700941743),
-    std::complex<double>(0.0251435557201315, -0.0005753615445942),
-    std::complex<double>(0.0004475745352473, 0.0102135659618127),
-    std::complex<double>(0.0090474375150397, -0.0032177128650026),
-    std::complex<double>(-0.0459124372023251, 0.0044990718645418),
-    std::complex<double>(0.0135433541303599, -0.0021789296923529),
-    std::complex<double>(-0.0306136798186735, 0.0064963361606382),
-    std::complex<double>(-0.0046440676338940, -0.0037281688158807),
-    std::complex<double>(-0.0006372791846825, 0.0008894047150233),
-    std::complex<double>(-0.0181611528840412, -0.0011106177431486),
-    std::complex<double>(0.0032325387394458, -0.0048123509184894),
-    std::complex<double>(-0.0136340313457176, 0.0021185000810664),
-    std::complex<double>(0.0001287985092565, -0.0032079544559908),
-    std::complex<double>(-0.0045503800737417, 0.0015366231416036)};
diff --git a/cpp/lobes/defaultcoefflba.cc b/cpp/lobes/defaultcoefflba.cc
deleted file mode 100644
index d8ea5404ece2ec3bf7da87f2adb21393ca021374..0000000000000000000000000000000000000000
--- a/cpp/lobes/defaultcoefflba.cc
+++ /dev/null
@@ -1,123 +0,0 @@
-// Beam model coefficients converted by convert_coeff.py.
-// Conversion performed on 2011/09/14/08:23:09 UTC using:
-//     convert_coeff.py element_beam_LBA.coeff defaultcoefflba.cc default_lba
-
-#include <complex>
-
-// Center frequency, and frequency range for which the beam model coefficients
-// are valid. The beam model is parameterized in terms of a normalized
-// frequency f' in the range [-1.0, 1.0]. The appropriate conversion is:
-//
-//     f' = (f - center) / range
-//
-const double default_lba_freq_center = 55e6;
-const double default_lba_freq_range = 45e6;
-
-// Shape of the coefficient array: 2x5x5x2 (the size of the last dimension is
-// implied, and always equal to 2).
-//
-const unsigned int default_lba_coeff_shape[3] = {2, 5, 5};
-
-// The array of coefficients in row-major order ("C"-order).
-//
-const std::complex<double> default_lba_coeff[100] = {
-    std::complex<double>(0.9982445079290715, 0.0000650863154389),
-    std::complex<double>(1.0006230902158257, -0.0022053287681416),
-    std::complex<double>(0.0001002692200362, 0.0006838211278268),
-    std::complex<double>(-0.0003660049052840, -0.0008418920419220),
-    std::complex<double>(-0.0010581424498791, 0.0015237878543047),
-    std::complex<double>(0.0007398729642721, 0.0028468649470433),
-    std::complex<double>(-0.0039458389254656, -0.0007048354913730),
-    std::complex<double>(0.0007040177887611, 0.0007856369612188),
-    std::complex<double>(-0.0031701591723043, -0.0010521154166512),
-    std::complex<double>(-0.0007213036752903, -0.0007227764008022),
-    std::complex<double>(0.0550606068782634, 0.0011958385659938),
-    std::complex<double>(-0.0160912944232080, 0.0703645376267940),
-    std::complex<double>(0.0033849565901213, -0.0244636379385135),
-    std::complex<double>(0.0234264238829944, 0.0084068836453700),
-    std::complex<double>(0.0557107413978542, -0.0634701730653090),
-    std::complex<double>(-0.0139549526991330, -0.1175401658864208),
-    std::complex<double>(0.1336911750356096, 0.0202651327657687),
-    std::complex<double>(-0.0113385668361727, -0.0339262369086247),
-    std::complex<double>(0.0962263571740972, 0.0440074333288440),
-    std::complex<double>(0.0313595045238824, 0.0230763038515351),
-    std::complex<double>(-0.8624889445327827, -0.1522883072804402),
-    std::complex<double>(-0.0386800869486029, -0.7569350701887934),
-    std::complex<double>(0.0891332399420108, 0.1876527151756476),
-    std::complex<double>(-0.1012363483900640, -0.1975118891151966),
-    std::complex<double>(-0.6404795825927633, 0.7568775384981410),
-    std::complex<double>(0.0767245154665722, 1.3441875993523555),
-    std::complex<double>(-0.8758406699506004, 0.3350237639226141),
-    std::complex<double>(0.2824832769101577, 0.6821307442669313),
-    std::complex<double>(-0.3144282315609649, -0.2763869580286276),
-    std::complex<double>(-0.1705959031354030, -0.0712085950559831),
-    std::complex<double>(0.4039567648146965, 0.0810473144253429),
-    std::complex<double>(-0.0350803390479135, 0.5214591717801087),
-    std::complex<double>(0.2232030356124932, -0.2248154851829713),
-    std::complex<double>(0.4704343293662089, -0.3552101485419532),
-    std::complex<double>(0.9646419509627557, -0.8095088593139815),
-    std::complex<double>(0.1635280638865702, -1.4854352979459096),
-    std::complex<double>(1.0331569921006993, 0.0509705885336283),
-    std::complex<double>(0.1501121326521990, -0.5193414816770609),
-    std::complex<double>(0.4715775965513117, 0.5077361528286819),
-    std::complex<double>(0.3847391427972284, 0.1136717951238837),
-    std::complex<double>(-0.0756250564248881, 0.0056622911723172),
-    std::complex<double>(-0.1267444401630109, -0.0349676272376008),
-    std::complex<double>(-0.1793752883639813, 0.0720222655359702),
-    std::complex<double>(-0.2678542619793421, 0.3152115802895427),
-    std::complex<double>(-0.3718069213271066, 0.2275266747872172),
-    std::complex<double>(-0.1372223722572021, 0.4314989948093362),
-    std::complex<double>(-0.3316657641578328, -0.1655909947939444),
-    std::complex<double>(-0.2158100484836540, 0.0614504774034524),
-    std::complex<double>(-0.1901597954359592, -0.2294955549701665),
-    std::complex<double>(-0.1864961465389693, -0.0486276177310768),
-    std::complex<double>(-0.0000762326746410, 0.0000118155774181),
-    std::complex<double>(0.0000118903581604, -0.0000251324432498),
-    std::complex<double>(-0.0002204197663391, -0.0000213776348027),
-    std::complex<double>(0.0001477083861977, 0.0000599750510518),
-    std::complex<double>(-0.0003281057522772, -0.0000770207588466),
-    std::complex<double>(0.0003478997686964, 0.0001481982639746),
-    std::complex<double>(-0.0000625695757282, 0.0000249138990722),
-    std::complex<double>(-0.0000960097542525, 0.0002521364065803),
-    std::complex<double>(0.0001275344578325, 0.0000652362392482),
-    std::complex<double>(-0.0003113309221942, 0.0001956734476566),
-    std::complex<double>(0.0029807707669629, -0.0003262084082071),
-    std::complex<double>(0.0001639620574332, -0.0000266272685197),
-    std::complex<double>(0.0076282580587895, 0.0026614359017468),
-    std::complex<double>(-0.0044850263974801, -0.0058337192660638),
-    std::complex<double>(0.0124258438959177, 0.0067985224235178),
-    std::complex<double>(-0.0126349778957970, -0.0100656881493938),
-    std::complex<double>(0.0059031372522229, 0.0008660479915339),
-    std::complex<double>(0.0039660364524413, -0.0100356333791398),
-    std::complex<double>(-0.0020520685193773, -0.0028564379463666),
-    std::complex<double>(0.0121039958869239, -0.0059701468961263),
-    std::complex<double>(-0.0229975846564195, 0.0010565261888195),
-    std::complex<double>(0.0019573207027441, 0.0050550600926414),
-    std::complex<double>(-0.0682274156850413, -0.0758159820140411),
-    std::complex<double>(0.0497303968865466, 0.1019681987654797),
-    std::complex<double>(-0.1757936183439326, -0.1363710820472197),
-    std::complex<double>(0.1765450269056824, 0.1555919358121995),
-    std::complex<double>(-0.1541299429420569, -0.0281422177614844),
-    std::complex<double>(0.0816399676454817, 0.0691599035109852),
-    std::complex<double>(-0.0235110916473515, 0.0306385386726702),
-    std::complex<double>(-0.0474273292450285, 0.0116831908947225),
-    std::complex<double>(0.0333560984394624, -0.0009767086536162),
-    std::complex<double>(0.0141704479374002, -0.0205386534626779),
-    std::complex<double>(0.0562541280098909, 0.0743149092143081),
-    std::complex<double>(-0.0226634801339250, -0.1439026188572270),
-    std::complex<double>(0.1238595124159999, 0.1766108700786397),
-    std::complex<double>(-0.1307647072780430, -0.2090615438301942),
-    std::complex<double>(0.1557916917691289, 0.0646351862895731),
-    std::complex<double>(0.0170294191358757, -0.1027926845803498),
-    std::complex<double>(0.0543537332385954, -0.0366524906364179),
-    std::complex<double>(0.1127180664279469, -0.0176607923511174),
-    std::complex<double>(-0.0126732549319889, 0.0002042370658763),
-    std::complex<double>(-0.0101360135082899, 0.0114084024114141),
-    std::complex<double>(-0.0102147881225462, -0.0176848554302252),
-    std::complex<double>(-0.0051268936720694, 0.0527621533959941),
-    std::complex<double>(-0.0110701836450407, -0.0593085026046026),
-    std::complex<double>(0.0140598301629874, 0.0738668439833535),
-    std::complex<double>(-0.0389912915621699, -0.0301364165752433),
-    std::complex<double>(-0.0462331759359031, 0.0405864871628086),
-    std::complex<double>(-0.0251598701859194, 0.0115712688652445),
-    std::complex<double>(-0.0563476280247398, 0.0079787883434624)};
diff --git a/cpp/lobes/lobeselementresponse.cc b/cpp/lobes/lobeselementresponse.cc
index a8081d0f5aa5e94e4e95d300ceb2af24b48eda2c..a1ec962fb7ebd8151c6648a2fefbe9cb37770c2d 100644
--- a/cpp/lobes/lobeselementresponse.cc
+++ b/cpp/lobes/lobeselementresponse.cc
@@ -1,14 +1,234 @@
+#define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1
+#include <cmath>
+
+#include "config.h"
 #include "lobeselementresponse.h"
 
 #include <map>
+#include <tuple>
+#include <complex>
+#include <H5Cpp.h>
+#include <iostream>
+
+// Anonymous namespace for evaluating the base functions
+namespace {
+double P(int m, int n, double x) {
+  double result = std::assoc_legendre(n, std::abs(m), x);
+
+  if (m < 0) {
+    result *=
+        std::pow(-1, -m) * std::tgamma(n + m + 1) / std::tgamma(n - m + 1);
+  }
+
+  return result;
+}
+
+double Pacc(int m, int n, double x) {
+  return (-(n + m) * (n - m + 1.0) * sqrt(1.0 - x * x) * P(m - 1, n, x) -
+          m * x * P(m, n, x)) /
+         (x * x - 1.0);
+}
+
+std::pair<std::complex<double>, std::complex<double>> F4far_new(int s, int m,
+                                                                int n,
+                                                                double theta,
+                                                                double phi) {
+  double C;
+  if (m) {
+    C = std::sqrt(60.0) * 1.0 / std::sqrt(n * (n + 1.0)) *
+        std::pow(-m / std::abs(m), m);
+  } else {
+    C = std::sqrt(60.0) * 1.0 / std::sqrt(n * (n + 1.0));
+  }
+
+  std::complex<double> q2;
+  std::complex<double> q3;
+
+  // From cpp >= cpp14, complex literals can be used
+  std::complex<double> i_neg = {0.0, -1.0};
+  std::complex<double> i_pos = {0.0, 1.0};
+
+  auto cos_theta = cos(theta);
+  auto sin_theta = sin(theta);
+  auto P_cos_theta = P(std::abs(m), n, cos_theta);
+  auto Pacc_cos_theta = Pacc(std::abs(m), n, cos_theta);
+  auto exp_i_m_phi = exp(i_pos * double(m) * phi);
+
+  if (s == 1) {
+    q2 = C * std::pow(i_neg, -n - 1) * i_pos * double(m) /
+         (sin_theta)*std::sqrt((2. * n + 1) / 2.0 *
+                               std::tgamma(n - std::abs(m) + 1) /
+                               std::tgamma(n + std::abs(m) + 1)) *
+         P_cos_theta * exp_i_m_phi;
+
+    q3 = C * std::pow(i_neg, -n - 1) *
+         std::sqrt((2. * n + 1) / 2.0 * std::tgamma(n - abs(m) + 1) /
+                   std::tgamma(n + abs(m) + 1)) *
+         Pacc_cos_theta * sin_theta * exp_i_m_phi;
+  } else if (s == 2) {
+    q2 = -C * std::pow(i_neg, -n) *
+         std::sqrt((2. * n + 1) / 2.0 * std::tgamma(n - abs(m) + 1) /
+                   std::tgamma(n + abs(m) + 1)) *
+         Pacc_cos_theta * sin_theta * exp_i_m_phi;
+
+    q3 = C * std::pow(i_neg, -n) * i_pos * double(m) / sin_theta *
+         std::sqrt((2. * n + 1) / 2.0 * std::tgamma(n - abs(m) + 1) /
+                   std::tgamma(n + abs(m) + 1)) *
+         P_cos_theta * exp_i_m_phi;
+  }
+
+  return std::make_pair(q2, q3);
+}
+}  // namespace
 
 namespace everybeam {
+LOBESElementResponse::LOBESElementResponse(std::string name) {
+  std::string file_name = std::string("LOBES_") + name + std::string(".h5");
+  std::string path_name = GetPath(file_name.c_str());
+  H5::H5File h5file;
+
+  try {
+    h5file.openFile(path_name.c_str(), H5F_ACC_RDONLY);
+  } catch (const H5::FileIException& e) {
+    throw std::runtime_error("Could not open LOBES coeffcients file: " +
+                             path_name);
+  }
+
+  H5::DataSet dataset = h5file.openDataSet("coefficients");
+  H5::DataSpace dataspace = dataset.getSpace();
+  int nr_elements = dataspace.getSimpleExtentNpoints();
+
+  const std::string REAL("r");
+  const std::string IMAG("i");
+
+  // Create HDF5 complex datatype
+  H5::CompType h5_dcomplex(sizeof(std::complex<double>));
+  h5_dcomplex.insertMember(REAL, 0, H5::PredType::NATIVE_DOUBLE);
+  h5_dcomplex.insertMember(IMAG, sizeof(double), H5::PredType::NATIVE_DOUBLE);
+
+  // Get the number of dimensions in the dataspace.
+  int ndims_coefficients = dataspace.getSimpleExtentNdims();
 
-LOBESElementResponse::LOBESElementResponse(std::string name) {}
+  // Get the dimension size of each dimension in the dataspace and display them.
+  hsize_t dims_coefficients[ndims_coefficients];
+  dataspace.getSimpleExtentDims(dims_coefficients, NULL);
+  coefficients_shape_ = std::vector<unsigned int>(
+      dims_coefficients, dims_coefficients + ndims_coefficients);
+
+  // Read coefficients
+  coefficients_.resize(coefficients_shape_[0], coefficients_shape_[1],
+                       coefficients_shape_[2], coefficients_shape_[3]);
+  dataset.read(coefficients_.data(), h5_dcomplex);
+
+  // Frequencies
+  dataset = h5file.openDataSet("frequencies");
+  dataspace = dataset.getSpace();
+  nr_elements = dataspace.getSimpleExtentNpoints();
+
+  frequencies_.resize(nr_elements);
+  dataset.read(frequencies_.data(), H5::PredType::NATIVE_DOUBLE);
+
+  // nms
+  dataset = h5file.openDataSet("nms");
+  dataspace = dataset.getSpace();
+  nr_elements = dataspace.getSimpleExtentNpoints();
+
+  // Get the number of dimensions in the dataspace.
+  int ndims_nms = dataspace.getSimpleExtentNdims();
+
+  // Get the dimension size of each dimension in the dataspace and display them.
+  hsize_t dims_nms[ndims_nms];
+  dataspace.getSimpleExtentDims(dims_nms, NULL);
+
+  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(
+    double theta, double phi) const {
+  LOBESElementResponse::BaseFunctions base_functions(nms_.size(), 2);
+  base_functions.setZero();
+
+  // Avoid singularities due to theta = 0
+  if (std::abs(theta) < 1e-6) {
+    // TODO: should be handled by a proper warning
+    std::cout << "LOBES is singular for zenith angle = 0rad. Will round "
+              << theta << " to a value of 1e-6" << std::endl;
+    theta = 1e-6;
+  }
+
+  for (size_t i = 0; i < nms_.size(); ++i) {
+    auto nms = nms_[i];
+    std::complex<double> q2, q3;
+    std::tie(q2, q3) = F4far_new(nms.s, nms.m, nms.n, theta, phi);
+    base_functions(i, 0) = q2;
+    base_functions(i, 1) = q3;
+  }
+  return base_functions;
+}
 
 void LOBESElementResponse::Response(
-    double freq, double theta, double phi,
-    std::complex<double> (&response)[2][2]) const {}
+    int element_id, double freq, double theta, double phi,
+    std::complex<double> (&response)[2][2]) const {
+  int freq_idx = FindFrequencyIdx(freq);
+  std::complex<double> xx = {0}, xy = {0}, yx = {0}, yy = {0};
+
+  int nr_rows = basefunctions_.rows();
+  if (nr_rows == 0) {
+    throw std::runtime_error(
+        "Number of rows in basefunctions_ member is 0. Did you run "
+        "SetFieldQuantities?");
+  }
+
+  for (int i = 0; i < nr_rows; ++i) {
+    std::complex<double> q2 = basefunctions_(i, 0);
+    std::complex<double> q3 = basefunctions_(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);
+  }
+
+  response[0][0] = xx;
+  response[0][1] = xy;
+  response[1][0] = yx;
+  response[1][1] = yy;
+}
 
 std::shared_ptr<LOBESElementResponse> LOBESElementResponse::GetInstance(
     std::string name) {
@@ -23,4 +243,11 @@ std::shared_ptr<LOBESElementResponse> LOBESElementResponse::GetInstance(
   return entry->second;
 };
 
+std::string LOBESElementResponse::GetPath(const char* filename) const {
+  std::stringstream ss;
+  ss << EVERYBEAM_DATA_DIR << "/lobes/";
+  ss << filename;
+  return ss.str();
+}
+
 }  // namespace everybeam
diff --git a/cpp/lobes/lobeselementresponse.h b/cpp/lobes/lobeselementresponse.h
index c0d9ee129039f3fb6ddc4db3bd79dc4c1fd6e76c..cf043086e138dadd8eed590ad40a0f8cf1d42c04 100644
--- a/cpp/lobes/lobeselementresponse.h
+++ b/cpp/lobes/lobeselementresponse.h
@@ -1,24 +1,103 @@
-#ifndef LOBES_ELEMENTRESPONSE_H
-#define LOBES_ELEMENTRESPONSE_H
+#ifndef EVERYBEAM_LOBES_ELEMENTRESPONSE_H
+#define EVERYBEAM_LOBES_ELEMENTRESPONSE_H
 
-#include "../elementresponse.h"
+#include <Eigen/Core>
+#include <unsupported/Eigen/CXX11/Tensor>
 
+#include "../fieldresponse.h"
+
+#include <vector>
 #include <memory>
+#include <iostream>
 
 namespace everybeam {
 
 //! Implementation of the Lobes response model
-class LOBESElementResponse : public ElementResponse {
+class LOBESElementResponse : public FieldResponse {
  public:
+  /**
+   * @brief Construct a new LOBESElementResponse object
+   *
+   * @param name (LOFAR) station name, i.e. CS302LBA
+   */
   LOBESElementResponse(std::string name);
 
+  /**
+   * @brief Stub override of the Response method, an element id
+   * is needed to compute the element response.
+   *
+   * @param freq
+   * @param theta
+   * @param phi
+   * @param response
+   */
   virtual void Response(
       double freq, double theta, double phi,
-      std::complex<double> (&response)[2][2]) const final override;
+      std::complex<double> (&response)[2][2]) const final override {
+    throw std::invalid_argument(
+        "LOBESElementResponse::response needs an element_id");
+  };
+
+  /**
+   * @brief Virtual implementation of Response method
+   *
+   * @param element_id ID of element
+   * @param freq Frequency of the plane wave (Hz).
+   * @param theta Angle wrt. z-axis (rad) (NOTE: parameter is just a stub to
+   * match override)
+   * @param phi Angle in the xy-plane wrt. x-axis  (rad) (NOTE: parameter is
+   * just a stub to match override)
+   * @param result Pointer to 2x2 array of Jones matrix
+   */
+  virtual void Response(int element_id, double freq, double theta, double phi,
+                        std::complex<double> (&response)[2][2]) const;
 
   static std::shared_ptr<LOBESElementResponse> GetInstance(std::string name);
-};
 
+  /**
+   * @brief Set field quantities (i.e. the basefunctions) for the LOBES element
+   * response given the direction of interest, specified in (theta, phi) angles.
+   * NOTE: this method overrides the "cached" basefunction_ member every time
+   * when called.
+   *
+   * @param theta Angle wrt. z-axis (rad)
+   * @param phi Angle in the xy-plane wrt. x-axis  (rad)
+   */
+  virtual void SetFieldQuantities(double theta, double phi) final override {
+    basefunctions_ = ComputeBaseFunctions(theta, phi);
+  };
+
+ private:
+  // Typdef of BaseFunctions as Eigen::Array type
+  typedef Eigen::Array<std::complex<double>, Eigen::Dynamic, 2> BaseFunctions;
+  BaseFunctions basefunctions_;
+
+  // Compose the path to a LOBES coefficient file
+  std::string GetPath(const char *) const;
+
+  // Find the closest frequency
+  size_t FindFrequencyIdx(double f) const {
+    auto is_closer = [f](int x, int y) { return abs(x - f) < abs(y - f); };
+    auto result =
+        std::min_element(frequencies_.begin(), frequencies_.end(), is_closer);
+    return std::distance(frequencies_.begin(), result);
+  }
+
+  // Compute the base functions given theta and phi angles
+  BaseFunctions ComputeBaseFunctions(double theta, double phi) const;
+
+  // Store h5 coefficients in coefficients_
+  Eigen::Tensor<std::complex<double>, 4> coefficients_;
+  std::vector<unsigned int> coefficients_shape_;
+
+  // Store h5 frequencies in frequencies_
+  std::vector<double> frequencies_;
+
+  struct nms_t {
+    int n, m, s;
+  };
+  std::vector<nms_t> nms_;
+};
 }  // namespace everybeam
 
 #endif
diff --git a/cpp/lofarreadutils.cc b/cpp/lofarreadutils.cc
index b473ee5bc9b98d1305e009dad3a3311adbbd1daf..ac732dbf76521c7f7b4693e8215789ae49767192 100644
--- a/cpp/lofarreadutils.cc
+++ b/cpp/lofarreadutils.cc
@@ -242,8 +242,15 @@ Antenna::Ptr ReadAntennaField(const Table &table, unsigned int id,
     } else {
       beam_former = std::make_shared<BeamFormerLofarLBA>(coordinate_system);
     }
+  } else if (element_response_model == ElementResponseModel::kLOBES) {
+    // BeamFormer is assigned a FieldResponse, for which common field quantities
+    // can be precomputed
+    beam_former = std::make_shared<BeamFormer>(
+        coordinate_system,
+        std::dynamic_pointer_cast<FieldResponse>(*element_response.get()));
   } else {
-    // Tiles / element should be kept unique, so work with generic BeamFormer
+    // All Tiles / elements should be kept unique, so work with generic
+    // BeamFormer
     beam_former = std::make_shared<BeamFormer>(coordinate_system);
   }
 
diff --git a/cpp/station.cc b/cpp/station.cc
index a2cc6a81ca0331bd042931b61729594baf00fdbc..07602dd15dfaa59f2b86a3f59af61395e97de8b2 100644
--- a/cpp/station.cc
+++ b/cpp/station.cc
@@ -23,8 +23,6 @@
 #include "station.h"
 #include "common/mathutils.h"
 #include "beamformerlofar.h"
-// #include "beamformerlofarhba.h"
-// #include "beamformerlofarlba.h"
 
 #include "hamaker/hamakerelementresponse.h"
 #include "oskar/oskarelementresponse.h"
@@ -40,11 +38,19 @@ Station::Station(const std::string &name, const vector3r_t &position,
       phase_reference_(position),
       element_response_model_(model),
       element_response_(nullptr) {
-  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}};
   ncp_pol0_.reset(new coords::ITRFDirection(ncppol0));
+
+  // LOBES is currently only supported for CS302LBA. Check this.
+  if (model == ElementResponseModel::kLOBES && name != "CS302LBA") {
+    std::cout << "Switching to Hamaker, LOBES is not yet "
+                 "supported for station "
+              << name << std::endl;
+    element_response_model_ = ElementResponseModel::kHamaker;
+  }
+  SetResponseModel(element_response_model_);
 }
 
 void Station::SetResponseModel(const ElementResponseModel model) {
diff --git a/cpp/test/tlofar_lba.cc b/cpp/test/tlofar_lba.cc
index a9c89b927ba2e41f6adf9c5a2911db92424ec5f6..71400d841e075f094533bac1de6a01c41d496b88 100644
--- a/cpp/test/tlofar_lba.cc
+++ b/cpp/test/tlofar_lba.cc
@@ -6,6 +6,7 @@
 #include "../elementresponse.h"
 #include "../station.h"
 #include "../common/types.h"
+#include "../../external/npy.hpp"
 
 #include "config.h"
 #include <complex>
@@ -25,8 +26,29 @@ using everybeam::telescope::Telescope;
 
 BOOST_AUTO_TEST_SUITE(tlofar_lba)
 
-BOOST_AUTO_TEST_CASE(load_lofar) {
+// Properties extracted from MS
+double time = 4.92183348e+09;
+double frequency = 57812500.;
+double ra(-1.44194878), dec(0.85078091);
+
+// Properties of grid
+std::size_t width(4), height(4);
+double dl(0.5 * M_PI / 180.), dm(0.5 * M_PI / 180.), shift_l(0.), shift_m(0.);
+
+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};
+
+BOOST_AUTO_TEST_CASE(test_hamaker) {
   Options options;
+  // Only checks if we are defaulting to LOBES. Effectively,
+  // all the computations will be done as if the Hamaker model was chosen
+  // except for station 20 (CS302LBA)
   options.element_response_model = ElementResponseModel::kHamaker;
 
   casacore::MeasurementSet ms(LOFAR_LBA_MOCK_MS);
@@ -42,22 +64,6 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
   const LOFAR& lofartelescope = static_cast<const LOFAR&>(*telescope.get());
   BOOST_CHECK_EQUAL(lofartelescope.GetStation(0)->GetName(), "CS001LBA");
 
-  // Properties extracted from MS
-  double time = 4.92183348e+09;
-  double frequency = 57812500.;
-  std::size_t width(4), height(4);
-  double ra(-1.44194878), dec(0.85078091), dl(0.5 * M_PI / 180.),
-      dm(0.5 * M_PI / 180.), shift_l(0.), shift_m(0.);
-
-  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};
-
   // Reference solution obtained with commit sha
   // 70a286e7dace4616417b0e973a624477f15c9ce3
   //
@@ -109,4 +115,49 @@ BOOST_AUTO_TEST_CASE(load_lofar) {
   }
 }
 
+BOOST_AUTO_TEST_CASE(test_lobes) {
+  Options options;
+  // Effectively, all the computations will be done as if the Hamaker model was
+  // chosen except for station 20 (CS302LBA)
+  options.element_response_model = ElementResponseModel::kLOBES;
+
+  casacore::MeasurementSet ms(LOFAR_LBA_MOCK_MS);
+
+  // Load LOFAR Telescope
+  std::unique_ptr<Telescope> telescope = Load(ms, options);
+
+  // Extract Station 20, should be station CS302LBA
+  const LOFAR& lofartelescope = static_cast<const LOFAR&>(*telescope.get());
+  BOOST_CHECK_EQUAL(lofartelescope.GetStation(20)->GetName(), "CS302LBA");
+
+  // Gridded response
+  std::unique_ptr<GriddedResponse> grid_response =
+      telescope->GetGriddedResponse(coord_system);
+
+  // Define buffer and get gridded responses
+  std::vector<std::complex<float>> antenna_buffer_single(
+      grid_response->GetStationBufferSize(1));
+
+  // Get the gridded response for station 20 (of course!)
+  grid_response->CalculateStation(antenna_buffer_single.data(), time, frequency,
+                                  20, 0);
+
+  // 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.6094082, 0.2714097},
+      {-0.9981958, 1.081614},
+      {-0.5575241, -0.3563573},
+      {-0.6945726, 0.1506443}};
+  std::size_t offset_13 = (3 + 1 * width) * 4;
+
+  for (std::size_t i = 0; i < 4; ++i) {
+    BOOST_CHECK(std::abs(antenna_buffer_single[offset_13 + i] -
+                         everybeam_ref_p13[i]) < 1e-6);
+  }
+  // const long unsigned leshape[] = {(long unsigned int)width, height, 2, 2};
+  // npy::SaveArrayAsNumpy("lobes_station_response.npy", false, 4, leshape,
+  //                       antenna_buffer_single);
+}
+
 BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file
diff --git a/python/lobes/CMakeLists.txt b/python/lobes/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..499176718cc52701401030bcb16e9ca62d3681f9
--- /dev/null
+++ b/python/lobes/CMakeLists.txt
@@ -0,0 +1,10 @@
+project(lobes_pybind11)
+
+# Create the binding library
+pybind11_add_module(pylobes 
+                    lobes.cc
+                    )
+
+
+target_include_directories(pylobes PUBLIC ${HDF5_LIBRARIES})
+target_link_libraries(pylobes PUBLIC ${HDF5_LIBRARIES})
diff --git a/python/lobes/lobes.cc b/python/lobes/lobes.cc
index 0aa5138fa9b13f385d2df9e95828fa1b819fc0c6..4dc8fd6606dec23daaaa8a5713c8f4ccf5e50e2a 100644
--- a/python/lobes/lobes.cc
+++ b/python/lobes/lobes.cc
@@ -1,19 +1,18 @@
 #define __STDCPP_WANT_MATH_SPEC_FUNCS__ 1
 #include <cmath>
 
+#include <Eigen/Core>
+#include <pybind11/pybind11.h>
+#include <pybind11/eigen.h>
+
 #include <iostream>
 #include <complex>
-
-#include "ElementResponse.h"
-
 #include <H5Cpp.h>
+#include <tuple>
 
-#include "lobes.h"
-
+namespace py = pybind11;
 using namespace std::complex_literals;
 
-#include <tuple>
-
 // python style enumerate function
 // To make it possible to write:
 //   for (auto [i, v] : enumerate(iterable)) {...}
@@ -129,202 +128,7 @@ std::pair<Eigen::VectorXcd, Eigen::VectorXcd> F4far_new(
   return std::make_pair(q2, q3);
 }
 
-Eigen::Array<std::complex<double>, 2, 2> ElementResponseLBA(double freq,
-                                                            double theta,
-                                                            double phi) {
-  Eigen::Array<std::complex<double>, 2, 2> response;
-
-  typedef std::complex<double> Response[2][2];
-  Response &response_ = *((Response *)response.data());
-
-  everybeam::ElementResponseLBA(freq, theta, phi, response_);
-  return std::move(response);
-}
-
-LobesBeamModel::LobesBeamModel(const std::string &data_file_name) {
-  H5::H5File h5file(data_file_name.c_str(), H5F_ACC_RDONLY);
-
-  H5::DataSet dataset = h5file.openDataSet("coefficients");
-  H5::DataSpace dataspace = dataset.getSpace();
-  int nr_elements = dataspace.getSimpleExtentNpoints();
-
-  const std::string REAL("r");
-  const std::string IMAG("i");
-
-  H5::CompType h5_dcomplex(sizeof(std::complex<double>));
-  h5_dcomplex.insertMember(REAL, 0, H5::PredType::NATIVE_DOUBLE);
-  h5_dcomplex.insertMember(IMAG, sizeof(double), H5::PredType::NATIVE_DOUBLE);
-
-  /*
-   * Get the number of dimensions in the dataspace.
-   */
-  int ndims_coefficients = dataspace.getSimpleExtentNdims();
-  /*
-   * Get the dimension size of each dimension in the dataspace and
-   * display them.
-   */
-  hsize_t dims_coefficients[ndims_coefficients];
-  dataspace.getSimpleExtentDims(dims_coefficients, NULL);
-
-  std::cout << "Dimensions of coefficients:" << std::endl << "    ";
-  for (int i = 0; i < ndims_coefficients; i++) {
-    std::cout << (size_t)(dims_coefficients[i]);
-    if (i < ndims_coefficients - 1)
-      std::cout << " x ";
-    else
-      std::cout << std::endl;
-  }
-
-  // Reshape coefficients
-  size_t nr_elements_first_dims = 1;
-  for (int i = 0; i < ndims_coefficients - 1; i++)
-    nr_elements_first_dims *= dims_coefficients[i];
-  std::cout << nr_elements_first_dims << " x "
-            << (size_t)(dims_coefficients[ndims_coefficients - 1]) << std::endl;
-
-  m_coefficients =
-      Eigen::ArrayXXcd(nr_elements_first_dims,
-                       (size_t)dims_coefficients[ndims_coefficients - 1]);
-  m_coefficients_shape = std::vector<size_t>(
-      dims_coefficients, dims_coefficients + ndims_coefficients);
-  dataset.read(m_coefficients.data(), h5_dcomplex);
-
-  //=====================================
-  dataset = h5file.openDataSet("frequencies");
-  dataspace = dataset.getSpace();
-  nr_elements = dataspace.getSimpleExtentNpoints();
-
-  m_frequencies = std::vector<double>(nr_elements);
-
-  dataset.read(m_frequencies.data(), H5::PredType::NATIVE_DOUBLE);
-  std::cout << "Frequencies:" << std::endl;
-  for (auto freq : m_frequencies) {
-    std::cout << freq << " ";
-  }
-  std::cout << std::endl << std::endl;
-
-  dataset = h5file.openDataSet("nms");
-  dataspace = dataset.getSpace();
-  nr_elements = dataspace.getSimpleExtentNpoints();
-
-  /*
-   * Get the number of dimensions in the dataspace.
-   */
-  int ndims_nms = dataspace.getSimpleExtentNdims();
-  /*
-   * Get the dimension size of each dimension in the dataspace and
-   * display them.
-   */
-  hsize_t dims_nms[ndims_nms];
-  dataspace.getSimpleExtentDims(dims_nms, NULL);
-
-  std::cout << "Dimensions of nms:" << ndims_nms << std::endl << "    ";
-
-  for (int i = 0; i < ndims_nms; i++) {
-    std::cout << (size_t)(dims_nms[i]);
-    if (i < ndims_nms - 1) std::cout << " x ";
-  }
-  std::cout << std::endl;
-
-  m_nms = py::array_t<int>({dims_nms[0], dims_nms[1]});
-  dataset.read(m_nms.mutable_data(), H5::PredType::NATIVE_INT);
-}
-
-py::array_t<std::complex<double>> LobesBeamModel::eval(
-    py::EigenDRef<const Eigen::ArrayXd> theta,
-    py::EigenDRef<const Eigen::ArrayXd> phi) {
-  double beta;
-  for (auto [i, freq] : enumerate(m_frequencies)) {
-    beta = 2.0 * M_PI * freq / 2.99792458e8;
-  }
-
-  beta = 2.0 * M_PI * m_frequencies[0] / 2.99792458e8;
-
-  int s, m, n;
-
-  auto nms = m_nms.unchecked<2>();
-
-  //     auto result = py::array_t<std::complex<double>>({2*theta.size(),
-  //     m_nms.shape(0)}); Eigen::Map<Eigen::MatrixXcd>
-  //     Base(result.mutable_data(), 2*theta.size(), m_nms.shape(0));
-
-  Eigen::MatrixXcd Base(2 * theta.size(), m_nms.shape(0));
-
-  for (int i = 0; i < m_nms.shape(0); i++) {
-    int n = nms(i, 0);
-    int m = nms(i, 1);
-    int s = nms(i, 2);
-    //         std::cout << i << ": (" << n << ", " << m << ", " << s << ")" <<
-    //         std::endl;
-    Eigen::VectorXcd e_phi, e_theta;
-    // Compute SWF
-    std::tie(e_theta, e_phi) = F4far_new(s, m, n, theta, phi, beta);
-    // Interleave
-    Base.col(i)(Eigen::seq(0, Eigen::last, 2)) = e_theta;
-    Base.col(i)(Eigen::seq(1, Eigen::last, 2)) = e_phi;
-  }
-
-  std::cout << Base.rows() << ", " << m_coefficients.rows() << std::endl;
-  std::cout << Base.cols() << ", " << m_coefficients.cols() << std::endl;
-
-  // Create array to return
-  auto result =
-      py::array_t<std::complex<double>>({Base.rows(), m_coefficients.rows()});
-
-  // Map the result array to an Eigen array, to be able to assign it the result
-  // of an Eigen matrix multiplication
-  Eigen::Map<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic,
-                           Eigen::RowMajor>>
-      result_matrix(result.mutable_data(), Base.rows(), m_coefficients.rows());
-
-  // Multiply base functions by coefficients
-  result_matrix = Base * m_coefficients.matrix().transpose();
-
-  return result;
-}
-
-// py::array_t<std::complex<double>> LobesBeamModel::eval(py::array_t<double>
-// theta, py::array_t<double> phi)
-// {
-//     py::buffer_info theta_buffer = theta.request();
-//     py::buffer_info phi_buffer = phi.request();
-//
-//     if (theta_buffer.ndim != 1 || phi_buffer.ndim != 1)
-//         throw std::runtime_error("Number of dimensions must be one");
-//
-//     if (theta_buffer.size != phi_buffer.size)
-//         throw std::runtime_error("Input shapes must match");
-//
-//     /* No pointer is passed, so NumPy will allocate the buffer */
-//     int size = theta_buffer.size;
-//     auto result = py::array_t<std::complex<double>>({size, 2, 2});
-//
-//     double beta = 1.0;
-//
-//     Eigen::Map<Eigen::ArrayXd> theta_((double *) theta_buffer.ptr, size);
-//     Eigen::Map<Eigen::ArrayXd> phi_(size, (double *) phi_buffer.ptr);
-//
-// //     std::pair<Eigen::VectorXcd,Eigen::VectorXcd> F4far_new(s, m, n, theta,
-// phi, beta)
-//
-//     int s,m,n;
-//
-// //     F4far_new(s, m, n,  theta_, phi_, beta);
-//
-// //     py::buffer_info buf3 = result.request();
-// //
-// //     double *ptr1 = (double *) buf1.ptr,
-// //            *ptr2 = (double *) buf2.ptr,
-// //            *ptr3 = (double *) buf3.ptr;
-// //
-// //     for (size_t idx = 0; idx < buf1.shape[0]; idx++)
-// //         ptr3[idx] = ptr1[idx] + ptr2[idx];
-//
-//     return result;
-//
-// }
-
-PYBIND11_MODULE(lobes, m) {
+PYBIND11_MODULE(pylobes, m) {
   m.doc() = "LOBES module";  // optional module docstring
 
   m.def("plustwo", &plustwo, "A function which adds two to a number");
@@ -339,9 +143,4 @@ PYBIND11_MODULE(lobes, m) {
   m.def("P", &P, "Legendre");
   m.def("Pacc", &Pacc, "Legendre");
   m.def("F4far_new", &F4far_new, "F4far_new");
-  m.def("ElementResponseLBA", &ElementResponseLBA, "ElementResponseLBA");
-
-  py::class_<LobesBeamModel>(m, "LobesBeamModel")
-      .def(py::init<const std::string &>())
-      .def("eval", &LobesBeamModel::eval);
 }
diff --git a/python/lobes/lobes.h b/python/lobes/lobes.h
deleted file mode 100644
index d16241526dd7db27ce3f2282893ceba5848c9109..0000000000000000000000000000000000000000
--- a/python/lobes/lobes.h
+++ /dev/null
@@ -1,42 +0,0 @@
-#include <Eigen/Core>
-#include <pybind11/pybind11.h>
-#include <pybind11/eigen.h>
-
-namespace py = pybind11;
-
-//! (Virtual) Beam model class
-class BeamModel {
- public:
-  BeamModel() {}
-
-  virtual py::array_t<std::complex<double>> eval(
-      py::EigenDRef<const Eigen::ArrayXd> theta,
-      py::EigenDRef<const Eigen::ArrayXd> phi) = 0;
-};
-
-//! Lobes beam model, wrapped with pybind11
-class LobesBeamModel : public BeamModel {
- public:
-  LobesBeamModel(const std::string &data_file_name);
-
-  py::array_t<std::complex<double>> eval(
-      py::EigenDRef<const Eigen::ArrayXd> theta,
-      py::EigenDRef<const Eigen::ArrayXd> phi) override;
-
- private:
-  Eigen::ArrayXXcd m_coefficients;
-  std::vector<size_t> m_coefficients_shape;
-
-  std::vector<double> m_frequencies;
-  py::array_t<int> m_nms;
-};
-
-//! Lobes beam model, not implemented!
-class HamakerBeamModel : public BeamModel {
- public:
-  HamakerBeamModel() {}
-
-  py::array_t<std::complex<double>> eval(
-      py::EigenDRef<const Eigen::ArrayXd> theta,
-      py::EigenDRef<const Eigen::ArrayXd> phi) override {}
-};
diff --git a/scripts/download_lobes_coeffs.sh b/scripts/download_lobes_coeffs.sh
new file mode 100755
index 0000000000000000000000000000000000000000..0f10eabd606e7c561203f1d9e6695045dddf321d
--- /dev/null
+++ b/scripts/download_lobes_coeffs.sh
@@ -0,0 +1,19 @@
+#!/bin/sh
+#
+# Author: Jakob Maljaars
+# Email: jakob.maljaars_@_stcorp.nl
+
+# Script for downloading and extracting the LOBES coefficients
+
+set -e
+
+SCRIPT_PATH=$(dirname "$0")
+cd $SCRIPT_PATH
+
+# Make directory coeffs/lobes in source directory
+cd ..
+mkdir -p coeffs/lobes
+cd coeffs/lobes
+
+# Download the h5 coefficient files in case they're not present (-nc option) in the coeffs/lobes directory
+wget -q -r -nc -nH --no-parent --cut-dirs=3 --accept-regex=LOBES_.*\.h5 --reject index.html -e robots=off https://www.astron.nl/citt/EveryBeam/lobes/
\ No newline at end of file