diff --git a/cpp/antenna.cc b/cpp/antenna.cc
index 6f100299751895270dd2ef148bd71fabb189d5ac..0d99cd30366a9250876e2cedc0736f594e48dea0 100644
--- a/cpp/antenna.cc
+++ b/cpp/antenna.cc
@@ -3,6 +3,7 @@
 #include "common/mathutils.h"
 
 namespace everybeam {
+
 vector3r_t Antenna::TransformToLocalDirection(const vector3r_t &direction) {
   vector3r_t local_direction{
       dot(coordinate_system_.axes.p, direction),
@@ -12,4 +13,34 @@ vector3r_t Antenna::TransformToLocalDirection(const vector3r_t &direction) {
 
   return local_direction;
 }
-}  // namespace everybeam
\ No newline at end of file
+
+void Antenna::Transform(const CoordinateSystem &coordinate_system) {
+  coordinate_system_.axes.p =
+      coordinate_system_.axes.p[0] * coordinate_system.axes.p +
+      coordinate_system_.axes.p[1] * coordinate_system.axes.q +
+      coordinate_system_.axes.p[2] * coordinate_system.axes.r;
+
+  coordinate_system_.axes.q =
+      coordinate_system_.axes.q[0] * coordinate_system.axes.p +
+      coordinate_system_.axes.q[1] * coordinate_system.axes.q +
+      coordinate_system_.axes.q[2] * coordinate_system.axes.r;
+
+  coordinate_system_.axes.r =
+      coordinate_system_.axes.r[0] * coordinate_system.axes.p +
+      coordinate_system_.axes.r[1] * coordinate_system.axes.q +
+      coordinate_system_.axes.r[2] * coordinate_system.axes.r;
+
+  coordinate_system_.origin =
+      coordinate_system.origin +
+      coordinate_system_.origin[0] * coordinate_system.axes.p +
+      coordinate_system_.origin[1] * coordinate_system.axes.q +
+      coordinate_system_.origin[2] * coordinate_system.axes.r;
+
+  phase_reference_position_ =
+      coordinate_system.origin +
+      phase_reference_position_[0] * coordinate_system.axes.p +
+      phase_reference_position_[1] * coordinate_system.axes.q +
+      phase_reference_position_[2] * coordinate_system.axes.r;
+}
+
+}  // namespace everybeam
diff --git a/cpp/antenna.h b/cpp/antenna.h
index ebc6658e933b84d360bdeaae3fb8351d17f98f05..e4a76d62eb4ac731a2c4c0176073ebb8eba1409b 100644
--- a/cpp/antenna.h
+++ b/cpp/antenna.h
@@ -115,6 +115,37 @@ class Antenna {
         phase_reference_position_(phase_reference_position),
         enabled_{true, true} {}
 
+  /**
+   * @brief Makes a copy of this Antenna object
+   *
+   * The method is virtual, so that copies can be created from a pointer
+   * to the base (Antenna) class.
+   * The original remains unchanged, therefore the method is const.
+   * The method has no implementation in the Antenna class, because
+   * Antenna is abstract, so no copy can be instantiated.
+   *
+   * This method is used by the ExtractAntenna method of the BeamFormer
+   * class to create a copy of one of the Antennas it contains.
+   */
+  virtual Ptr Clone() const = 0;
+
+  /**
+   * @brief Transform internal coordinate systems and positions
+   *
+   * @param coordinate_system to apply in the transformation
+   *
+   * This method is used by BeamFormer::ExtractAntenna to lift
+   * an antenna out of the beamformer.
+   *
+   * The transformation is needed because the coordinate system of
+   * an antenna in a beamformer is expressed in terms of
+   * the coordinate system of the beamformer.
+   * To turn an embedded antenna into a stand-alone antenna,
+   * the coordinate system of the beamformer needs to be
+   * applied to the coordinate system of the antenna
+   */
+  void Transform(const CoordinateSystem &coordinate_system);
+
   /**
    * @brief Compute the %Antenna Response
    *
diff --git a/cpp/beamformer.cc b/cpp/beamformer.cc
index daf3048c00ccee7b6f88067b620a5da9cdc6c653..8d887ffe1284a738b9abe4593ea130bab954afb7 100644
--- a/cpp/beamformer.cc
+++ b/cpp/beamformer.cc
@@ -6,6 +6,25 @@
 #include <cmath>
 
 namespace everybeam {
+
+Antenna::Ptr BeamFormer::Clone() const {
+  auto beamformer_clone = BeamFormer::Ptr(
+      new BeamFormer(coordinate_system_, phase_reference_position_));
+
+  // antennas_ is a vector of pointers to Antennas, so
+  // this creates a shallow copy, in the sense that
+  // the antennas are not copied, only the pointers.
+  beamformer_clone->antennas_ = antennas_;
+
+  return beamformer_clone;
+}
+
+Antenna::Ptr BeamFormer::ExtractAntenna(size_t antenna_index) const {
+  Antenna::Ptr antenna = antennas_[antenna_index]->Clone();
+  antenna->Transform(coordinate_system_);
+  return antenna;
+}
+
 vector3r_t BeamFormer::TransformToLocalPosition(const vector3r_t &position) {
   // Get antenna position relative to coordinate system origin
   vector3r_t dposition{position[0] - coordinate_system_.origin[0],
diff --git a/cpp/beamformer.h b/cpp/beamformer.h
index 6bd52d91bce5e077aab5e90109b0fffc0ba85d56..a9811e25468531edfa773097bafa98e4c58aca9c 100644
--- a/cpp/beamformer.h
+++ b/cpp/beamformer.h
@@ -52,6 +52,8 @@ class BeamFormer : public Antenna {
         TransformToLocalPosition(phase_reference_position_);
   }
 
+  Antenna::Ptr Clone() const override;
+
   /**
    * @brief Add an antenna to the m_antenna array.
    *
@@ -59,6 +61,24 @@ class BeamFormer : public Antenna {
    */
   void AddAntenna(Antenna::Ptr antenna) { antennas_.push_back(antenna); }
 
+  /**
+   * @brief Extracts an antenna from the beamformer
+   *
+   * @param antenna_index index of antenna to extact
+   * @returns pointer to a copy of antenna with index antenna_index
+   *
+   * The antenna is extracted such that it can be used stand-alone,
+   * independent of the beamformer. The coordinate system of the extracted
+   * antenna is transformed from internal representation to external
+   * representation by application of the beamformer coordinate system to
+   * the antenna coordinate system.
+   *
+   * The returned antenna can be either an Element or a BeamFormer.
+   *
+   * The beamformer itself remains unchanged.
+   */
+  Antenna::Ptr ExtractAntenna(size_t antenna_index) const;
+
  protected:
   vector3r_t
       local_phase_reference_position_;  // in coordinate system of Antenna
diff --git a/cpp/beamformeridenticalantennas.cc b/cpp/beamformeridenticalantennas.cc
index 3aa01a7a8722639d0cfa56ea2c13cfafedde2400..75b3ba99552f9979d320a6746cb06646594ae7af 100644
--- a/cpp/beamformeridenticalantennas.cc
+++ b/cpp/beamformeridenticalantennas.cc
@@ -7,6 +7,13 @@
 
 namespace everybeam {
 
+Antenna::Ptr BeamFormerIdenticalAntennas::Clone() const {
+  auto beamformer_clone = std::make_shared<BeamFormerIdenticalAntennas>(
+      coordinate_system_, phase_reference_position_);
+  beamformer_clone->antennas_ = antennas_;
+  return beamformer_clone;
+}
+
 matrix22c_t BeamFormerIdenticalAntennas::LocalResponse(
     real_t time, real_t freq, const vector3r_t &direction,
     const Options &options) const {
diff --git a/cpp/beamformeridenticalantennas.h b/cpp/beamformeridenticalantennas.h
index 31185fd9abb3cb46bd3e46bed97ac63f0ccb9615..26aedab2d48ff2c16423d3c72965b5a2bf0a4ce7 100644
--- a/cpp/beamformeridenticalantennas.h
+++ b/cpp/beamformeridenticalantennas.h
@@ -35,6 +35,8 @@ class BeamFormerIdenticalAntennas : public BeamFormer {
   BeamFormerIdenticalAntennas(vector3r_t phase_reference_position)
       : BeamFormer(phase_reference_position) {}
 
+  Antenna::Ptr Clone() const override;
+
  private:
   // Compute the BeamFormer response in certain direction of arrival (ITRF, m)
   // and return (Jones) matrix of response
diff --git a/cpp/element.cc b/cpp/element.cc
index f3c119aea66fa4a3ecc300411912feefef55941e..5803d1ce03fad5d33da7e32d0bd27cc8ffea7e7b 100644
--- a/cpp/element.cc
+++ b/cpp/element.cc
@@ -2,6 +2,15 @@
 #include "common/mathutils.h"
 
 namespace everybeam {
+
+Antenna::Ptr Element::Clone() const {
+  auto element_clone =
+      Element::Ptr(new Element(coordinate_system_, element_response_, id_));
+  element_clone->enabled_[0] = enabled_[0];
+  element_clone->enabled_[1] = enabled_[1];
+  return element_clone;
+}
+
 matrix22c_t Element::LocalResponse(real_t time, real_t freq,
                                    const vector3r_t &direction, size_t id,
                                    const Options &options) const {
diff --git a/cpp/element.h b/cpp/element.h
index 39c0ec32df013f7437325f70b6c6ac366eddce2f..6e8e07223d14424ab4ffba7f8938d57c5aee701a 100644
--- a/cpp/element.h
+++ b/cpp/element.h
@@ -32,6 +32,8 @@ class Element : public Antenna {
         id_(id),
         element_response_(element_response) {}
 
+  Antenna::Ptr Clone() const override;
+
   /**
    * @brief Compute the local response of the element.
    *
diff --git a/cpp/lofarreadutils.cc b/cpp/lofarreadutils.cc
index 2fa2b397d181d28af8657b47becec9b1dd52abd7..9344c16bad0a9a8e7c8214d2f14cee329aafeb8f 100644
--- a/cpp/lofarreadutils.cc
+++ b/cpp/lofarreadutils.cc
@@ -313,15 +313,6 @@ Station::Ptr ReadLofarStation(const MeasurementSet &ms, unsigned int id,
 
     station->SetAntenna(beam_former);
 
-    size_t field_id = 0;
-    size_t element_id = 0;
-    Antenna::CoordinateSystem coordinate_system =
-        common::ReadCoordinateSystem(tab_field, field_id);
-    auto model = station->GetElementResponse();
-    // TODO: rotate coordinate system for antenna
-    auto element =
-        Element::Ptr(new Element(coordinate_system, model, element_id));
-    station->SetElement(element);
   } else if (telescope_name == "AARTFAAC") {
     ROScalarColumn<String> ant_type_col(common::GetSubTable(ms, "OBSERVATION"),
                                         "AARTFAAC_ANTENNA_TYPE");
diff --git a/cpp/msv3readutils.cc b/cpp/msv3readutils.cc
index ddcfe8519bdeae8677b11a02a93a5e39cd536265..4e34394993a2d44ca7ee13d064d19ef9be2566fa 100644
--- a/cpp/msv3readutils.cc
+++ b/cpp/msv3readutils.cc
@@ -137,15 +137,6 @@ Station::Ptr ReadMSv3Station(const MeasurementSet &ms, unsigned int id,
 
   station->SetAntenna(beam_former);
 
-  size_t field_id = 0;
-  size_t element_id = 0;
-  Antenna::CoordinateSystem coordinate_system =
-      common::ReadCoordinateSystem(tab_phased_array, field_id);
-  auto element_response = station->GetElementResponse();
-  auto element = Element::Ptr(
-      new Element(coordinate_system, element_response, element_id));
-  station->SetElement(element);
-
   return station;
 }
 
diff --git a/cpp/station.cc b/cpp/station.cc
index 38cd35a62d9ae43092170329f82508357dfbd93d..332d97341b1f029151aaf3aa6b5ac0fe16890a73 100644
--- a/cpp/station.cc
+++ b/cpp/station.cc
@@ -80,6 +80,21 @@ const vector3r_t &Station::GetPhaseReference() const {
   return phase_reference_;
 }
 
+void Station::SetAntenna(Antenna::Ptr antenna) {
+  antenna_ = antenna;
+
+  // The antenna can be either an Element or a BeamFormer
+  // If it is a BeamFormer we recursively extract the first antenna
+  // until we have an Element.
+  // The extraction returns copies so antenna_ remains unchanged.
+  // The element that is found is used in ComputeElementResponse to
+  // compute the element response.
+  while (auto beamformer = std::dynamic_pointer_cast<BeamFormer>(antenna)) {
+    antenna = beamformer->ExtractAntenna(0);
+  }
+  element_ = std::dynamic_pointer_cast<Element>(antenna);
+}
+
 // ========================================================
 matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
                                             const vector3r_t &direction,
@@ -102,13 +117,20 @@ matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
 matrix22c_t Station::ComputeElementResponse(real_t time, real_t freq,
                                             const vector3r_t &direction,
                                             const bool rotate) const {
-  //     if (rotate)
-  //       return itsElement->response(time, freq, direction)
-  //           * Rotation(time, direction);
-  //     else
-  //       return itsElement->response(time, freq, direction);
+  Antenna::Options options;
+  options.rotate = rotate;
 
-  return element_->Response(time, freq, direction);
+  if (options.rotate) {
+    vector3r_t ncp_ = NCP(time);
+    vector3r_t east = normalize(cross(ncp_, direction));
+    vector3r_t north = cross(direction, east);
+    options.east = east;
+    options.north = north;
+  }
+
+  matrix22c_t response;
+  response = element_->Response(time, freq, direction, options);
+  return response;
 }
 
 matrix22c_t Station::Response(real_t time, real_t freq,
@@ -129,14 +151,6 @@ matrix22c_t Station::Response(real_t time, real_t freq,
 
   matrix22c_t response = antenna_->Response(time, freq, direction, options);
 
-  //     if (rotate) {
-  //         std::cout << "rotate" << std::endl;
-  //         auto r = Rotation(time, direction);
-  //         std::cout << r[0][0] << ", " << r[0][1] << std::endl;
-  //         std::cout << r[1][0] << ", " << r[1][1] << std::endl;
-  //         response = response * r;
-  //         std::cout << response[0][0] << std::endl;
-  //     }
   return response;
 }
 
diff --git a/cpp/station.h b/cpp/station.h
index 85197a712e06364efe6d1f4ddd3ce0c3983223ad..105b58421727bf3c983b27f831e0f534c7f50063 100644
--- a/cpp/station.h
+++ b/cpp/station.h
@@ -333,12 +333,9 @@ class Station {
   }
 
   //! Set antenna attribute, usually a BeamFormer, but can also be an Element
-  void SetAntenna(Antenna::Ptr antenna) { antenna_ = antenna; }
+  void SetAntenna(Antenna::Ptr antenna);
 
-  Antenna::Ptr GetAntenna() { return antenna_; }
-
-  //! Set Element attribute
-  void SetElement(Element::Ptr element) { element_ = element; }
+  Antenna::Ptr GetAntenna() const { return antenna_; }
 
  private:
   vector3r_t NCP(real_t time) const;