diff --git a/StationResponse/CMakeLists.txt b/StationResponse/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..16350a206a86c1e877f491e6a2e865ddebdf9386
--- /dev/null
+++ b/StationResponse/CMakeLists.txt
@@ -0,0 +1,14 @@
+# $Id$
+
+lofar_package(StationResponse 0.1 DEPENDS Common ElementResponse)
+
+include(LofarFindPackage)
+lofar_find_package(Boost REQUIRED)
+lofar_find_package(Casacore REQUIRED COMPONENTS casa measures ms tables images coordinates)
+
+# Uncomment to check for unsafe conversions (gcc), for example conversion of
+# size_t to unsigned int (truncation).
+#add_definitions(-Wconversion)
+
+add_subdirectory(include/StationResponse)
+add_subdirectory(src)
diff --git a/StationResponse/include/StationResponse/AntennaField.h b/StationResponse/include/StationResponse/AntennaField.h
new file mode 100644
index 0000000000000000000000000000000000000000..934e0b2247251b8628449391c62277ffb526a319
--- /dev/null
+++ b/StationResponse/include/StationResponse/AntennaField.h
@@ -0,0 +1,261 @@
+//# AntennaField.h: Representation of a LOFAR antenna field, with methods to
+//# compute its response to incoming radiation.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_ARRAYFIELD_H
+#define LOFAR_STATIONRESPONSE_ARRAYFIELD_H
+
+// \file
+// Representation of a LOFAR antenna field, with methods to compute its response
+// to incoming radiation.
+
+#include <Common/lofar_smartptr.h>
+#include <Common/lofar_string.h>
+#include <Common/lofar_vector.h>
+#include <StationResponse/Constants.h>
+#include <StationResponse/Types.h>
+#include <StationResponse/ITRFDirection.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+/**
+ * \brief Base class that represents a generic LOFAR antenna field.
+ */
+class AntennaField
+{
+public:
+ typedef shared_ptr<AntennaField> Ptr;
+ typedef shared_ptr<const AntennaField> ConstPtr;
+
+ /**
+ * \brief Antenna field coordinate system.
+ *
+ * A right handed, cartesian, local coordinate system with coordinate axes
+ * \p p, \p q, and \p r is associated with each antenna field.
+ *
+ * The r-axis is orthogonal to the antenna field, and points towards the
+ * local pseudo zenith.
+ *
+ * The q-axis is the northern bisector of the \p X and \p Y dipoles, i.e.
+ * it is the reference direction from which the orientation of the dual
+ * dipole antennae is determined. The q-axis points towards the North at
+ * the core. At remote sites it is defined as the intersection of the
+ * antenna field plane and a plane parallel to the meridian plane at the
+ * core. This ensures the reference directions at all sites are similar.
+ *
+ * The p-axis is orthogonal to both other axes, and points towards the East
+ * at the core.
+ *
+ * The axes and origin of the anntena field coordinate system are expressed
+ * as vectors in the geocentric, cartesian, ITRF coordinate system, in
+ * meters.
+ *
+ * \sa "LOFAR Reference Plane and Reference Direction", M.A. Brentjens,
+ * LOFAR-ASTRON-MEM-248.
+ */
+ struct CoordinateSystem
+ {
+ struct Axes
+ {
+ vector3r_t p;
+ vector3r_t q;
+ vector3r_t r;
+ };
+
+ vector3r_t origin;
+ Axes axes;
+ };
+
+ /** A single antenna. */
+ struct Antenna
+ {
+ /**
+ * \brief Position of the antenna relative to the antenna field center
+ * (origin). This is a vector expressed in the geocentric, cartesian,
+ * ITRF coordinate system, in meters.
+ */
+ vector3r_t position;
+
+ /**
+ * \brief Status of the \p X and \p Y signal paths of the antenna,
+ * respectively.
+ */
+ bool enabled[2];
+ };
+
+ typedef vector<Antenna> AntennaList;
+
+
+ AntennaField(const string &name, const CoordinateSystem &coordinates);
+
+ virtual ~AntennaField();
+
+ /** Return the name of the antenna field. */
+ const string &name() const;
+
+ /** Return the phase reference position of the antenna field. */
+ const vector3r_t &position() const;
+
+ /** Return the antenna field coordinate system. */
+ const CoordinateSystem &coordinates() const;
+
+ /** Add an antenna to the antenna field. */
+ void addAntenna(const Antenna &antenna);
+
+ /** Return the number of antennae in the antenna field. */
+ size_t nAntennae() const;
+
+ /*!
+ * \name Antennae accessors
+ * These member functions provide access to the antennae that are part of
+ * the antenna field.
+ */
+ // @{
+
+ /** Return a read-only reference to the antenna with the requested index. */
+ const Antenna &antenna(size_t n) const;
+
+ /** Return a writeable reference to the antenna with the requested index. */
+ Antenna &antenna(size_t n);
+
+ /** Return a read-only iterator that points to the first antenna of the
+ * antenna field.
+ */
+ AntennaList::const_iterator beginAntennae() const;
+
+ /** Return a read-only iterator that points one position past the last
+ * antenna of the antenna field.
+ */
+ AntennaList::const_iterator endAntennae() const;
+
+ // @}
+
+ /*!
+ * \brief Compute the response of the antenna field for a plane wave of
+ * frequency \p freq, arriving from direction \p direction, with the analog
+ * %tile beam former steered towards \p direction0. For LBA antenna fields,
+ * \p direction0 has no effect.
+ *
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Frequency of the plane wave (Hz).
+ * \param direction Direction of arrival (ITRF, m).
+ * \param direction0 Tile beam former reference direction (ITRF, m).
+ * \return Jones matrix that represents the response of the antenna field.
+ *
+ * The directions \p direction, and \p direction0 are vectors that
+ * represent a direction of \e arrival. These vectors have unit length and
+ * point \e from the ground \e towards the direction from which the plane
+ * wave arrives.
+ */
+ virtual matrix22c_t response(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ /*!
+ * \brief Compute the array factor of the antenna field for a plane wave of
+ * frequency \p freq, arriving from direction \p direction, analog %tile
+ * beam former steered towards \p direction0. For LBA antenna fields,
+ * \p direction0 has no effect.
+ *
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Frequency of the plane wave (Hz).
+ * \param direction Direction of arrival (ITRF, m).
+ * \param direction0 Tile beam former reference direction (ITRF, m).
+ * \return A diagonal matrix with the array factor of the X and Y antennae.
+ *
+ * The directions \p direction, and \p direction0 are vectors that
+ * represent a direction of \e arrival. These vectors have unit length and
+ * point \e from the ground \e towards the direction from which the plane
+ * wave arrives.
+ */
+ virtual diag22c_t arrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ /*!
+ * \brief Compute the response of the antenna field for a plane wave of
+ * frequency \p freq, arriving from direction \p direction, with the analog
+ * %tile beam former steered towards \p direction0. For LBA antenna fields,
+ * \p direction0 has no effect.
+ *
+ * This method returns the non-normalized (raw) response. This allows the
+ * response of several antenna fields to be summed together, followed by
+ * normalization of the sum.
+ *
+ * \see response(real_t time, real_t freq, const vector3r_t &direction,
+ * const vector3r_t &direction0) const
+ */
+ virtual raw_response_t rawResponse(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ /*!
+ * \brief Compute the array factor of the antenna field for a plane wave of
+ * frequency \p freq, arriving from direction \p direction, analog %tile
+ * beam former steered towards \p direction0. For LBA antenna fields,
+ * \p direction0 has no effect.
+ *
+ * This method returns the non-normalized (raw) array factor. This allows
+ * the array factor of several antenna fields to be summed together,
+ * followed by normalization of the sum.
+ *
+ * \see diag22c_t arrayFactor(real_t time, real_t freq,
+ * const vector3r_t &direction, const vector3r_t &direction0) const
+ */
+ virtual raw_array_factor_t rawArrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const = 0;
+
+ /*!
+ * \brief Compute the response of a single antenna for a plane wave of
+ * frequency \p freq, arriving from direction \p direction.
+ *
+ */
+ virtual matrix22c_t elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction) const = 0;
+
+protected:
+ /** Compute the parallactic rotation. */
+ matrix22r_t rotation(real_t time, const vector3r_t &direction) const;
+
+ /** Transform a vector from ITRF to antenna field coordinates. */
+ vector3r_t itrf2field(const vector3r_t &itrf) const;
+
+private:
+ vector3r_t ncp(real_t time) const;
+
+ string itsName;
+ CoordinateSystem itsCoordinateSystem;
+ AntennaList itsAntennae;
+ ITRFDirection::Ptr itsNCP;
+ mutable real_t itsNCPCacheTime;
+ mutable vector3r_t itsNCPCacheDirection;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/AntennaFieldHBA.h b/StationResponse/include/StationResponse/AntennaFieldHBA.h
new file mode 100644
index 0000000000000000000000000000000000000000..6caafee496441b05685dafa17c97a54e51bdd7a4
--- /dev/null
+++ b/StationResponse/include/StationResponse/AntennaFieldHBA.h
@@ -0,0 +1,73 @@
+//# AntennaFieldHBA.h: Representation of an HBA antenna field.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_ANTENNAFIELDHBA_H
+#define LOFAR_STATIONRESPONSE_ANTENNAFIELDHBA_H
+
+// \file
+// Representation of an HBA antenna field.
+
+#include <StationResponse/AntennaField.h>
+#include <StationResponse/AntennaModelHBA.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class AntennaFieldHBA: public AntennaField
+{
+public:
+ typedef shared_ptr<AntennaFieldHBA> Ptr;
+ typedef shared_ptr<const AntennaFieldHBA> ConstPtr;
+
+ AntennaFieldHBA(const string &name, const CoordinateSystem &coordinates,
+ const AntennaModelHBA::ConstPtr &model);
+
+ virtual matrix22c_t response(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ virtual diag22c_t arrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ virtual raw_response_t rawResponse(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ virtual raw_array_factor_t rawArrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ virtual matrix22c_t elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction) const;
+
+private:
+ AntennaModelHBA::ConstPtr itsAntennaModel;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/AntennaFieldLBA.h b/StationResponse/include/StationResponse/AntennaFieldLBA.h
new file mode 100644
index 0000000000000000000000000000000000000000..f4c082542910e7bf94109bed6ab3ff3c828d8c94
--- /dev/null
+++ b/StationResponse/include/StationResponse/AntennaFieldLBA.h
@@ -0,0 +1,64 @@
+//# AntennaFieldLBA.h: Representation of an LBA antenna field.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_ANTENNAFIELDLBA_H
+#define LOFAR_STATIONRESPONSE_ANTENNAFIELDLBA_H
+
+// \file
+// Representation of an LBA antenna field.
+
+#include <StationResponse/AntennaField.h>
+#include <StationResponse/AntennaModelLBA.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class AntennaFieldLBA: public AntennaField
+{
+public:
+ typedef shared_ptr<AntennaFieldLBA> Ptr;
+ typedef shared_ptr<const AntennaFieldLBA> ConstPtr;
+
+ AntennaFieldLBA(const string &name, const CoordinateSystem &coordinates,
+ const AntennaModelLBA::ConstPtr &model);
+
+ virtual raw_array_factor_t rawArrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ virtual matrix22c_t elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction) const;
+
+private:
+ AntennaModelLBA::ConstPtr itsAntennaModel;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/AntennaModelHBA.h b/StationResponse/include/StationResponse/AntennaModelHBA.h
new file mode 100644
index 0000000000000000000000000000000000000000..9ef452c0d5be39b56da9b69004a737018dfdca8d
--- /dev/null
+++ b/StationResponse/include/StationResponse/AntennaModelHBA.h
@@ -0,0 +1,69 @@
+//# AntennaModelHBA.h: HBA antenna model interface definitions.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_ANTENNAMODELHBA_H
+#define LOFAR_STATIONRESPONSE_ANTENNAMODELHBA_H
+
+// \file
+// HBA antenna model interface definitions.
+
+#include <StationResponse/Types.h>
+#include <Common/lofar_smartptr.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class AntennaModelHBA
+{
+public:
+ typedef shared_ptr<AntennaModelHBA> Ptr;
+ typedef shared_ptr<const AntennaModelHBA> ConstPtr;
+
+ virtual ~AntennaModelHBA();
+
+ virtual matrix22c_t response(real_t freq, const vector3r_t &direction,
+ const vector3r_t &direction0) const;
+
+ virtual diag22c_t arrayFactor(real_t freq, const vector3r_t &direction,
+ const vector3r_t &direction0) const;
+
+ virtual raw_response_t rawResponse(real_t freq, const vector3r_t &direction,
+ const vector3r_t &direction0) const;
+
+ virtual raw_array_factor_t rawArrayFactor(real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const = 0;
+
+ virtual matrix22c_t elementResponse(real_t freq,
+ const vector3r_t &direction) const = 0;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/AntennaModelLBA.h b/StationResponse/include/StationResponse/AntennaModelLBA.h
new file mode 100644
index 0000000000000000000000000000000000000000..7b6bfa543c633284dd52813866bbe6fa9c2ac9d4
--- /dev/null
+++ b/StationResponse/include/StationResponse/AntennaModelLBA.h
@@ -0,0 +1,57 @@
+//# AntennaModelLBA.h: LBA antenna model interface definitions.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_ANTENNAMODELLBA_H
+#define LOFAR_STATIONRESPONSE_ANTENNAMODELLBA_H
+
+// \file
+// LBA antenna model interface definitions.
+
+#include <StationResponse/Types.h>
+#include <Common/lofar_smartptr.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class AntennaModelLBA
+{
+public:
+ typedef shared_ptr<AntennaModelLBA> Ptr;
+ typedef shared_ptr<const AntennaModelLBA> ConstPtr;
+
+ virtual ~AntennaModelLBA();
+
+ virtual matrix22c_t response(real_t freq, const vector3r_t &direction)
+ const = 0;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/CMakeLists.txt b/StationResponse/include/StationResponse/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b4dc8dff3c3c708ec21a53a78bc8c087148742f7
--- /dev/null
+++ b/StationResponse/include/StationResponse/CMakeLists.txt
@@ -0,0 +1,23 @@
+# $Id$
+
+# Create symbolic link to include directory.
+execute_process(COMMAND ${CMAKE_COMMAND} -E create_symlink
+ ${CMAKE_CURRENT_SOURCE_DIR}
+ ${CMAKE_BINARY_DIR}/include/${PACKAGE_NAME})
+
+# Install header files.
+install(FILES
+ AntennaField.h
+ AntennaFieldHBA.h
+ AntennaFieldLBA.h
+ AntennaModelHBA.h
+ AntennaModelLBA.h
+ Constants.h
+ DualDipoleAntenna.h
+ ITRFDirection.h
+ LofarMetaDataUtil.h
+ MathUtil.h
+ Station.h
+ TileAntenna.h
+ Types.h
+ DESTINATION include/${PACKAGE_NAME})
diff --git a/StationResponse/include/StationResponse/Constants.h b/StationResponse/include/StationResponse/Constants.h
new file mode 100644
index 0000000000000000000000000000000000000000..c77be6d22cdbe63369c4f9cc3a6c9d45745d62a4
--- /dev/null
+++ b/StationResponse/include/StationResponse/Constants.h
@@ -0,0 +1,60 @@
+//# Constants.h: %Constants used in this library.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_CONSTANTS_H
+#define LOFAR_STATIONRESPONSE_CONSTANTS_H
+
+// \file
+// %Constants used in this library.
+
+#include <StationResponse/Types.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+/** %Constants used in this library. */
+namespace Constants
+{
+/** 2.0 * pi */
+const real_t _2pi = 6.283185307179586476925286;
+
+/** pi / 2.0 */
+const real_t pi_2 = 1.570796326794896619231322;
+
+/** pi / 4.0 */
+const real_t pi_4 = 0.7853981633974483096156608;
+
+/** Speed of light (m/s) */
+const real_t c = 2.99792458e+08;
+} //# namespace Constants
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/DualDipoleAntenna.h b/StationResponse/include/StationResponse/DualDipoleAntenna.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ac6439fc6e7572464a4f29cb72362b2d3d5df96
--- /dev/null
+++ b/StationResponse/include/StationResponse/DualDipoleAntenna.h
@@ -0,0 +1,55 @@
+//# DualDipoleAntenna.h: Semi-analytical model of a LOFAR LBA dual dipole
+//# antenna.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_DUALDIPOLEANTENNA_H
+#define LOFAR_STATIONRESPONSE_DUALDIPOLEANTENNA_H
+
+// \file
+// Semi-analytical model of a LOFAR LBA dual dipole antenna.
+
+#include <StationResponse/AntennaModelLBA.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class DualDipoleAntenna: public AntennaModelLBA
+{
+public:
+ typedef shared_ptr<DualDipoleAntenna> Ptr;
+ typedef shared_ptr<const DualDipoleAntenna> ConstPtr;
+
+ virtual matrix22c_t response(real_t freq, const vector3r_t &direction)
+ const;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/ITRFDirection.h b/StationResponse/include/StationResponse/ITRFDirection.h
new file mode 100644
index 0000000000000000000000000000000000000000..280a6111a99443ddbe180c8eb45bdd5a72906d45
--- /dev/null
+++ b/StationResponse/include/StationResponse/ITRFDirection.h
@@ -0,0 +1,65 @@
+//# ITRFDirection.h: Functor that maps time to an ITRF direction.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_ITRFDIRECTION_H
+#define LOFAR_STATIONRESPONSE_ITRFDIRECTION_H
+
+// \file
+// Functor that maps time to an ITRF direction.
+
+#include <StationResponse/Types.h>
+#include <Common/lofar_smartptr.h>
+
+#include <casacore/measures/Measures/MeasFrame.h>
+#include <casacore/measures/Measures/MeasConvert.h>
+#include <casacore/measures/Measures/MCDirection.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class ITRFDirection
+{
+public:
+ typedef shared_ptr<ITRFDirection> Ptr;
+ typedef shared_ptr<const ITRFDirection> ConstPtr;
+
+ ITRFDirection(const vector3r_t &position, const vector2r_t &direction);
+ ITRFDirection(const vector3r_t &position, const vector3r_t &direction);
+
+ vector3r_t at(real_t time) const;
+
+private:
+ mutable casacore::MeasFrame itsFrame;
+ mutable casacore::MDirection::Convert itsConverter;
+};
+
+// @}
+
+} //# namespace BBS
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/LofarMetaDataUtil.h b/StationResponse/include/StationResponse/LofarMetaDataUtil.h
new file mode 100644
index 0000000000000000000000000000000000000000..57327b0ed0eb3c0bfbdd7857ef71b4335822a401
--- /dev/null
+++ b/StationResponse/include/StationResponse/LofarMetaDataUtil.h
@@ -0,0 +1,65 @@
+//# LofarMetaDataUtil.h: Utility functions to read the meta data relevant for
+//# simulating the beam from LOFAR observations stored in MS format.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_LOFARMETADATAUTIL_H
+#define LOFAR_STATIONRESPONSE_LOFARMETADATAUTIL_H
+
+// \file
+// Utility functions to read the meta data relevant for simulating the beam from
+// LOFAR observations stored in MS format.
+
+#include <StationResponse/Station.h>
+#include <casacore/ms/MeasurementSets/MeasurementSet.h>
+#include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
+#include <casacore/measures/Measures/MDirection.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+Station::Ptr readStation(const casacore::MeasurementSet &ms, unsigned int id);
+
+template <typename T>
+void readStations(const casacore::MeasurementSet &ms, T out_it)
+{
+ casacore::ROMSAntennaColumns antenna(ms.antenna());
+ for(unsigned int i = 0; i < antenna.nrow(); ++i)
+ {
+ *out_it++ = readStation(ms, i);
+ }
+}
+
+// Read the tile beam direction from a LOFAR MS. If it is not defined,
+// this function returns the delay center.
+casacore::MDirection readTileBeamDirection(const casacore::MeasurementSet &ms);
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/MathUtil.h b/StationResponse/include/StationResponse/MathUtil.h
new file mode 100644
index 0000000000000000000000000000000000000000..cb16fb898c2e0f9e229b7253823d53525fcb502a
--- /dev/null
+++ b/StationResponse/include/StationResponse/MathUtil.h
@@ -0,0 +1,172 @@
+//# MathUtil.h: Various mathematical operations on vectors and matrices.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_MATHUTIL_H
+#define LOFAR_STATIONRESPONSE_MATHUTIL_H
+
+// \file
+// Various mathematical operations on vectors and matrices.
+
+#include <StationResponse/Constants.h>
+#include <StationResponse/Types.h>
+#include <Common/lofar_math.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+inline real_t dot(const vector3r_t &arg0, const vector3r_t &arg1)
+{
+ return arg0[0] * arg1[0] + arg0[1] * arg1[1] + arg0[2] * arg1[2];
+}
+
+inline double norm(const vector3r_t &arg0)
+{
+ return sqrt(dot(arg0, arg0));
+}
+
+inline vector3r_t operator*(real_t arg0, const vector3r_t arg1)
+{
+ vector3r_t result = {{arg0 * arg1[0], arg0 * arg1[1], arg0 * arg1[2]}};
+ return result;
+}
+
+inline vector3r_t normalize(const vector3r_t &arg0)
+{
+ return 1.0 / norm(arg0) * arg0;
+}
+
+inline vector2r_t cart2thetaphi(const vector3r_t &cart)
+{
+ real_t r = sqrt(cart[0] * cart[0] + cart[1] * cart[1]);
+ vector2r_t thetaphi = {{Constants::pi_2 - atan2(cart[2], r), atan2(cart[1],
+ cart[0])}};
+ return thetaphi;
+}
+
+inline vector3r_t thetaphi2cart(const vector2r_t &thetaphi)
+{
+ real_t r = sin(thetaphi[0]);
+ vector3r_t cart = {{r * cos(thetaphi[1]), r * sin(thetaphi[1]),
+ cos(thetaphi[0])}};
+ return cart;
+}
+
+// returns az, el, r.
+inline vector3r_t cart2sph(const vector3r_t &cart)
+{
+ real_t r = sqrt(cart[0] * cart[0] + cart[1] * cart[1]);
+
+ vector3r_t sph;
+ sph[0] = atan2(cart[1], cart[0]);
+ sph[1] = atan2(cart[2], r);
+ sph[2] = norm(cart);
+ return sph;
+}
+
+// expects az, el, r.
+inline vector3r_t sph2cart(const vector3r_t &sph)
+{
+ vector3r_t cart = {{sph[2] * cos(sph[1]) * cos(sph[0]), sph[2] * cos(sph[1])
+ * sin(sph[0]), sph[2] * sin(sph[1])}};
+ return cart;
+}
+
+inline matrix22c_t operator*(const matrix22c_t &arg0, const matrix22r_t &arg1)
+{
+ matrix22c_t result;
+ result[0][0] = arg0[0][0] * arg1[0][0] + arg0[0][1] * arg1[1][0];
+ result[0][1] = arg0[0][0] * arg1[0][1] + arg0[0][1] * arg1[1][1];
+ result[1][0] = arg0[1][0] * arg1[0][0] + arg0[1][1] * arg1[1][0];
+ result[1][1] = arg0[1][0] * arg1[0][1] + arg0[1][1] * arg1[1][1];
+ return result;
+}
+
+inline vector3r_t cross(const vector3r_t &arg0, const vector3r_t &arg1)
+{
+ vector3r_t result;
+ result[0] = arg0[1] * arg1[2] - arg0[2] * arg1[1];
+ result[1] = arg0[2] * arg1[0] - arg0[0] * arg1[2];
+ result[2] = arg0[0] * arg1[1] - arg0[1] * arg1[0];
+ return result;
+}
+
+inline vector3r_t operator+(const vector3r_t &arg0, const vector3r_t &arg1)
+{
+ vector3r_t result = {{arg0[0] + arg1[0], arg0[1] + arg1[1],
+ arg0[2] + arg1[2]}};
+ return result;
+}
+
+inline vector3r_t operator-(const vector3r_t &arg0, const vector3r_t &arg1)
+{
+ vector3r_t result = {{arg0[0] - arg1[0], arg0[1] - arg1[1],
+ arg0[2] - arg1[2]}};
+ return result;
+}
+
+inline matrix22c_t normalize(const raw_response_t &raw)
+{
+ matrix22c_t response = {{ {{}}, {{}} }};
+
+ if(raw.weight[0] != 0.0)
+ {
+ response[0][0] = raw.response[0][0] / raw.weight[0];
+ response[0][1] = raw.response[0][1] / raw.weight[0];
+ }
+
+ if(raw.weight[1] != 0.0)
+ {
+ response[1][0] = raw.response[1][0] / raw.weight[1];
+ response[1][1] = raw.response[1][1] / raw.weight[1];
+ }
+
+ return response;
+}
+
+inline diag22c_t normalize(const raw_array_factor_t &raw)
+{
+ diag22c_t af = {{}};
+
+ if(raw.weight[0] != 0.0)
+ {
+ af[0] = raw.factor[0] / raw.weight[0];
+ }
+
+ if(raw.weight[1] != 0.0)
+ {
+ af[1] = raw.factor[1] / raw.weight[1];
+ }
+
+ return af;
+}
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/Station.h b/StationResponse/include/StationResponse/Station.h
new file mode 100644
index 0000000000000000000000000000000000000000..0019c33dcdf0c9577e4b0e9ec9184096812ab710
--- /dev/null
+++ b/StationResponse/include/StationResponse/Station.h
@@ -0,0 +1,361 @@
+//# Station.h: Representation of the station beam former.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_STATION_H
+#define LOFAR_STATIONRESPONSE_STATION_H
+
+// \file
+// Representation of the station beam former.
+
+#include <Common/lofar_smartptr.h>
+#include <Common/lofar_string.h>
+#include <Common/lofar_vector.h>
+#include <StationResponse/AntennaField.h>
+#include <StationResponse/Types.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class Station
+{
+public:
+ typedef shared_ptr<Station> Ptr;
+ typedef shared_ptr<const Station> ConstPtr;
+ typedef vector<AntennaField::ConstPtr> FieldList;
+
+ /*!
+ * \brief Construct a new Station instance.
+ *
+ * \param name Name of the station.
+ * \param position Position of the station (ITRF, m).
+ */
+ Station(const string &name, const vector3r_t &position);
+
+ /*!
+ * \brief Return the name of the station.
+ */
+ const string &name() const;
+
+ /*!
+ * \brief Return the position of the station (ITRF, m).
+ */
+ const vector3r_t &position() const;
+
+ /*!
+ * \brief Set the phase reference position. This is the position where the
+ * delay of the incoming plane wave is assumed to be zero.
+ *
+ * \param reference Phase reference position (ITRF, m).
+ *
+ * By default, it is assumed the position of the station is also the phase
+ * reference position. Use this method to set the phase reference position
+ * explicitly when this assumption is false.
+ */
+ void setPhaseReference(const vector3r_t &reference);
+
+ /*!
+ * \brief Return the phase reference position (ITRF, m).
+ *
+ * \see Station::setPhaseReference()
+ */
+ const vector3r_t &phaseReference() const;
+
+ /*!
+ * \brief Add an antenna field to the station.
+ *
+ * Physical %LOFAR stations consist of an LBA field, and either one (remote
+ * and international stations) or two (core stations) HBA fields. Virtual
+ * %LOFAR stations can consist of a combination of the antenna fields of
+ * several physical stations.
+ *
+ * Use this method to add the appropriate antenna fields to the station.
+ */
+ void addField(const AntennaField::ConstPtr &field);
+
+
+ /*!
+ * \brief Return the number of available antenna fields.
+ */
+ size_t nFields() const;
+
+ /*!
+ * \brief Return the requested antenna field.
+ *
+ * \param i Antenna field number (0-based).
+ * \return An AntennaField::ConstPtr to the requested AntennaField
+ * instance, or an empty AntennaField::ConstPtr if \p i is out of bounds.
+ */
+ AntennaField::ConstPtr field(size_t i) const;
+
+ /*!
+ * \brief Return an iterator that points to the beginning of the list of
+ * antenna fields.
+ */
+ FieldList::const_iterator beginFields() const;
+
+ /*!
+ * \brief Return an iterator that points to the end of the list of antenna
+ * fields.
+ */
+ FieldList::const_iterator endFields() const;
+
+ /*!
+ * \brief Compute the response of the station for a plane wave of frequency
+ * \p freq, arriving from direction \p direction, with the %station beam
+ * former steered towards \p station0, and, for HBA stations, the analog
+ * %tile beam former steered towards \p tile0. For LBA stations, \p tile0
+ * has no effect.
+ *
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Frequency of the plane wave (Hz).
+ * \param direction Direction of arrival (ITRF, m).
+ * \param freq0 %Station beam former reference frequency (Hz).
+ * \param station0 %Station beam former reference direction (ITRF, m).
+ * \param tile0 Tile beam former reference direction (ITRF, m).
+ * \return Jones matrix that represents the %station response.
+ *
+ * For any given sub-band, the %LOFAR station beam former computes weights
+ * for a single reference frequency. Usually, this reference frequency is
+ * the center frequency of the sub-band. For any frequency except the
+ * reference frequency, these weights are an approximation. This aspect of
+ * the system is taken into account in the computation of the response.
+ * Therefore, both the frequency of interest \p freq and the reference
+ * frequency \p freq0 need to be specified.
+ *
+ * The directions \p direction, \p station0, and \p tile0 are vectors that
+ * represent a direction of \e arrival. These vectors have unit length and
+ * point \e from the ground \e towards the direction from which the plane
+ * wave arrives.
+ */
+ matrix22c_t response(real_t time, real_t freq, const vector3r_t &direction,
+ real_t freq0, const vector3r_t &station0, const vector3r_t &tile0)
+ const;
+
+ /*!
+ * \brief Compute the array factor of the station for a plane wave of
+ * frequency \p freq, arriving from direction \p direction, with the
+ * %station beam former steered towards \p station0, and, for HBA stations
+ * the analog %tile beam former steered towards \p tile0. For LBA stations,
+ * \p tile0 has no effect.
+ *
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Frequency of the plane wave (Hz).
+ * \param direction Direction of arrival (ITRF, m).
+ * \param freq0 %Station beam former reference frequency (Hz).
+ * \param station0 %Station beam former reference direction (ITRF, m).
+ * \param tile0 Tile beam former reference direction (ITRF, m).
+ * \return A diagonal matrix with the array factor of the X and Y antennae.
+ *
+ * For any given sub-band, the %LOFAR station beam former computes weights
+ * for a single reference frequency. Usually, this reference frequency is
+ * the center frequency of the sub-band. For any frequency except the
+ * reference frequency, these weights are an approximation. This aspect of
+ * the system is taken into account in the computation of the response.
+ * Therefore, both the frequency of interest \p freq and the reference
+ * frequency \p freq0 need to be specified.
+ *
+ * The directions \p direction, \p station0, and \p tile0 are vectors that
+ * represent a direction of \e arrival. These vectors have unit length and
+ * point \e from the ground \e towards the direction from which the plane
+ * wave arrives.
+ */
+ diag22c_t arrayFactor(real_t time, real_t freq, const vector3r_t &direction,
+ real_t freq0, const vector3r_t &station0, const vector3r_t &tile0)
+ const;
+
+ /*!
+ * \name Convenience member functions
+ * These member functions perform the same function as the corresponding
+ * non-template member functions, for a list of frequencies or (frequency,
+ * reference frequency) pairs.
+ */
+ // @{
+
+ /*!
+ * \brief Convenience method to compute the response of the station for a
+ * list of frequencies, and a fixed reference frequency.
+ *
+ * \param count Number of frequencies.
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Input iterator for a list of frequencies (Hz) of length
+ * \p count.
+ * \param direction Direction of arrival (ITRF, m).
+ * \param freq0 %Station beam former reference frequency (Hz).
+ * \param station0 %Station beam former reference direction (ITRF, m).
+ * \param tile0 Tile beam former reference direction (ITRF, m).
+ * \param buffer Output iterator with room for \p count instances of type
+ * ::matrix22c_t.
+ *
+ * \see response(real_t time, real_t freq, const vector3r_t &direction,
+ * real_t freq0, const vector3r_t &station0, const vector3r_t &tile0) const
+ */
+ template <typename T, typename U>
+ void response(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const;
+
+ /*!
+ * \brief Convenience method to compute the array factor of the station for
+ * a list of frequencies, and a fixed reference frequency.
+ *
+ * \param count Number of frequencies.
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Input iterator for a list of frequencies (Hz) of length
+ * \p count.
+ * \param direction Direction of arrival (ITRF, m).
+ * \param freq0 %Station beam former reference frequency (Hz).
+ * \param station0 %Station beam former reference direction (ITRF, m).
+ * \param tile0 Tile beam former reference direction (ITRF, m).
+ * \param buffer Output iterator with room for \p count instances of type
+ * ::diag22c_t.
+ *
+ * \see arrayFactor(real_t time, real_t freq, const vector3r_t &direction,
+ * real_t freq0, const vector3r_t &station0, const vector3r_t &tile0) const
+ */
+ template <typename T, typename U>
+ void arrayFactor(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const;
+
+ /*!
+ * \brief Convenience method to compute the response of the station for a
+ * list of (frequency, reference frequency) pairs.
+ *
+ * \param count Number of frequencies.
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Input iterator for a list of frequencies (Hz) of length
+ * \p count.
+ * \param direction Direction of arrival (ITRF, m).
+ * \param freq0 Input iterator for a list of %Station beam former reference
+ * frequencies (Hz) of length \p count.
+ * \param station0 %Station beam former reference direction (ITRF, m).
+ * \param tile0 Tile beam former reference direction (ITRF, m).
+ * \param buffer Output iterator with room for \p count instances of type
+ * ::matrix22c_t.
+ *
+ * \see response(real_t time, real_t freq, const vector3r_t &direction,
+ * real_t freq0, const vector3r_t &station0, const vector3r_t &tile0) const
+ */
+ template <typename T, typename U>
+ void response(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, T freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const;
+
+ /*!
+ * \brief Convenience method to compute the array factor of the station for
+ * list of (frequency, reference frequency) pairs.
+ *
+ * \param count Number of frequencies.
+ * \param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * \param freq Input iterator for a list of frequencies (Hz) of length
+ * \p count.
+ * \param direction Direction of arrival (ITRF, m).
+ * \param freq0 %Station beam former reference frequency (Hz).
+ * \param station0 %Station beam former reference direction (ITRF, m).
+ * \param tile0 Tile beam former reference direction (ITRF, m).
+ * \param buffer Output iterator with room for \p count instances of type
+ * ::diag22c_t.
+ *
+ * \see arrayFactor(real_t time, real_t freq, const vector3r_t &direction,
+ * real_t freq0, const vector3r_t &station0, const vector3r_t &tile0) const
+ */
+ template <typename T, typename U>
+ void arrayFactor(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, T freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const;
+
+ // @}
+
+private:
+ raw_array_factor_t fieldArrayFactor(const AntennaField::ConstPtr &field,
+ real_t time, real_t freq, const vector3r_t &direction, real_t freq0,
+ const vector3r_t &position0, const vector3r_t &direction0) const;
+
+private:
+ string itsName;
+ vector3r_t itsPosition;
+ vector3r_t itsPhaseReference;
+ FieldList itsFields;
+};
+
+// @}
+
+//# ------------------------------------------------------------------------- //
+//# - Implementation: Station - //
+//# ------------------------------------------------------------------------- //
+
+template <typename T, typename U>
+void Station::response(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const
+{
+ for(unsigned int i = 0; i < count; ++i)
+ {
+ *buffer++ = response(time, *freq++, direction, freq0, station0, tile0);
+ }
+}
+
+template <typename T, typename U>
+void Station::arrayFactor(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const
+{
+ for(unsigned int i = 0; i < count; ++i)
+ {
+ *buffer++ = arrayFactor(time, *freq++, direction, freq0, station0,
+ tile0);
+ }
+}
+
+template <typename T, typename U>
+void Station::response(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, T freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const
+{
+ for(unsigned int i = 0; i < count; ++i)
+ {
+ *buffer++ = response(time, *freq++, direction, *freq0++, station0,
+ tile0);
+ }
+}
+
+template <typename T, typename U>
+void Station::arrayFactor(unsigned int count, real_t time, T freq,
+ const vector3r_t &direction, T freq0, const vector3r_t &station0,
+ const vector3r_t &tile0, U buffer) const
+{
+ for(unsigned int i = 0; i < count; ++i)
+ {
+ *buffer++ = arrayFactor(time, *freq++, direction, *freq0++, station0,
+ tile0);
+ }
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/TileAntenna.h b/StationResponse/include/StationResponse/TileAntenna.h
new file mode 100644
index 0000000000000000000000000000000000000000..3695aa3e81657a3bcdbe5bfb91726287d5cf9040
--- /dev/null
+++ b/StationResponse/include/StationResponse/TileAntenna.h
@@ -0,0 +1,68 @@
+//# TileAntenna.h: Semi-analytical model of a LOFAR HBA tile.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_TILEANTENNA_H
+#define LOFAR_STATIONRESPONSE_TILEANTENNA_H
+
+// \file
+// Semi-analytical model of a LOFAR HBA tile.
+
+#include <StationResponse/AntennaModelHBA.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+class TileAntenna: public AntennaModelHBA
+{
+public:
+ typedef shared_ptr<TileAntenna> Ptr;
+ typedef shared_ptr<const TileAntenna> ConstPtr;
+
+ typedef static_array<vector3r_t, 16> TileConfig;
+
+ explicit TileAntenna(const TileConfig &config);
+
+ void setConfig(const TileConfig &config);
+
+ const TileConfig &config() const;
+
+ virtual raw_array_factor_t rawArrayFactor(real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const;
+
+ virtual matrix22c_t elementResponse(real_t freq,
+ const vector3r_t &direction) const;
+
+private:
+ TileConfig itsConfig;
+};
+
+// @}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/include/StationResponse/Types.h b/StationResponse/include/StationResponse/Types.h
new file mode 100644
index 0000000000000000000000000000000000000000..8565e127528dc692ae7750e4b9c4ee09c7cad037
--- /dev/null
+++ b/StationResponse/include/StationResponse/Types.h
@@ -0,0 +1,234 @@
+//# Types.h: Types used in this library.
+//#
+//# Copyright (C) 2013
+//# 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_STATIONRESPONSE_TYPES_H
+#define LOFAR_STATIONRESPONSE_TYPES_H
+
+// \file
+// Types used in this library.
+
+#include <Common/lofar_complex.h>
+#include <Common/StreamUtil.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+// \addtogroup StationResponse
+// @{
+
+/**
+ * \brief Array of static size.
+ *
+ * This struct wraps a (unidimensional) C array. Like C arrays, it is possible
+ * to use initializers, for example:
+ * \code
+ * int foo[3] = {1, 2, 3};
+ * static_array<int, 3> bar = {{1, 2, 3}};
+ * \endcode
+ * Note the \e double curly braces.
+ *
+ * Unlike C arrays, it is possible to return instances of this struct.
+ * \code
+ * static_array<int, 3> foo()
+ * {
+ * static_array<int, 3> bar = {{1, 2, 3}};
+ * return bar;
+ * }
+ * \endcode
+ * \see boost::array
+ */
+template <typename T, size_t N>
+struct static_array
+{
+ typedef T *iterator;
+ typedef const T *const_iterator;
+
+ T __data[N];
+
+ /** Returns the size of the array. */
+ static size_t size();
+
+ /**
+ * \brief Read-only access to a specific array element.
+ * \param n Index (0-based) of the element to access.
+ * \return A constant reference to the element at index \p n.
+ *
+ * This operator provides array-style element access. No range check is
+ * performed on the index \p n. The result of this operator is undefined
+ * for out of range indices.
+ */
+ const T &operator[](size_t n) const;
+
+ /**
+ * \brief Access to a specific array element.
+ * \param n Index (0-based) of the element to access.
+ * \return A non-constant reference to the element at index \p n.
+ *
+ * This operator provides array-style element access. No range check is
+ * performed on the index \p n. The result of this operator is undefined
+ * for out of range indices.
+ */
+ T &operator[](size_t n);
+
+ /** Returns an iterator that points to the first element in the array. */
+ iterator begin();
+
+ /** Returns an iterator that points one past the last element in the
+ * array.
+ */
+ iterator end();
+
+ /** Returns a read-only iterator that points to the first element in the
+ * array.
+ */
+ const_iterator begin() const;
+
+ /** Returns a read-only iterator that points one past the last element in
+ * the array.
+ */
+ const_iterator end() const;
+};
+
+/** Print the contents of a static array. */
+template <typename T, size_t N>
+ostream &operator<<(ostream &out, const static_array<T,N> &obj);
+
+/** Type used for real scalars. */
+typedef double real_t;
+
+/** Type used for complex scalars. */
+typedef dcomplex complex_t;
+
+/** Type used for 2-dimensional real vectors. */
+typedef static_array<real_t, 2> vector2r_t;
+
+/** Type used for 3-dimensional real vectors. */
+typedef static_array<real_t, 3> vector3r_t;
+
+/** Type used for 2x2 real diagonal matrices. */
+typedef static_array<real_t, 2> diag22r_t;
+
+/** Type used for 2x2 complex diagonal matrices. */
+typedef static_array<complex_t, 2> diag22c_t;
+
+/** Type used for 2x2 real matrices. */
+typedef static_array<static_array<real_t, 2>, 2> matrix22r_t;
+
+/** Type used for 2x2 complex matrices. */
+typedef static_array<static_array<complex_t, 2>, 2> matrix22c_t;
+
+/** Response of an array of antenna elements. */
+struct raw_response_t
+{
+ /** Combined response of all (enabled) antenna elements in the array. */
+ matrix22c_t response;
+
+ /** Number of antenna elements contributing to the combined response, per
+ * polarization.
+ */
+ diag22r_t weight;
+};
+
+/** Array factor of an array of antenna elements. A wave front of an incident
+ * plane wave will arrive at each antenna element at a potentially different
+ * time. The time of arrival depends on the location of the antenna element and
+ * the direction of arrival of the plane wave. With respect to a pre-defined
+ * phase reference location, there is a (possibly negative) delay between the
+ * arrival of a wave front at a given antenna element and the arrival of the
+ * same wave front at the phase reference location. The array factor is the sum
+ * of the phase shifts due to these delays. It describes the "sensitivity" of
+ * the array as a function of direction.
+ */
+struct raw_array_factor_t
+{
+ /** Array factor due to all (enabled) antenna elements in the array. */
+ diag22c_t factor;
+
+ /** Number of antenna elements contributing to the array factor, per
+ * polarization.
+ */
+ diag22r_t weight;
+};
+
+// @}
+
+//# ------------------------------------------------------------------------- //
+//# - Implementation: static_array - //
+//# ------------------------------------------------------------------------- //
+
+template <typename T, size_t N>
+inline const T &static_array<T, N>::operator[](size_t n) const
+{
+ return __data[n];
+}
+
+template <typename T, size_t N>
+inline T &static_array<T, N>::operator[](size_t n)
+{
+ return __data[n];
+}
+
+template <typename T, size_t N>
+inline typename static_array<T, N>::iterator static_array<T, N>::begin()
+{
+ return __data;
+}
+
+template <typename T, size_t N>
+inline typename static_array<T, N>::iterator static_array<T, N>::end()
+{
+ return __data + N;
+}
+
+template <typename T, size_t N>
+inline typename static_array<T, N>::const_iterator static_array<T, N>::begin()
+ const
+{
+ return __data;
+}
+
+template <typename T, size_t N>
+inline typename static_array<T, N>::const_iterator static_array<T, N>::end()
+ const
+{
+ return __data + N;
+}
+
+template <typename T, size_t N>
+inline size_t static_array<T, N>::size()
+{
+ return N;
+}
+
+template <typename T, size_t N>
+ostream &operator<<(ostream &out, const static_array<T,N> &obj)
+{
+ print(out, obj.begin(), obj.end());
+ return out;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
+
+#endif
diff --git a/StationResponse/src/AntennaField.cc b/StationResponse/src/AntennaField.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d3ca5879353d523a70109c631d1046055f08b713
--- /dev/null
+++ b/StationResponse/src/AntennaField.cc
@@ -0,0 +1,233 @@
+//# AntennaField.cc: Representation of a LOFAR antenna field, with methods to
+//# compute its response to incoming radiation.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/AntennaField.h>
+#include <StationResponse/Constants.h>
+#include <StationResponse/MathUtil.h>
+#include <ElementResponse/ElementResponse.h>
+#include <casacore/measures/Measures/MeasFrame.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+AntennaField::AntennaField(const string &name,
+ const CoordinateSystem &coordinates)
+ : itsName(name),
+ itsCoordinateSystem(coordinates),
+ itsNCPCacheTime(-1)
+{
+ vector3r_t ncp = {{0.0, 0.0, 1.0}};
+ itsNCP.reset(new ITRFDirection(position(), ncp));
+}
+
+AntennaField::~AntennaField()
+{
+}
+
+const string &AntennaField::name() const
+{
+ return itsName;
+}
+
+const vector3r_t &AntennaField::position() const
+{
+ return itsCoordinateSystem.origin;
+}
+
+const AntennaField::CoordinateSystem &AntennaField::coordinates() const
+{
+ return itsCoordinateSystem;
+}
+
+void AntennaField::addAntenna(const Antenna &antenna)
+{
+ itsAntennae.push_back(antenna);
+}
+
+size_t AntennaField::nAntennae() const
+{
+ return itsAntennae.size();
+}
+
+const AntennaField::Antenna &AntennaField::antenna(size_t n) const
+{
+ return itsAntennae[n];
+}
+
+AntennaField::Antenna &AntennaField::antenna(size_t n)
+{
+ return itsAntennae[n];
+}
+
+AntennaField::AntennaList::const_iterator AntennaField::beginAntennae() const
+{
+ return itsAntennae.begin();
+}
+
+AntennaField::AntennaList::const_iterator AntennaField::endAntennae() const
+{
+ return itsAntennae.end();
+}
+
+vector3r_t AntennaField::ncp(real_t time) const
+{
+ if(time != itsNCPCacheTime)
+ {
+ itsNCPCacheDirection = itsNCP->at(time);
+ itsNCPCacheTime = time;
+ }
+
+ return itsNCPCacheDirection;
+}
+
+vector3r_t AntennaField::itrf2field(const vector3r_t &itrf) const
+{
+ const CoordinateSystem::Axes &axes = itsCoordinateSystem.axes;
+ vector3r_t station = {{dot(axes.p, itrf), dot(axes.q, itrf),
+ dot(axes.r, itrf)}};
+ return station;
+}
+
+matrix22r_t AntennaField::rotation(real_t time, const vector3r_t &direction)
+ const
+{
+ // Compute the cross product of the NCP and the target direction. This
+ // yields a vector tangent to the celestial sphere at the target
+ // direction, pointing towards the East (the direction of +Y in the IAU
+ // definition, or positive right ascension).
+ // Test if the direction is equal to the NCP. If it is, take a random
+ // vector orthogonal to v1 (the east is not defined here).
+ vector3r_t v1;
+ if (abs(ncp(time)[0]-direction[0])<1e-9 &&
+ abs(ncp(time)[1]-direction[1])<1e-9 &&
+ abs(ncp(time)[2]-direction[2])<1e-9) {
+ // Make sure v1 is orthogonal to ncp(time). The first two components
+ // of v1 are arbitrary.
+ v1[0]=1.;
+ v1[1]=0.;
+ v1[2]=-(ncp(time)[0]*v1[0]+ncp(time)[1]*v1[1])/ncp(time)[2];
+ v1=normalize(v1);
+ } else {
+ v1 = normalize(cross(ncp(time), direction));
+ }
+
+ // Compute the cross product of the antenna field normal (R) and the
+ // target direction. This yields a vector tangent to the topocentric
+ // spherical coordinate system at the target direction, pointing towards
+ // the direction of positive phi (which runs East over North around the
+ // pseudo zenith).
+ // Test if the normal is equal to the target direction. If it is, take
+ // a random vector orthogonal to the normal.
+ vector3r_t v2;
+ if (abs(itsCoordinateSystem.axes.r[0]-direction[0])<1e-9 &&
+ abs(itsCoordinateSystem.axes.r[1]-direction[1])<1e-9 &&
+ abs(itsCoordinateSystem.axes.r[2]-direction[2])<1e-9) {
+ // Make sure v2 is orthogonal to r. The first two components
+ // of v2 are arbitrary.
+ v2[0]=1.;
+ v2[1]=0.;
+ v2[2]=-(itsCoordinateSystem.axes.r[0]*v2[0]+
+ itsCoordinateSystem.axes.r[1]*v2[1])/
+ itsCoordinateSystem.axes.r[2];
+ v2=normalize(v2);
+ } else {
+ v2 = normalize(cross(itsCoordinateSystem.axes.r, direction));
+ }
+
+ // Compute the cosine and sine of the parallactic angle, i.e. the angle
+ // between v1 and v2, both tangent to a latitude circle of their
+ // respective spherical coordinate systems.
+ real_t coschi = dot(v1, v2);
+ real_t sinchi = dot(cross(v1, v2), direction);
+
+ // The input coordinate system is a right handed system with its third
+ // axis along the direction of propagation (IAU +Z). The output
+ // coordinate system is right handed as well, but its third axis points
+ // in the direction of arrival (i.e. exactly opposite).
+ //
+ // Because the electromagnetic field is always perpendicular to the
+ // direction of propagation, we only need to relate the (X, Y) axes of
+ // the input system to the corresponding (theta, phi) axes of the output
+ // system.
+ //
+ // To this end, we first rotate the input system around its third axis
+ // to align the Y axis with the phi axis. The X and theta axis are
+ // parallel after this rotation, but point in opposite directions. To
+ // align the X axis with the theta axis, we flip it.
+ //
+ // The Jones matrix to align the Y axis with the phi axis when these are
+ // separated by an angle phi (measured counter-clockwise around the
+ // direction of propagation, looking towards the origin), is given by:
+ //
+ // [ cos(phi) sin(phi)]
+ // [-sin(phi) cos(phi)]
+ //
+ // Here, cos(phi) and sin(phi) can be computed directly, without having
+ // to compute phi first (see the computation of coschi and sinchi
+ // above).
+ //
+ // Now, sinchi as computed above is opposite to sin(phi), because the
+ // direction used in the computation is the direction of arrival instead
+ // of the direction of propagation. Therefore, the sign of sinchi needs
+ // to be reversed. Furthermore, as explained above, the X axis has to be
+ // flipped to align with the theta axis. The Jones matrix returned from
+ // this function is therefore given by:
+ //
+ // [-coschi sinchi]
+ // [ sinchi coschi]
+ matrix22r_t rotation = {{{{-coschi, sinchi}}, {{sinchi, coschi}}}};
+ return rotation;
+}
+
+matrix22c_t AntennaField::response(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ return normalize(rawResponse(time, freq, direction, direction0));
+}
+
+diag22c_t AntennaField::arrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ return normalize(rawArrayFactor(time, freq, direction, direction0));
+}
+
+raw_response_t AntennaField::rawResponse(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ raw_array_factor_t af = rawArrayFactor(time, freq, direction, direction0);
+
+ raw_response_t result;
+ result.response = elementResponse(time, freq, direction);
+ result.response[0][0] *= af.factor[0];
+ result.response[0][1] *= af.factor[0];
+ result.response[1][0] *= af.factor[1];
+ result.response[1][1] *= af.factor[1];
+ result.weight = af.weight;
+ return result;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/AntennaFieldHBA.cc b/StationResponse/src/AntennaFieldHBA.cc
new file mode 100644
index 0000000000000000000000000000000000000000..ff099fb04a22ea2987f0f8c397f353bb69b66e6a
--- /dev/null
+++ b/StationResponse/src/AntennaFieldHBA.cc
@@ -0,0 +1,77 @@
+//# AntennaFieldHBA.cc: Representation of an HBA antenna field.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/AntennaFieldHBA.h>
+#include <StationResponse/MathUtil.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+AntennaFieldHBA::AntennaFieldHBA(const string &name,
+ const CoordinateSystem &coordinates, const AntennaModelHBA::ConstPtr &model)
+ : AntennaField(name, coordinates),
+ itsAntennaModel(model)
+{
+}
+
+matrix22c_t AntennaFieldHBA::response(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ return itsAntennaModel->response(freq, itrf2field(direction),
+ itrf2field(direction0)) * rotation(time, direction);
+}
+
+diag22c_t AntennaFieldHBA::arrayFactor(real_t, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ return itsAntennaModel->arrayFactor(freq, itrf2field(direction),
+ itrf2field(direction0));
+}
+
+raw_response_t AntennaFieldHBA::rawResponse(real_t time, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ raw_response_t result = itsAntennaModel->rawResponse(freq,
+ itrf2field(direction), itrf2field(direction0));
+ result.response = result.response * rotation(time, direction);
+ return result;
+}
+
+raw_array_factor_t AntennaFieldHBA::rawArrayFactor(real_t, real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ return itsAntennaModel->rawArrayFactor(freq, itrf2field(direction),
+ itrf2field(direction0));
+}
+
+matrix22c_t AntennaFieldHBA::elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction) const
+{
+ return itsAntennaModel->elementResponse(freq, itrf2field(direction))
+ * rotation(time, direction);
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/AntennaFieldLBA.cc b/StationResponse/src/AntennaFieldLBA.cc
new file mode 100644
index 0000000000000000000000000000000000000000..90838c135e6790e2a9d757bf6a108a24d3933059
--- /dev/null
+++ b/StationResponse/src/AntennaFieldLBA.cc
@@ -0,0 +1,54 @@
+//# AntennaFieldLBA.cc: Representation of an LBA antenna field.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/AntennaFieldLBA.h>
+#include <StationResponse/MathUtil.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+AntennaFieldLBA::AntennaFieldLBA(const string &name,
+ const CoordinateSystem &coordinates, const AntennaModelLBA::ConstPtr &model)
+ : AntennaField(name, coordinates),
+ itsAntennaModel(model)
+{
+}
+
+raw_array_factor_t AntennaFieldLBA::rawArrayFactor(real_t, real_t,
+ const vector3r_t&, const vector3r_t&) const
+{
+ raw_array_factor_t af = {{{1.0, 1.0}}, {{1.0, 1.0}}};
+ return af;
+}
+
+matrix22c_t AntennaFieldLBA::elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction) const
+{
+ return itsAntennaModel->response(freq, itrf2field(direction))
+ * rotation(time, direction);
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/AntennaModelHBA.cc b/StationResponse/src/AntennaModelHBA.cc
new file mode 100644
index 0000000000000000000000000000000000000000..52ca0accc5428c38aa73964bd74046d36dcda223
--- /dev/null
+++ b/StationResponse/src/AntennaModelHBA.cc
@@ -0,0 +1,64 @@
+//# AntennaModelHBA.cc: HBA antenna model interface definitions.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/AntennaModelHBA.h>
+#include <StationResponse/MathUtil.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+AntennaModelHBA::~AntennaModelHBA()
+{
+}
+
+matrix22c_t AntennaModelHBA::response(real_t freq, const vector3r_t &direction,
+ const vector3r_t &direction0) const
+{
+ return normalize(rawResponse(freq, direction, direction0));
+}
+
+diag22c_t AntennaModelHBA::arrayFactor(real_t freq, const vector3r_t &direction,
+ const vector3r_t &direction0) const
+{
+ return normalize(rawArrayFactor(freq, direction, direction0));
+}
+
+raw_response_t AntennaModelHBA::rawResponse(real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ raw_array_factor_t af = rawArrayFactor(freq, direction, direction0);
+
+ raw_response_t result;
+ result.response = elementResponse(freq, direction);
+ result.response[0][0] *= af.factor[0];
+ result.response[0][1] *= af.factor[0];
+ result.response[1][0] *= af.factor[1];
+ result.response[1][1] *= af.factor[1];
+ result.weight = af.weight;
+ return result;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/AntennaModelLBA.cc b/StationResponse/src/AntennaModelLBA.cc
new file mode 100644
index 0000000000000000000000000000000000000000..87d02b945013406efe05d18c4481dc98a2651d42
--- /dev/null
+++ b/StationResponse/src/AntennaModelLBA.cc
@@ -0,0 +1,36 @@
+//# AntennaModelLBA.cc: LBA antenna model interface definitions.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/AntennaModelLBA.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+AntennaModelLBA::~AntennaModelLBA()
+{
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/CMakeLists.txt b/StationResponse/src/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e31a7a54574b0e908c6ba2760cf49626696f1850
--- /dev/null
+++ b/StationResponse/src/CMakeLists.txt
@@ -0,0 +1,20 @@
+# $Id$
+
+include(LofarPackageVersion)
+
+lofar_add_library(stationresponse
+ Package__Version.cc
+ AntennaField.cc
+ AntennaFieldHBA.cc
+ AntennaFieldLBA.cc
+ AntennaModelHBA.cc
+ AntennaModelLBA.cc
+ DualDipoleAntenna.cc
+ ITRFDirection.cc
+ LofarMetaDataUtil.cc
+ MathUtil.cc
+ Station.cc
+ TileAntenna.cc
+ Types.cc)
+
+lofar_add_bin_program(makeresponseimage makeresponseimage.cc)
diff --git a/StationResponse/src/DualDipoleAntenna.cc b/StationResponse/src/DualDipoleAntenna.cc
new file mode 100644
index 0000000000000000000000000000000000000000..659d9f45ddb16ad04a66678d36fad174439fa152
--- /dev/null
+++ b/StationResponse/src/DualDipoleAntenna.cc
@@ -0,0 +1,51 @@
+//# DualDipoleAntenna.cc: Semi-analytical model of a LOFAR LBA dual dipole
+//# antenna.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/DualDipoleAntenna.h>
+#include <StationResponse/Constants.h>
+#include <StationResponse/MathUtil.h>
+#include <ElementResponse/ElementResponse.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+matrix22c_t DualDipoleAntenna::response(real_t freq,
+ const vector3r_t &direction) const
+{
+ // The positive X dipole direction is SW of the reference orientation,
+ // which translates to a phi coordinate of 5/4*pi in the topocentric
+ // spherical coordinate system. The phi coordinate is corrected for this
+ // offset before evaluating the antenna model.
+ vector2r_t thetaphi = cart2thetaphi(direction);
+ thetaphi[1] -= 5.0 * Constants::pi_4;
+ matrix22c_t response;
+ element_response_lba(freq, thetaphi[0], thetaphi[1],
+ reinterpret_cast<std::complex<double> (&)[2][2]>(response));
+ return response;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/ITRFDirection.cc b/StationResponse/src/ITRFDirection.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f7f9c49e02551bfac1b6afbdeb0c14acbeaa449d
--- /dev/null
+++ b/StationResponse/src/ITRFDirection.cc
@@ -0,0 +1,77 @@
+//# ITRFDirection.cc: Functor that maps time to an ITRF direction.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/ITRFDirection.h>
+
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MEpoch.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+ITRFDirection::ITRFDirection(const vector3r_t &position,
+ const vector2r_t &direction)
+{
+ casacore::MVPosition mvPosition(position[0], position[1], position[2]);
+ casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
+ itsFrame = casacore::MeasFrame(casacore::MEpoch(), mPosition);
+
+ // Order of angles seems to be longitude (along the equator), lattitude
+ // (towards the pole).
+ casacore::MVDirection mvDirection(direction[0], direction[1]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+ itsConverter = casacore::MDirection::Convert(mDirection,
+ casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
+}
+
+ITRFDirection::ITRFDirection(const vector3r_t &position,
+ const vector3r_t &direction)
+{
+ casacore::MVPosition mvPosition(position[0], position[1], position[2]);
+ casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
+ itsFrame = casacore::MeasFrame(casacore::MEpoch(), mPosition);
+
+ casacore::MVDirection mvDirection(direction[0], direction[1], direction[2]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+ itsConverter = casacore::MDirection::Convert(mDirection,
+ casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
+}
+
+vector3r_t ITRFDirection::at(real_t time) const
+{
+ // Cannot use MeasFrame::resetEpoch(Double), because that assumes the
+ // argument is UTC in (fractional) days (MJD).
+ itsFrame.resetEpoch(casacore::Quantity(time, "s"));
+
+ const casacore::MDirection &mITRF = itsConverter();
+ const casacore::MVDirection &mvITRF = mITRF.getValue();
+
+ vector3r_t itrf = {{mvITRF(0), mvITRF(1), mvITRF(2)}};
+ return itrf;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/LofarMetaDataUtil.cc b/StationResponse/src/LofarMetaDataUtil.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4d96af8dfe4f3f0c1b7b945d625e66bc5d977282
--- /dev/null
+++ b/StationResponse/src/LofarMetaDataUtil.cc
@@ -0,0 +1,335 @@
+//# LofarMetaDataUtil.cc: Utility functions to read the meta data relevant for
+//# simulating the beam from LOFAR observations stored in MS format.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/LofarMetaDataUtil.h>
+#include <StationResponse/AntennaFieldLBA.h>
+#include <StationResponse/AntennaFieldHBA.h>
+#include <StationResponse/MathUtil.h>
+#include <StationResponse/TileAntenna.h>
+#include <StationResponse/DualDipoleAntenna.h>
+
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/measures/Measures/MCDirection.h>
+#include <casacore/measures/Measures/MCPosition.h>
+#include <casacore/measures/Measures/MeasTable.h>
+#include <casacore/measures/Measures/MeasConvert.h>
+#include <casacore/measures/TableMeasures/ScalarMeasColumn.h>
+
+#include <stdexcept>
+
+#include <casacore/ms/MeasurementSets/MSAntenna.h>
+#include <casacore/ms/MSSel/MSSelection.h>
+#include <casacore/ms/MSSel/MSAntennaParse.h>
+#include <casacore/ms/MeasurementSets/MSAntennaColumns.h>
+#include <casacore/ms/MeasurementSets/MSDataDescription.h>
+#include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
+#include <casacore/ms/MeasurementSets/MSField.h>
+#include <casacore/ms/MeasurementSets/MSFieldColumns.h>
+#include <casacore/ms/MeasurementSets/MSObservation.h>
+#include <casacore/ms/MeasurementSets/MSObsColumns.h>
+#include <casacore/ms/MeasurementSets/MSPolarization.h>
+#include <casacore/ms/MeasurementSets/MSPolColumns.h>
+#include <casacore/ms/MeasurementSets/MSSpectralWindow.h>
+#include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+using namespace casacore;
+
+bool hasColumn(const Table &table, const string &column)
+{
+ return table.tableDesc().isColumn(column);
+}
+
+bool hasSubTable(const Table &table, const string &name)
+{
+ return table.keywordSet().isDefined(name);
+}
+
+Table getSubTable(const Table &table, const string &name)
+{
+ return table.keywordSet().asTable(name);
+}
+
+TileAntenna::TileConfig readTileConfig(const Table &table, unsigned int row)
+{
+ ROArrayQuantColumn<Double> c_tile_offset(table, "TILE_ELEMENT_OFFSET", "m");
+
+ // Read tile configuration for HBA antenna fields.
+ Matrix<Quantity> aips_offset = c_tile_offset(row);
+ assert(aips_offset.ncolumn() == TileAntenna::TileConfig::size());
+
+ TileAntenna::TileConfig config;
+ for(unsigned int i = 0; i < config.size(); ++i)
+ {
+ config[i][0] = aips_offset(0, i).getValue();
+ config[i][1] = aips_offset(1, i).getValue();
+ config[i][2] = aips_offset(2, i).getValue();
+ }
+ return config;
+}
+
+void transformToFieldCoordinates(TileAntenna::TileConfig &config,
+ const AntennaField::CoordinateSystem::Axes &axes)
+{
+ for(unsigned int i = 0; i < config.size(); ++i)
+ {
+ const vector3r_t position = config[i];
+ config[i][0] = dot(position, axes.p);
+ config[i][1] = dot(position, axes.q);
+ config[i][2] = dot(position, axes.r);
+ }
+}
+
+AntennaField::CoordinateSystem readCoordinateSystemAartfaac(
+ const Table &table, unsigned int id)
+{
+ ROArrayQuantColumn<Double> c_position(table, "POSITION", "m");
+
+ // Read antenna field center (ITRF).
+ Vector<Quantity> aips_position = c_position(id);
+ assert(aips_position.size() == 3);
+
+ vector3r_t position = {{aips_position(0).getValue(),
+ aips_position(1).getValue(), aips_position(2).getValue()}};
+
+ TableRecord keywordset = table.keywordSet();
+ Matrix<double> aips_axes;
+ keywordset.get("AARTFAAC_COORDINATE_AXES", aips_axes);
+ assert(aips_axes.shape().isEqual(IPosition(2, 3, 3)));
+
+ vector3r_t p = {{aips_axes(0, 0), aips_axes(1, 0), aips_axes(2, 0)}};
+ vector3r_t q = {{aips_axes(0, 1), aips_axes(1, 1), aips_axes(2, 1)}};
+ vector3r_t r = {{aips_axes(0, 2), aips_axes(1, 2), aips_axes(2, 2)}};
+
+ AntennaField::CoordinateSystem system = {position, {p, q, r}};
+
+ return system;
+}
+
+AntennaField::CoordinateSystem readCoordinateSystem(const Table &table,
+ unsigned int id)
+{
+ ROArrayQuantColumn<Double> c_position(table, "POSITION", "m");
+ ROArrayQuantColumn<Double> c_axes(table, "COORDINATE_AXES", "m");
+
+ // Read antenna field center (ITRF).
+ Vector<Quantity> aips_position = c_position(id);
+ assert(aips_position.size() == 3);
+
+ vector3r_t position = {{aips_position(0).getValue(),
+ aips_position(1).getValue(), aips_position(2).getValue()}};
+
+ // Read antenna field coordinate axes (ITRF).
+ Matrix<Quantity> aips_axes = c_axes(id);
+ assert(aips_axes.shape().isEqual(IPosition(2, 3, 3)));
+
+ vector3r_t p = {{aips_axes(0, 0).getValue(), aips_axes(1, 0).getValue(),
+ aips_axes(2, 0).getValue()}};
+ vector3r_t q = {{aips_axes(0, 1).getValue(), aips_axes(1, 1).getValue(),
+ aips_axes(2, 1).getValue()}};
+ vector3r_t r = {{aips_axes(0, 2).getValue(), aips_axes(1, 2).getValue(),
+ aips_axes(2, 2).getValue()}};
+
+ AntennaField::CoordinateSystem system = {position, {p, q, r}};
+ return system;
+}
+
+void readAntennae(const Table &table, unsigned int id,
+ const AntennaField::Ptr &field)
+{
+ ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
+ ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG");
+
+ // Read element offsets and flags.
+ Matrix<Quantity> aips_offset = c_offset(id);
+ assert(aips_offset.shape().isEqual(IPosition(2, 3, aips_offset.ncolumn())));
+
+ Matrix<Bool> aips_flag = c_flag(id);
+ assert(aips_flag.shape().isEqual(IPosition(2, 2, aips_offset.ncolumn())));
+
+ for(size_t i = 0; i < aips_offset.ncolumn(); ++i)
+ {
+ AntennaField::Antenna antenna;
+ antenna.position[0] = aips_offset(0, i).getValue();
+ antenna.position[1] = aips_offset(1, i).getValue();
+ antenna.position[2] = aips_offset(2, i).getValue();
+ antenna.enabled[0] = !aips_flag(0, i);
+ antenna.enabled[1] = !aips_flag(1, i);
+ field->addAntenna(antenna);
+ }
+}
+
+AntennaField::Ptr readAntennaField(const Table &table, unsigned int id)
+{
+ AntennaField::Ptr field;
+ AntennaField::CoordinateSystem system = readCoordinateSystem(table, id);
+
+ ROScalarColumn<String> c_name(table, "NAME");
+ const string &name = c_name(id);
+
+ if(name == "LBA")
+ {
+ DualDipoleAntenna::Ptr model(new DualDipoleAntenna());
+ field = AntennaField::Ptr(new AntennaFieldLBA(name, system, model));
+ }
+ else // HBA, HBA0, HBA1
+ {
+ TileAntenna::TileConfig config = readTileConfig(table, id);
+ transformToFieldCoordinates(config, system.axes);
+
+ TileAntenna::Ptr model(new TileAntenna(config));
+ field = AntennaField::Ptr(new AntennaFieldHBA(name, system, model));
+ }
+
+ readAntennae(table, id, field);
+ return field;
+}
+
+AntennaField::Ptr readAntennaFieldAartfaac(const Table &table, const string &ant_type,
+ unsigned int id)
+{
+ AntennaField::Ptr field;
+ AntennaField::CoordinateSystem system = readCoordinateSystemAartfaac(table, id);
+
+ if (ant_type == "LBA")
+ {
+ DualDipoleAntenna::Ptr model(new DualDipoleAntenna());
+ field = AntennaField::Ptr(new AntennaFieldLBA(ant_type, system, model));
+ }
+ else // HBA
+ {
+ // TODO: implement this
+ throw std::runtime_error("HBA for Aartfaac is not implemented yet.");
+ }
+
+ // Add only one antenna to the field (no offset, always enabled)
+ AntennaField::Antenna antenna;
+ antenna.position[0] = 0.;
+ antenna.position[1] = 0.;
+ antenna.position[2] = 0.;
+ antenna.enabled[0] = true;
+ antenna.enabled[1] = true;
+
+ field->addAntenna(antenna);
+
+ return field;
+}
+
+void readStationPhaseReference(const Table &table, unsigned int id,
+ const Station::Ptr &station)
+{
+ const string columnName("LOFAR_PHASE_REFERENCE");
+ if(hasColumn(table, columnName))
+ {
+ ROScalarMeasColumn<MPosition> c_reference(table, columnName);
+ MPosition mReference = MPosition::Convert(c_reference(id),
+ MPosition::ITRF)();
+ MVPosition mvReference = mReference.getValue();
+ vector3r_t reference = {{mvReference(0), mvReference(1),
+ mvReference(2)}};
+
+ station->setPhaseReference(reference);
+ }
+}
+
+Station::Ptr readStation(const MeasurementSet &ms, unsigned int id)
+{
+ ROMSAntennaColumns antenna(ms.antenna());
+ assert(antenna.nrow() > id && !antenna.flagRow()(id));
+
+ // Get station name.
+ const string name(antenna.name()(id));
+
+ // Get station position (ITRF).
+ MPosition mPosition = MPosition::Convert(antenna.positionMeas()(id),
+ MPosition::ITRF)();
+ MVPosition mvPosition = mPosition.getValue();
+ const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
+
+ // Create station.
+ Station::Ptr station(new Station(name, position));
+
+ // Read phase reference position (if available).
+ readStationPhaseReference(ms.antenna(), id, station);
+
+ // Read antenna field information.
+ ROScalarColumn<String> telescope_name_col(getSubTable(ms, "OBSERVATION"),
+ "TELESCOPE_NAME");
+ string telescope_name = telescope_name_col(0);
+
+ if (telescope_name == "LOFAR")
+ {
+ Table tab_field = getSubTable(ms, "LOFAR_ANTENNA_FIELD");
+ tab_field = tab_field(tab_field.col("ANTENNA_ID") == static_cast<Int>(id));
+
+ for(size_t i = 0; i < tab_field.nrow(); ++i)
+ {
+ station->addField(readAntennaField(tab_field, i));
+ }
+ }
+ else if (telescope_name == "AARTFAAC")
+ {
+ ROScalarColumn<String> ant_type_col(getSubTable(ms, "OBSERVATION"),
+ "AARTFAAC_ANTENNA_TYPE");
+ string ant_type = ant_type_col(0);
+
+ Table tab_field = getSubTable(ms, "ANTENNA");
+ station -> addField(readAntennaFieldAartfaac(tab_field, ant_type, id));
+ }
+
+ return station;
+}
+
+MDirection readTileBeamDirection(const casacore::MeasurementSet &ms) {
+ MDirection tileBeamDir;
+
+ Table fieldTable = getSubTable(ms, "FIELD");
+
+ if (fieldTable.nrow() != 1) {
+ throw std::runtime_error("MS has multiple fields, this does not work with the LOFAR beam library.");
+ }
+
+ if (hasColumn(fieldTable, "LOFAR_TILE_BEAM_DIR"))
+ {
+ ROArrayMeasColumn<MDirection> tileBeamCol(fieldTable,
+ "LOFAR_TILE_BEAM_DIR");
+ tileBeamDir = *(tileBeamCol(0).data());
+ }
+ else
+ {
+ ROArrayMeasColumn<MDirection> tileBeamCol(fieldTable,
+ "DELAY_DIR");
+ tileBeamDir = *(tileBeamCol(0).data());
+ }
+
+ return tileBeamDir;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/MathUtil.cc b/StationResponse/src/MathUtil.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1ce4c5af2ee7bf1169d138edacbd2beb6657c881
--- /dev/null
+++ b/StationResponse/src/MathUtil.cc
@@ -0,0 +1,32 @@
+//# MathUtil.cc: Various mathematical operations on vectors and matrices.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/MathUtil.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/Station.cc b/StationResponse/src/Station.cc
new file mode 100644
index 0000000000000000000000000000000000000000..7c4e21045a43e83bd37767c818b8f980fcd93c5c
--- /dev/null
+++ b/StationResponse/src/Station.cc
@@ -0,0 +1,176 @@
+//# Station.cc: Representation of the station beam former.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/Station.h>
+#include <StationResponse/MathUtil.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+Station::Station(const string &name, const vector3r_t &position)
+ : itsName(name),
+ itsPosition(position),
+ itsPhaseReference(position)
+{
+}
+
+const string &Station::name() const
+{
+ return itsName;
+}
+
+const vector3r_t &Station::position() const
+{
+ return itsPosition;
+}
+
+void Station::setPhaseReference(const vector3r_t &reference)
+{
+ itsPhaseReference = reference;
+}
+
+const vector3r_t &Station::phaseReference() const
+{
+ return itsPhaseReference;
+}
+
+void Station::addField(const AntennaField::ConstPtr &field)
+{
+ itsFields.push_back(field);
+}
+
+size_t Station::nFields() const
+{
+ return itsFields.size();
+}
+
+AntennaField::ConstPtr Station::field(size_t i) const
+{
+ return (i < itsFields.size() ? itsFields[i] : AntennaField::ConstPtr());
+}
+
+Station::FieldList::const_iterator Station::beginFields() const
+{
+ return itsFields.begin();
+}
+
+Station::FieldList::const_iterator Station::endFields() const
+{
+ return itsFields.end();
+}
+
+matrix22c_t Station::response(real_t time, real_t freq,
+ const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
+ const vector3r_t &tile0) const
+{
+ raw_response_t result = {{{{{}}, {{}}}}, {{}}};
+ for(FieldList::const_iterator field_it = beginFields(),
+ field_end = endFields(); field_it != field_end; ++field_it)
+ {
+ raw_array_factor_t field = fieldArrayFactor(*field_it, time, freq,
+ direction, freq0, phaseReference(), station0);
+
+ raw_response_t antenna = (*field_it)->rawResponse(time, freq,
+ direction, tile0);
+
+ result.response[0][0] += field.factor[0] * antenna.response[0][0];
+ result.response[0][1] += field.factor[0] * antenna.response[0][1];
+ result.response[1][0] += field.factor[1] * antenna.response[1][0];
+ result.response[1][1] += field.factor[1] * antenna.response[1][1];
+
+ result.weight[0] += field.weight[0] * antenna.weight[0];
+ result.weight[1] += field.weight[1] * antenna.weight[1];
+ }
+
+ return normalize(result);
+}
+
+diag22c_t Station::arrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
+ const vector3r_t &tile0) const
+{
+ raw_array_factor_t af = {{{}}, {{}}};
+ for(FieldList::const_iterator field_it = beginFields(),
+ field_end = endFields(); field_it != field_end; ++field_it)
+ {
+ raw_array_factor_t field = fieldArrayFactor(*field_it, time, freq,
+ direction, freq0, phaseReference(), station0);
+
+ raw_array_factor_t antenna = (*field_it)->rawArrayFactor(time, freq,
+ direction, tile0);
+
+ af.factor[0] += field.factor[0] * antenna.factor[0];
+ af.factor[1] += field.factor[1] * antenna.factor[1];
+ af.weight[0] += field.weight[0] * antenna.weight[0];
+ af.weight[1] += field.weight[1] * antenna.weight[1];
+ }
+
+ return normalize(af);
+}
+
+raw_array_factor_t
+Station::fieldArrayFactor(const AntennaField::ConstPtr &field,
+ real_t, real_t freq, const vector3r_t &direction, real_t freq0,
+ const vector3r_t &position0, const vector3r_t &direction0) const
+{
+ real_t k = Constants::_2pi * freq / Constants::c;
+ real_t k0 = Constants::_2pi * freq0 / Constants::c;
+
+ vector3r_t offset = field->position() - position0;
+
+ raw_array_factor_t af = {{{}}, {{}}};
+ typedef AntennaField::AntennaList AntennaList;
+ for(AntennaList::const_iterator antenna_it = field->beginAntennae(),
+ antenna_end = field->endAntennae(); antenna_it != antenna_end;
+ ++antenna_it)
+ {
+ if(!antenna_it->enabled[0] && !antenna_it->enabled[1])
+ {
+ continue;
+ }
+
+ vector3r_t position = offset + antenna_it->position;
+ real_t phase = k * dot(position, direction) - k0 * dot(position,
+ direction0);
+ complex_t shift = complex_t(cos(phase), sin(phase));
+
+ if(antenna_it->enabled[0])
+ {
+ af.factor[0] += shift;
+ ++af.weight[0];
+ }
+
+ if(antenna_it->enabled[1])
+ {
+ af.factor[1] += shift;
+ ++af.weight[1];
+ }
+ }
+
+ return af;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/TileAntenna.cc b/StationResponse/src/TileAntenna.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2813563ec96cad64914c102168e2470aef1854b6
--- /dev/null
+++ b/StationResponse/src/TileAntenna.cc
@@ -0,0 +1,100 @@
+//# TileAntenna.cc: Semi-analytical model of a LOFAR HBA tile.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/TileAntenna.h>
+#include <StationResponse/Constants.h>
+#include <StationResponse/MathUtil.h>
+#include <ElementResponse/ElementResponse.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+TileAntenna::TileAntenna(const TileConfig &config)
+ : itsConfig(config)
+{
+}
+
+void TileAntenna::setConfig(const TileConfig &config)
+{
+ itsConfig = config;
+}
+
+const TileAntenna::TileConfig &TileAntenna::config() const
+{
+ return itsConfig;
+}
+
+raw_array_factor_t TileAntenna::rawArrayFactor(real_t freq,
+ const vector3r_t &direction, const vector3r_t &direction0) const
+{
+ // Angular wave number.
+ real_t k = Constants::_2pi * freq / Constants::c;
+
+ // We need to compute the phase difference between a signal arriving from
+ // the target direction and a signal arriving from the reference direction,
+ // both relative to the center of the tile.
+ //
+ // This phase difference can be computed using a single dot product per
+ // dipole element by exploiting the fact that the dot product is
+ // distributive over vector addition:
+ //
+ // a . b + a . c = a . (b + c)
+ //
+ vector3r_t difference = direction - direction0;
+
+ complex_t af(0.0, 0.0);
+ for(TileConfig::const_iterator element_it = itsConfig.begin(),
+ element_end = itsConfig.end(); element_it != element_end; ++element_it)
+ {
+ // Compute the effective delay for a plane wave approaching from the
+ // direction of interest with respect to the position of element i
+ // when beam forming in the reference direction using time delays.
+ real_t shift = k * dot(difference, *element_it);
+ af += complex_t(cos(shift), sin(shift));
+ }
+
+ real_t size = itsConfig.size();
+ raw_array_factor_t result = {{{af, af}}, {{size, size}}};
+ return result;
+}
+
+matrix22c_t TileAntenna::elementResponse(real_t freq,
+ const vector3r_t &direction) const
+{
+ // The positive X dipole direction is SW of the reference orientation,
+ // which translates to a phi coordinate of 5/4*pi in the topocentric
+ // spherical coordinate system. The phi coordinate is corrected for this
+ // offset before evaluating the antenna model.
+ vector2r_t thetaphi = cart2thetaphi(direction);
+ thetaphi[1] -= 5.0 * Constants::pi_4;
+
+ matrix22c_t response;
+ element_response_hba(freq, thetaphi[0], thetaphi[1],
+ reinterpret_cast<std::complex<double> (&)[2][2]>(response));
+ return response;
+}
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/Types.cc b/StationResponse/src/Types.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5267ae75691117365e166f513d915a879f48ee8a
--- /dev/null
+++ b/StationResponse/src/Types.cc
@@ -0,0 +1,32 @@
+//# Types.cc: Declaration of types used in this library.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+#include <StationResponse/Types.h>
+
+namespace LOFAR
+{
+namespace StationResponse
+{
+
+} //# namespace StationResponse
+} //# namespace LOFAR
diff --git a/StationResponse/src/makeresponseimage.cc b/StationResponse/src/makeresponseimage.cc
new file mode 100644
index 0000000000000000000000000000000000000000..aa545ed0ca194e826fd40a80e320758ca6855a62
--- /dev/null
+++ b/StationResponse/src/makeresponseimage.cc
@@ -0,0 +1,636 @@
+//# makeresponseimage.cc: Generate images of the station response for a given
+//# MeasurementSet.
+//#
+//# Copyright (C) 2013
+//# 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 <lofar_config.h>
+
+#include <StationResponse/Package__Version.h>
+#include <StationResponse/LofarMetaDataUtil.h>
+#include <Common/InputParSet.h>
+#include <Common/lofar_sstream.h>
+#include <Common/LofarLogger.h>
+#include <Common/SystemUtil.h>
+#include <Common/Version.h>
+#include <casacore/coordinates/Coordinates/CoordinateSystem.h>
+#include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
+#include <casacore/coordinates/Coordinates/StokesCoordinate.h>
+#include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
+#include <casacore/images/Images/PagedImage.h>
+#include <casacore/measures/Measures/MCDirection.h>
+#include <casacore/measures/Measures/MCPosition.h>
+#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/measures/Measures/MeasConvert.h>
+#include <casacore/measures/Measures/MeasTable.h>
+#include <casacore/measures/Measures/MEpoch.h>
+#include <casacore/measures/Measures/MPosition.h>
+#include <casacore/ms/MeasurementSets/MeasurementSet.h>
+#include <casacore/ms/MeasurementSets/MSDataDescription.h>
+#include <casacore/ms/MeasurementSets/MSDataDescColumns.h>
+#include <casacore/ms/MeasurementSets/MSField.h>
+#include <casacore/ms/MeasurementSets/MSFieldColumns.h>
+#include <casacore/ms/MeasurementSets/MSObservation.h>
+#include <casacore/ms/MeasurementSets/MSObsColumns.h>
+#include <casacore/ms/MeasurementSets/MSSpectralWindow.h>
+#include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
+#include <casacore/tables/Tables/ExprNode.h>
+
+// There is no wrapped include file lofar_iterator.h.
+#include <iterator>
+
+using namespace casacore;
+using namespace LOFAR;
+using namespace LOFAR::StationResponse;
+using LOFAR::operator<<;
+
+namespace
+{
+ /*!
+ * \brief Map of ITRF directions required to compute an image of the
+ * station beam.
+ *
+ * The station beam library uses the ITRF coordinate system to express
+ * station positions and source directions. Since the Earth moves with
+ * respect to the sky, the ITRF coordinates of a source vary with time.
+ * This structure stores the ITRF coordinates for the station and tile beam
+ * former reference directions, as well as for a grid of points on the sky,
+ * along with the time for which these ITRF coordinates are valid.
+ */
+ struct ITRFDirectionMap
+ {
+ /*!
+ * \brief The time for which this ITRF direction map is valid (MJD(UTC)
+ * in seconds).
+ */
+ double_t time0;
+
+ /*!
+ * \brief Station beam former reference direction expressed in ITRF
+ * coordinates.
+ */
+ vector3r_t station0;
+
+ /*!
+ * \brief Tile beam former reference direction expressed in ITRF
+ * coordinates.
+ */
+ vector3r_t tile0;
+
+ /*!
+ * \brief ITRF coordinates for a grid of points on the sky.
+ */
+ Matrix<vector3r_t> directions;
+ };
+
+ /*!
+ * \brief Create an ITRFDirectionMap.
+ *
+ * \param coordinates Sky coordinate system definition.
+ * \param shape Number of points along the RA and DEC axis.
+ * \param epoch Time for which to compute the ITRF coordinates.
+ * \param position0 Station beam former reference position (phase
+ * reference).
+ * \param station0 Station beam former reference direction (pointing).
+ * \param tile0 Tile beam former reference direction (pointing).
+ */
+ ITRFDirectionMap makeDirectionMap(const DirectionCoordinate &coordinates,
+ const IPosition &shape, const MEpoch &epoch, const MPosition &position0,
+ const MDirection &station0, const MDirection &tile0);
+
+ /*!
+ * \brief Create a DirectionCoordinate instance that defines an image
+ * coordinate system on the sky (J2000).
+ *
+ * \param reference Direction that corresponds to the origin of the
+ * coordinate system.
+ * \param size Number of points along each axis (RA, DEC). The index of the
+ * origin of the coordinate system is (size / 2, size / 2).
+ * \param delta Angular step size in radians (assumed to be the same for
+ * both axes).
+ */
+ DirectionCoordinate makeCoordinates(const MDirection &reference,
+ unsigned int size, double delta);
+
+ /*!
+ * \brief Convert an ITRF position given as a StationResponse::vector3r_t
+ * instance to a casacore::MPosition.
+ */
+ MPosition toMPositionITRF(const vector3r_t &position);
+
+ /*!
+ * \brief Convert a casacore::MPosition instance to a
+ * StationResponse::vector3r_t instance.
+ */
+ vector3r_t fromMPosition(const MPosition &position);
+
+ /*!
+ * \brief Convert a casacore::MDirection instance to a
+ * StationResponse::vector3r_t instance.
+ */
+ vector3r_t fromMDirection(const MDirection &direction);
+
+ /*!
+ * \brief Check if the specified column exists as a column of the specified
+ * table.
+ *
+ * \param table The Table instance to check.
+ * \param column The name of the column.
+ */
+ bool hasColumn(const Table &table, const string &column);
+
+ /*!
+ * \brief Check if the specified sub-table exists as a sub-table of the
+ * specified table.
+ *
+ * \param table The Table instance to check.
+ * \param name The name of the sub-table.
+ */
+ bool hasSubTable(const Table &table, const string &name);
+
+ /*!
+ * \brief Provide access to a sub-table by name.
+ *
+ * \param table The Table instance to which the sub-table is associated.
+ * \param name The name of the sub-table.
+ */
+ Table getSubTable(const Table &table, const string &name);
+
+ /*!
+ * \brief Attempt to read the position of the observatory. If the
+ * observatory position is unknown, the specified default position is
+ * returned.
+ *
+ * \param ms MeasurementSet to read the observatory position from.
+ * \param idObservation Identifier that determines of which observation the
+ * observatory position should be read.
+ * \param defaultPosition The position that will be returned if the
+ * observatory position is unknown.
+ */
+ MPosition readObservatoryPosition(const MeasurementSet &ms,
+ unsigned int idObservation, const MPosition &defaultPosition);
+
+ /*!
+ * \brief Read the list of unique timestamps.
+ *
+ * \param ms MeasurementSet to read the list of unique timestamps from.
+ */
+ Vector<Double> readUniqueTimes(const MeasurementSet &ms);
+
+ /*!
+ * \brief Read the reference frequency of the subband associated to the
+ * specified data description identifier.
+ *
+ * \param ms MeasurementSet to read the reference frequency from.
+ * \param idDataDescription Identifier that determines of which subband the
+ * reference frequency should be read.
+ */
+ double readFreqReference(const MeasurementSet &ms,
+ unsigned int idDataDescription);
+
+ /*!
+ * \brief Read the phase reference direction.
+ *
+ * \param ms MeasurementSet to read the phase reference direction from.
+ * \param idField Identifier of the field of which the phase reference
+ * direction should be read.
+ */
+ MDirection readPhaseReference(const MeasurementSet &ms,
+ unsigned int idField);
+
+ /*!
+ * \brief Read the station beam former reference direction.
+ *
+ * \param ms MeasurementSet to read the station beam former reference
+ * direction from.
+ * \param idField Identifier of the field of which the station beam former
+ * reference direction should be read.
+ */
+ MDirection readDelayReference(const MeasurementSet &ms,
+ unsigned int idField);
+
+ /*!
+ * \brief Read the station beam former reference direction.
+ *
+ * \param ms MeasurementSet to read the tile beam former reference
+ * direction from.
+ * \param idField Identifier of the field of which the tile beam former
+ * reference direction should be read.
+ */
+ MDirection readTileReference(const MeasurementSet &ms,
+ unsigned int idField);
+
+ /*!
+ * \brief Store the specified cube of pixel data as a CASA image.
+ *
+ * \param data Pixel data. The third dimension is assumed to be of length
+ * 4, referring to the correlation products XX, XY, YX, YY (in this order).
+ * \param coordinates Sky coordinate system definition.
+ * \param frequency Frequency for which the pixel data is valid (Hz).
+ * \param name File name of the output image.
+ */
+ template <class T>
+ void store(const Cube<T> &data, const DirectionCoordinate &coordinates,
+ double frequency, const string &name);
+
+ /*!
+ * \brief Convert a string to a CASA Quantity (value with unit).
+ */
+ Quantity readQuantity(const String &in);
+
+ /*!
+ * \brief Remove all elements from the range [first, last) that fall
+ * outside the interval [min, max].
+ *
+ * This function returns an iterator new_last such that the range [first,
+ * new_last) contains no elements that fall outside the interval [min,
+ * max]. The iterators in the range [new_last, last) are all still
+ * dereferenceable, but the elements that they point to are unspecified.
+ * The order of the elements that are not removed is unchanged.
+ */
+ template <typename T>
+ T filter(T first, T last, int min, int max);
+} //# unnamed namespace
+
+
+int main(int argc, char *argv[])
+{
+ TEST_SHOW_VERSION(argc, argv, StationResponse);
+ INIT_LOGGER(basename(string(argv[0])));
+ Version::show<StationResponseVersion>(cout);
+
+ // Parse inputs.
+ LOFAR::InputParSet inputs;
+ inputs.create("ms", "", "Name of input MeasurementSet", "string");
+ inputs.create("stations", "0", "IDs of stations to process", "int vector");
+ inputs.create("cellsize", "60arcsec", "Angular pixel size",
+ "quantity string");
+ inputs.create("size", "256", "Number of pixels along each axis", "int");
+ inputs.create("offset", "0s", "Time offset from the start of the MS",
+ "quantity string");
+ inputs.create("frames", "0", "Number of images that will be generated for"
+ " each station (equally spaced over the duration of the MS)", "int");
+ inputs.create("abs", "false", "If set to true, store the absolute value of"
+ " the beam response instead of the complex value (intended for use with"
+ " older versions of casaviewer)", "bool");
+ inputs.readArguments(argc, argv);
+
+ vector<int> stationIDs(inputs.getIntVector("stations"));
+ Quantity cellsize = readQuantity(inputs.getString("cellsize"));
+ unsigned int size = max(inputs.getInt("size"), 1);
+ Quantity offset = readQuantity(inputs.getString("offset"));
+ unsigned int nFrames = max(inputs.getInt("frames"), 1);
+ Bool abs = inputs.getBool("abs");
+
+ // Open MS.
+ MeasurementSet ms(inputs.getString("ms"));
+ unsigned int idObservation = 0, idField = 0, idDataDescription = 0;
+
+ // Read station meta-data.
+ vector<Station::Ptr> stations;
+ readStations(ms, std::back_inserter(stations));
+
+ // Remove illegal station indices.
+ stationIDs.erase(filter(stationIDs.begin(), stationIDs.end(), 0,
+ static_cast<int>(stations.size()) - 1), stationIDs.end());
+
+ // Read unique timestamps
+ Table selection =
+ ms(ms.col("OBSERVATION_ID") == static_cast<Int>(idObservation)
+ && ms.col("FIELD_ID") == static_cast<Int>(idField)
+ && ms.col("DATA_DESC_ID") == static_cast<Int>(idDataDescription));
+ Vector<Double> time = readUniqueTimes(selection);
+
+ // Read reference frequency.
+ double refFrequency = readFreqReference(ms, idDataDescription);
+
+ // Use the position of the first selected station as the array reference
+ // position if the observatory position cannot be found.
+ MPosition refPosition = readObservatoryPosition(ms, idField,
+ toMPositionITRF(stations.front()->position()));
+
+ // Read phase reference direction.
+ MDirection refPhase = readPhaseReference(ms, idField);
+
+ // Read delay reference direction.
+ MDirection refDelay = readDelayReference(ms, idField);
+
+ // Read tile reference direction.
+ MDirection refTile = readTileReference(ms, idField);
+
+ // Create image coordinate system.
+ IPosition shape(2, size, size);
+ DirectionCoordinate coordinates = makeCoordinates(refPhase, size,
+ cellsize.getValue("rad"));
+
+ // Compute station response images.
+ Cube<Complex> response(size, size, 4);
+
+ MEpoch refEpoch;
+ Quantity refTime(time(0), "s");
+ refTime = refTime + offset;
+ Quantity deltaTime((time(time.size() - 1) - time(0) - offset.getValue("s"))
+ / (nFrames - 1), "s");
+
+ for(size_t j = 0; j < nFrames; ++j)
+ {
+ cout << "[ Frame: " << (j + 1) << "/" << nFrames << " Offset: +"
+ << refTime.getValue() - time(0) << " s ]" << endl;
+
+ // Update reference epoch.
+ refEpoch.set(refTime);
+
+ cout << "Creating ITRF direction map... " << flush;
+ ITRFDirectionMap directionMap = makeDirectionMap(coordinates, shape,
+ refEpoch, refPosition, refDelay, refTile);
+ cout << "done." << endl;
+
+ cout << "Computing response images... " << flush;
+ for(vector<int>::const_iterator it = stationIDs.begin(),
+ end = stationIDs.end(); it != end; ++it)
+ {
+ Station::ConstPtr station = stations[*it];
+ cout << *it << ":" << station->name() << " " << flush;
+
+ for(size_t y = 0; y < size; ++y)
+ {
+ for(size_t x = 0; x < size; ++x)
+ {
+ matrix22c_t E = station->response(directionMap.time0,
+ refFrequency, directionMap.directions(x, y),
+ refFrequency, directionMap.station0,
+ directionMap.tile0);
+
+ response(x, y, 0) = E[0][0];
+ response(x, y, 1) = E[0][1];
+ response(x, y, 2) = E[1][0];
+ response(x, y, 3) = E[1][1];
+ }
+ }
+
+ std::ostringstream oss;
+ oss << "response-" << station->name() << "-frame-" << (j + 1)
+ << ".img";
+
+ if(abs)
+ {
+ store(Cube<Float>(amplitude(response)), coordinates,
+ refFrequency, oss.str());
+ }
+ else
+ {
+ store(response, coordinates, refFrequency, oss.str());
+ }
+ }
+ cout << endl;
+
+ refTime = refTime + deltaTime;
+ }
+
+ return 0;
+}
+
+namespace
+{
+ ITRFDirectionMap makeDirectionMap(const DirectionCoordinate &coordinates,
+ const IPosition &shape, const MEpoch &epoch, const MPosition &position0,
+ const MDirection &station0, const MDirection &tile0)
+ {
+ ITRFDirectionMap map;
+
+ // Convert from MEpoch to a time in MJD(UTC) in seconds.
+ MEpoch mEpochUTC = MEpoch::Convert(epoch, MEpoch::Ref(MEpoch::UTC))();
+ MVEpoch mvEpochUTC = mEpochUTC.getValue();
+ Quantity qEpochUTC = mvEpochUTC.getTime();
+ map.time0 = qEpochUTC.getValue("s");
+
+ // Create conversion engine J2000 => ITRF at the specified epoch.
+ MDirection::Convert convertor = MDirection::Convert(MDirection::J2000,
+ MDirection::Ref(MDirection::ITRF, MeasFrame(epoch, position0)));
+
+ // Compute station and tile beam former reference directions in ITRF at
+ // the specified epoch.
+ map.station0 = fromMDirection(convertor(station0));
+ map.tile0 = fromMDirection(convertor(tile0));
+
+ // Pre-allocate space for the grid of ITRF directions.
+ map.directions.resize(shape);
+
+ // Compute ITRF directions.
+ MDirection world;
+ Vector<Double> pixel = coordinates.referencePixel();
+ for(pixel(1) = 0.0; pixel(1) < shape(1); ++pixel(1))
+ {
+ for(pixel(0) = 0.0; pixel(0) < shape(0); ++pixel(0))
+ {
+ // CoordinateSystem::toWorld(): RA range [-pi,pi], DEC range
+ // [-pi/2,pi/2].
+ if(coordinates.toWorld(world, pixel))
+ {
+ map.directions(pixel(0), pixel(1)) =
+ fromMDirection(convertor(world));
+ }
+ }
+ }
+
+ return map;
+ }
+
+ DirectionCoordinate makeCoordinates(const MDirection &reference,
+ unsigned int size, double delta)
+ {
+ MDirection referenceJ2000 = MDirection::Convert(reference,
+ MDirection::J2000)();
+ Quantum<Vector<Double> > referenceAngles = referenceJ2000.getAngle();
+ double ra = referenceAngles.getBaseValue()(0);
+ double dec = referenceAngles.getBaseValue()(1);
+
+ Matrix<Double> xform(2,2);
+ xform = 0.0; xform.diagonal() = 1.0;
+ return DirectionCoordinate(MDirection::J2000,
+ Projection(Projection::SIN), ra, dec, -delta, delta, xform,
+ size / 2, size / 2);
+ }
+
+ MPosition toMPositionITRF(const vector3r_t &position)
+ {
+ MVPosition mvITRF(position[0], position[1], position[2]);
+ return MPosition(mvITRF, MPosition::ITRF);
+ }
+
+ vector3r_t fromMPosition(const MPosition &position)
+ {
+ MVPosition mvPosition = position.getValue();
+ vector3r_t result = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
+ return result;
+ }
+
+ vector3r_t fromMDirection(const MDirection &direction)
+ {
+ MVDirection mvDirection = direction.getValue();
+ vector3r_t result = {{mvDirection(0), mvDirection(1), mvDirection(2)}};
+ return result;
+ }
+
+ bool hasColumn(const Table &table, const string &column)
+ {
+ return table.tableDesc().isColumn(column);
+ }
+
+ bool hasSubTable(const Table &table, const string &name)
+ {
+ return table.keywordSet().isDefined(name);
+ }
+
+ Table getSubTable(const Table &table, const string &name)
+ {
+ return table.keywordSet().asTable(name);
+ }
+
+ MPosition readObservatoryPosition(const MeasurementSet &ms,
+ unsigned int idObservation, const MPosition &defaultPosition)
+ {
+ // Get the instrument position in ITRF coordinates, or use the centroid
+ // of the station positions if the instrument position is unknown.
+ ROMSObservationColumns observation(ms.observation());
+ ASSERT(observation.nrow() > idObservation);
+ ASSERT(!observation.flagRow()(idObservation));
+
+ // Read observatory name and try to look-up its position.
+ const string observatory = observation.telescopeName()(idObservation);
+
+ // Look-up observatory position, default to specified default position.
+ MPosition position(defaultPosition);
+ MeasTable::Observatory(position, observatory);
+ return position;
+ }
+
+ Vector<Double> readUniqueTimes(const MeasurementSet &ms)
+ {
+ Table tab_sorted = ms.sort("TIME", Sort::Ascending, Sort::HeapSort
+ | Sort::NoDuplicates);
+
+ ROScalarColumn<Double> c_time(tab_sorted, "TIME");
+ return c_time.getColumn();
+ }
+
+ double readFreqReference(const MeasurementSet &ms,
+ unsigned int idDataDescription)
+ {
+ ROMSDataDescColumns desc(ms.dataDescription());
+ ASSERT(desc.nrow() > idDataDescription);
+ ASSERT(!desc.flagRow()(idDataDescription));
+ uInt idWindow = desc.spectralWindowId()(idDataDescription);
+
+ ROMSSpWindowColumns window(ms.spectralWindow());
+ ASSERT(window.nrow() > idWindow);
+ ASSERT(!window.flagRow()(idWindow));
+
+ return window.refFrequency()(idWindow);
+ }
+
+ MDirection readPhaseReference(const MeasurementSet &ms,
+ unsigned int idField)
+ {
+ ROMSFieldColumns field(ms.field());
+ ASSERT(field.nrow() > idField);
+ ASSERT(!field.flagRow()(idField));
+
+ return field.phaseDirMeas(idField);
+ }
+
+ MDirection readDelayReference(const MeasurementSet &ms,
+ unsigned int idField)
+ {
+ ROMSFieldColumns field(ms.field());
+ ASSERT(field.nrow() > idField);
+ ASSERT(!field.flagRow()(idField));
+
+ return field.delayDirMeas(idField);
+ }
+
+ MDirection readTileReference(const MeasurementSet &ms, unsigned int idField)
+ {
+ // The MeasurementSet class does not support LOFAR specific columns, so
+ // we use ROArrayMeasColumn to read the tile beam reference direction.
+ Table tab_field = getSubTable(ms, "FIELD");
+
+ static const String columnName = "LOFAR_TILE_BEAM_DIR";
+ if(hasColumn(tab_field, columnName))
+ {
+ ROArrayMeasColumn<MDirection> c_direction(tab_field, columnName);
+ if(c_direction.isDefined(idField))
+ {
+ return c_direction(idField)(IPosition(1, 0));
+ }
+ }
+
+ // By default, the tile beam reference direction is assumed to be equal
+ // to the station beam reference direction (for backward compatibility,
+ // and for non-HBA measurements).
+ return readDelayReference(ms, idField);
+ }
+
+ template <class T>
+ void store(const Cube<T> &data, const DirectionCoordinate &coordinates,
+ double frequency, const string &name)
+ {
+ ASSERT(data.shape()(2) == 4);
+
+ Vector<Int> stokes(4);
+ stokes(0) = Stokes::XX;
+ stokes(1) = Stokes::XY;
+ stokes(2) = Stokes::YX;
+ stokes(3) = Stokes::YY;
+
+ CoordinateSystem csys;
+ csys.addCoordinate(coordinates);
+ csys.addCoordinate(StokesCoordinate(stokes));
+ csys.addCoordinate(SpectralCoordinate(MFrequency::TOPO, frequency, 0.0,
+ 0.0));
+
+ PagedImage<T> im(TiledShape(IPosition(4, data.shape()(0),
+ data.shape()(1), 4, 1)), csys, name);
+ im.putSlice(data, IPosition(4, 0, 0, 0, 0));
+ }
+
+ Quantity readQuantity(const String &in)
+ {
+ Quantity result;
+ bool status = Quantity::read(result, in);
+ ASSERT(status);
+ return result;
+ }
+
+ template <typename T>
+ T filter(T first, T last, int min, int max)
+ {
+ T new_last = first;
+ for(; first != last; ++first)
+ {
+ if(*first >= min && *first <= max)
+ {
+ *new_last++ = *first;
+ }
+ }
+
+ return new_last;
+ }
+} // unnamed namespace.