diff --git a/CEP/Calibration/pystationresponse/src/__init__.py b/CEP/Calibration/pystationresponse/src/__init__.py
index de29c1e1bd5402e5f631f07deaf20e6aa4ae0224..a3e5a59bbd6a680a804bc8102892ca96e6f8495b 100755
--- a/CEP/Calibration/pystationresponse/src/__init__.py
+++ b/CEP/Calibration/pystationresponse/src/__init__.py
@@ -82,6 +82,15 @@ class stationresponse(object):
         """
         self._response.setRefDelay(ra, dec)
 
+    def getRefDelay (self, time):
+        """Get the reference direction used by the station beamformer.
+        Returns an ITRF vector in meters (numpy array of 3 floats).
+ 
+        `time`
+          Time at which to evaluate the direction
+        """
+        return self._response.getRefDelay(time)
+
     def setRefTile (self, ra, dec):
         """Set the reference direction used by the analog tile beamformer
         (relevant for HBA observations only). By default, LOFAR_TILE_BEAM_DIR
@@ -95,6 +104,16 @@ class stationresponse(object):
         """
         self._response.setRefTile(ra, dec)
 
+    def getRefTile (self, time):
+        """Get the reference direction used by the analog tile beamformer
+        (relevant for HBA observations only).
+        Returns an ITRF vector in meters (numpy array of 3 floats).
+ 
+        `time`
+          Time at which to evaluate the direction
+        """
+        return self._response.getRefTile(time)
+
     def setDirection (self, ra, dec):
         """Set the direction of interest (can be and often will be different
         from the pointing). By default, PHASE_DIR of field 0 is used.
@@ -106,6 +125,15 @@ class stationresponse(object):
         """
         self._response.setDirection(ra, dec)
 
+    def getDirection (self, time):
+        """Get the direction of interest.
+        Returns an ITRF vector in meters (numpy array of 3 floats).
+
+        `time`
+          Time at which to evaluate the direction
+        """
+        return self._response.getDirection(time)
+
     def evaluate (self, time):
         """Compute the beam Jones matrix for all stations and channels at the
         given time. The result is returned as a 4-dim complex numpy array with
@@ -158,3 +186,25 @@ class stationresponse(object):
           Frequency to compute beam at (in Hz)
         """
         return self._response.evaluate3(time, station, freq)
+
+    def evaluateFreqITRF (self, time, station, freq, direction, station0, tile0):
+        """Compute the beam Jones matrix for the given time, station, and
+        frequency, with the given ITRF directions.
+        The result is returned as a 2-dim complex numpy array with
+        shape: 2 x 2.
+
+        `time`
+          Time (MJD in seconds).
+        `station`
+          Station number (as defined in the ANTENNA table of the Measurement
+          Set).
+        `frequency`
+          Frequency to compute beam at (in Hz)
+        `direction`
+          ITRF direction to compute beam at (numpy array with 3 floats)
+        `station0`
+          ITRF direction of the station beamformer (numpy array with 3 floats)
+        `tile0`
+          ITRF direction of the tile beamformer (numpy array with 3 floats)
+        """
+        return self._response.evaluate4(time, station, freq, direction, station0, tile0)
diff --git a/CEP/Calibration/pystationresponse/src/pystationresponse.cc b/CEP/Calibration/pystationresponse/src/pystationresponse.cc
index 6a8b3f032e8c1ac93c9c748e44202bebc42e807a..e8f29bbc21c1cd59445496438ab5aa59b0826d33 100755
--- a/CEP/Calibration/pystationresponse/src/pystationresponse.cc
+++ b/CEP/Calibration/pystationresponse/src/pystationresponse.cc
@@ -171,23 +171,41 @@ namespace BBS
     // direction is the direction used by the station beamformer.
     void setRefDelay(double ra, double dec);
 
+    // Get the delay reference direction in meters, ITRF. The delay reference
+    // direction is the direction used by the station beamformer.
+    ValueHolder getRefDelay(real_t time);
+
     // Set the tile reference direction in radians, J2000. The tile reference
     // direction is the direction used by the analog tile beamformer and is
     // relevant only for HBA observations.
     void setRefTile(double ra, double dec);
 
+    // Get the tile reference direction in meters, ITRF. The delay reference
+    // direction is the direction used by the analog tile beamformer and is
+    // relevant only for HBA observations.
+    ValueHolder getRefTile(real_t time);
+
     // Set the direction of interest in radians, J2000. Can and often will be
     // different than the delay and/or tile reference direction.
     void setDirection(double ra, double dec);
 
+    // Get the direction of intereset in meters, ITRF.
+    ValueHolder getDirection(real_t time);
+
     // Compute the LOFAR beam Jones matrices for the given time, station, and/or
     // channel.
     ValueHolder evaluate0(double time);
     ValueHolder evaluate1(double time, int station);
     ValueHolder evaluate2(double time, int station, int channel);
     ValueHolder evaluate3(double time, int station, double freq);
+    ValueHolder evaluate4(double time, int station, double freq, const ValueHolder& direction, const ValueHolder& station0, const ValueHolder& tile0); 
 
   private:
+    Matrix<DComplex> evaluate_itrf(
+      const Station::ConstPtr &station, double time, double freq, double freq0,
+      const vector3r_t &direction, const vector3r_t &station0, 
+      const vector3r_t &tile0) const;
+
     Matrix<DComplex> evaluate(const Station::ConstPtr &station, double time,
       double freq,  double freq0) const;
 
@@ -280,18 +298,45 @@ namespace BBS
     itsRefDelay.reset(new ITRFDirection(itsRefPosition, direction));
   }
 
+  ValueHolder PyStationResponse::getRefDelay(real_t time)
+  {
+    vector3r_t refDelay=itsRefDelay->at(time);
+    Vector<Double> result(3);
+    result(0)=refDelay[0]; result(1)=refDelay[1]; result(2)=refDelay[2];
+
+    return ValueHolder(result);
+  }
+
   void PyStationResponse::setRefTile(double ra, double dec)
   {
     vector2r_t direction = {{ra, dec}};
     itsRefTile.reset(new ITRFDirection(itsRefPosition, direction));
   }
 
+  ValueHolder PyStationResponse::getRefTile(real_t time)
+  {
+    vector3r_t refTile=itsRefTile->at(time);
+    Vector<Double> result(3);
+    result(0)=refTile[0]; result(1)=refTile[1]; result(2)=refTile[2];
+
+    return ValueHolder(result);
+  }
+
   void PyStationResponse::setDirection(double ra, double dec)
   {
     vector2r_t direction = {{ra, dec}};
     itsDirection.reset(new ITRFDirection(itsRefPosition, direction));
   }
 
+  ValueHolder PyStationResponse::getDirection(real_t time)
+  {
+    vector3r_t direction=itsDirection->at(time);
+    Vector<Double> result(3);
+    result(0)=direction[0]; result(1)=direction[1]; result(2)=direction[2];
+
+    return ValueHolder(result);
+  }
+
   ValueHolder PyStationResponse::evaluate0(double time)
   {
     Array<DComplex> result(IPosition(4, 2, 2, itsChanFreq.size(),
@@ -364,6 +409,29 @@ namespace BBS
     return ValueHolder(evaluate(itsStations(station), time, freq, freq0));
   }
 
+  ValueHolder PyStationResponse::evaluate4(double time, int station, double freq, const ValueHolder& vh_direction, const ValueHolder& vh_station0, const ValueHolder& vh_tile0)
+  {
+    ASSERT (vh_direction.dataType() == TpArrayDouble);
+    ASSERT (vh_station0.dataType() == TpArrayDouble);
+    ASSERT (vh_tile0.dataType() == TpArrayDouble);
+    Array<Double> arr_dir(vh_direction.asArrayDouble());
+    Array<Double> st0_dir(vh_direction.asArrayDouble());
+    Array<Double> tile_dir(vh_direction.asArrayDouble());
+    vector3r_t direction={{arr_dir.data()[0],arr_dir.data()[1],arr_dir.data()[2]}};
+    vector3r_t station0 ={{st0_dir.data()[0],st0_dir.data()[1],st0_dir.data()[2]}};
+    vector3r_t tile0    ={{tile_dir.data()[0],tile_dir.data()[1],tile_dir.data()[2]}};
+
+    if(itsUseChanFreq)
+    {
+      return ValueHolder(evaluate_itrf(itsStations(station), time, freq, freq, 
+                                       direction, station0, tile0));
+    }
+
+    double freq0 = itsRefFreq(0);
+    return ValueHolder(evaluate_itrf(itsStations(station), time, freq, freq0,
+                                     direction, station0, tile0));
+  }
+
   Cube<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
     double time, const Vector<Double> &freq, const Vector<Double> &freq0) const
   {
@@ -445,16 +513,14 @@ namespace BBS
     return result;
   }
 
-  Matrix<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
-    double time, double freq,  double freq0) const
+  Matrix<DComplex> PyStationResponse::evaluate_itrf(
+    const Station::ConstPtr &station, double time, double freq, double freq0,
+    const vector3r_t &direction, const vector3r_t &station0, 
+    const vector3r_t &tile0) const
   {
     Matrix<DComplex> result(2, 2, 0.0);
     if(itsUseArrayFactor)
     {
-      vector3r_t direction = itsDirection->at(time);
-      vector3r_t station0 = itsRefDelay->at(time);
-      vector3r_t tile0 = itsRefTile->at(time);
-
       if(itsUseElementResponse)
       {
         matrix22c_t response = station->response(time, freq, direction, freq0,
@@ -491,7 +557,6 @@ namespace BBS
       // the station is always selected.
       AntennaField::ConstPtr field = *station->beginFields();
 
-      vector3r_t direction = itsDirection->at(time);
       matrix22c_t response = field->elementResponse(time, freq,
         direction);
 
@@ -514,6 +579,22 @@ namespace BBS
     return result;
   }
 
+  Matrix<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
+    double time, double freq,  double freq0) const
+  {
+    vector3r_t direction;
+    vector3r_t station0;
+    vector3r_t tile0;
+    if (itsUseArrayFactor) {
+      direction = itsDirection->at(time);
+      station0 = itsRefDelay->at(time);
+      tile0 = itsRefTile->at(time);
+    } else if (itsUseElementResponse) {
+      direction = itsDirection->at(time);
+    }
+    return evaluate_itrf(station, time, freq, freq0, direction, station0, tile0);
+  }
+
   void PyStationResponse::invert(matrix22c_t &in) const
   {
     complex_t invDet = 1.0 / (in[0][0] * in[1][1] - in[0][1] * in[1][0]);
@@ -548,10 +629,16 @@ namespace BBS
         (boost::python::arg("type")="other"))
       .def ("setRefDelay", &PyStationResponse::setRefDelay,
         (boost::python::arg("ra"), boost::python::arg("dec")))
+      .def ("getRefDelay", &PyStationResponse::getRefDelay,
+        (boost::python::arg("time")))
       .def ("setRefTile", &PyStationResponse::setRefTile,
         (boost::python::arg("ra"), boost::python::arg("dec")))
+      .def ("getRefTile", &PyStationResponse::getRefTile,
+        (boost::python::arg("time")))
       .def ("setDirection", &PyStationResponse::setDirection,
         (boost::python::arg("ra"), boost::python::arg("dec")))
+      .def ("getDirection", &PyStationResponse::getDirection,
+        (boost::python::arg("time")))
       .def ("evaluate0", &PyStationResponse::evaluate0,
         (boost::python::arg("time")))
       .def ("evaluate1", &PyStationResponse::evaluate1,
@@ -562,6 +649,9 @@ namespace BBS
       .def ("evaluate3", &PyStationResponse::evaluate3,
       (boost::python::arg("time"), boost::python::arg("station"),
          boost::python::arg("freq")))
+      .def ("evaluate4", &PyStationResponse::evaluate4,
+      (boost::python::arg("time"), boost::python::arg("station"),
+         boost::python::arg("freq"), boost::python::arg("direction")))
       ;
   }