diff --git a/Antenna.h b/Antenna.h
index 087ed2b93de9b1f55f9c5012e302716475e4d44c..0d1f69ded126671ba757b7810369703d6c9aeec7 100644
--- a/Antenna.h
+++ b/Antenna.h
@@ -10,194 +10,167 @@
namespace everybeam {
/**
- * @brief (Virtual) class describing an antenna, and computing the corresponding response()
- * and arrayFactor().
- *
+ * @brief (Virtual) class describing an antenna, and computing the corresponding
+ * response() and arrayFactor().
+ *
*/
-class Antenna
-{
-public:
-
- /**
- * \brief %Station 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;
-
- constexpr static Axes identity_axes = {
- {1.0, 0.0, 0.0},
- {0.0, 1.0, 0.0},
- {0.0, 0.0, 1.0}
- };
-
- constexpr static vector3r_t zero_origin = {0.0, 0.0, 0.0};
+class Antenna {
+ public:
+ /**
+ * \brief %Station 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;
};
-
- constexpr static CoordinateSystem identity_coordinate_system{
- CoordinateSystem::zero_origin,
- CoordinateSystem::identity_axes
- };
-
- typedef std::shared_ptr<Antenna> Ptr;
-
- /**
- * @brief Struct containing antenna options
- *
- */
- struct Options
- {
- real_t freq0; //!< %Antenna reference frequency (Hz).
- vector3r_t station0; //!< Reference direction (ITRF, m)
- vector3r_t tile0; //!< Tile beam former reference direction (ITRF, m).
- bool rotate; //!< Boolean deciding if paralactic rotation should be applied.
- vector3r_t east; //!< Eastward pointing unit vector
- vector3r_t north; //!< Northward pointing unit vector
- };
-
- /**
- * @brief Construct a new %Antenna object
- *
- */
- Antenna() :
- // default coordinate system
- // no shift of origin, no rotation
- Antenna({
- CoordinateSystem::zero_origin, // origin
- CoordinateSystem::identity_axes
- })
- {}
-
- /**
- * @brief Construct a new %Antenna object, given a coordinate system
- *
- * @param coordinate_system
- */
- Antenna(const CoordinateSystem &coordinate_system) :
- // default phase reference system is the origin of the coordinate system
- Antenna(coordinate_system, coordinate_system.origin)
- {}
-
- /**
- * @brief Construct a new %Antenna object, given a coordinate system and a
- * phase reference position.
- *
- * @param coordinate_system Coordinate system
- * @param phase_reference_position Phase reference position
- */
- Antenna(const CoordinateSystem &coordinate_system, const vector3r_t &phase_reference_position) :
- m_coordinate_system(coordinate_system),
+ vector3r_t origin;
+ Axes axes;
+
+ constexpr static Axes identity_axes = {
+ {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}};
+
+ constexpr static vector3r_t zero_origin = {0.0, 0.0, 0.0};
+ };
+
+ constexpr static CoordinateSystem identity_coordinate_system{
+ CoordinateSystem::zero_origin, CoordinateSystem::identity_axes};
+
+ typedef std::shared_ptr<Antenna> Ptr;
+
+ /**
+ * @brief Struct containing antenna options
+ *
+ */
+ struct Options {
+ real_t freq0; //!< %Antenna reference frequency (Hz).
+ vector3r_t station0; //!< Reference direction (ITRF, m)
+ vector3r_t tile0; //!< Tile beam former reference direction (ITRF, m).
+ bool
+ rotate; //!< Boolean deciding if paralactic rotation should be applied.
+ vector3r_t east; //!< Eastward pointing unit vector
+ vector3r_t north; //!< Northward pointing unit vector
+ };
+
+ /**
+ * @brief Construct a new %Antenna object
+ *
+ */
+ Antenna()
+ : // default coordinate system
+ // no shift of origin, no rotation
+ Antenna({CoordinateSystem::zero_origin, // origin
+ CoordinateSystem::identity_axes}) {}
+
+ /**
+ * @brief Construct a new %Antenna object, given a coordinate system
+ *
+ * @param coordinate_system
+ */
+ Antenna(const CoordinateSystem &coordinate_system)
+ : // default phase reference system is the origin of the coordinate
+ // system
+ Antenna(coordinate_system, coordinate_system.origin) {}
+
+ /**
+ * @brief Construct a new %Antenna object, given a coordinate system and a
+ * phase reference position.
+ *
+ * @param coordinate_system Coordinate system
+ * @param phase_reference_position Phase reference position
+ */
+ Antenna(const CoordinateSystem &coordinate_system,
+ const vector3r_t &phase_reference_position)
+ : m_coordinate_system(coordinate_system),
m_phase_reference_position(phase_reference_position),
- m_enabled{true, true}
- {
- }
-
- /**
- * @brief Compute the %Antenna response
- *
- * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
- * @param freq Frequency of the plane wave (Hz).
- * @param direction Direction of arrival (ITRF, m).
- * @param options
- * @return matrix22c_t Jones matrix
- */
- matrix22c_t response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options = {})
- {
- // Transform direction and directions in options to local coordinatesystem
- vector3r_t local_direction = transform_to_local_direction(direction);
- Options local_options = {
- .freq0 = options.freq0,
- .station0 = transform_to_local_direction(options.station0),
- .tile0 = transform_to_local_direction(options.tile0),
- .rotate = options.rotate,
- .east = transform_to_local_direction(options.east),
- .north = transform_to_local_direction(options.north)
- };
- matrix22c_t response = local_response(time, freq, local_direction, local_options);
- return response;
- }
-
- /**
- * @brief Compute the array factor of the antenna
- *
- * @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 options
- * @return diag22c_t
- */
- diag22c_t arrayFactor(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options = {})
- {
- // Transform direction and directions in options to local coordinatesystem
- vector3r_t local_direction = transform_to_local_direction(direction);
- Options local_options = {
- .freq0 = options.freq0,
- .station0 = transform_to_local_direction(options.station0),
- .tile0 = transform_to_local_direction(options.tile0)
- };
- return local_arrayFactor(time, freq, local_direction, local_options);
- }
-
- CoordinateSystem m_coordinate_system;
- vector3r_t m_phase_reference_position;
- bool m_enabled[2];
-
-private:
-
- virtual matrix22c_t local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const = 0;
-
- virtual diag22c_t local_arrayFactor(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const
- { return {1.0, 1.0}; }
-
- vector3r_t transform_to_local_direction(const vector3r_t &direction);
-
+ m_enabled{true, true} {}
+
+ /**
+ * @brief Compute the %Antenna response
+ *
+ * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * @param freq Frequency of the plane wave (Hz).
+ * @param direction Direction of arrival (ITRF, m).
+ * @param options
+ * @return matrix22c_t Jones matrix
+ */
+ matrix22c_t response(real_t time, real_t freq, const vector3r_t &direction,
+ const Options &options = {}) {
+ // Transform direction and directions in options to local coordinatesystem
+ vector3r_t local_direction = transform_to_local_direction(direction);
+ Options local_options = {
+ .freq0 = options.freq0,
+ .station0 = transform_to_local_direction(options.station0),
+ .tile0 = transform_to_local_direction(options.tile0),
+ .rotate = options.rotate,
+ .east = transform_to_local_direction(options.east),
+ .north = transform_to_local_direction(options.north)};
+ matrix22c_t response =
+ local_response(time, freq, local_direction, local_options);
+ return response;
+ }
+
+ /**
+ * @brief Compute the array factor of the antenna
+ *
+ * @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 options
+ * @return diag22c_t
+ */
+ diag22c_t arrayFactor(real_t time, real_t freq, const vector3r_t &direction,
+ const Options &options = {}) {
+ // Transform direction and directions in options to local coordinatesystem
+ vector3r_t local_direction = transform_to_local_direction(direction);
+ Options local_options = {
+ .freq0 = options.freq0,
+ .station0 = transform_to_local_direction(options.station0),
+ .tile0 = transform_to_local_direction(options.tile0)};
+ return local_arrayFactor(time, freq, local_direction, local_options);
+ }
+
+ CoordinateSystem m_coordinate_system;
+ vector3r_t m_phase_reference_position;
+ bool m_enabled[2];
+
+ private:
+ virtual matrix22c_t local_response(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const Options &options) const = 0;
+
+ virtual diag22c_t local_arrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const Options &options) const {
+ return {1.0, 1.0};
+ }
+
+ vector3r_t transform_to_local_direction(const vector3r_t &direction);
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/BeamFormer.h b/BeamFormer.h
index 8fe72d32db622d05c06fcb0f4a1468916f392038..38b1bb1517b59a9634976db706d27ef30c9df6b8 100644
--- a/BeamFormer.h
+++ b/BeamFormer.h
@@ -8,88 +8,87 @@
#include "Types.h"
namespace everybeam {
-class BeamFormer : public Antenna
-{
-public:
+class BeamFormer : public Antenna {
+ public:
+ typedef std::shared_ptr<BeamFormer> Ptr;
- typedef std::shared_ptr<BeamFormer> Ptr;
+ /**
+ * @brief Construct a new BeamFormer object
+ *
+ */
+ BeamFormer() : Antenna() {
+ m_local_phase_reference_position =
+ transform_to_local_position(m_phase_reference_position);
+ }
- /**
- * @brief Construct a new BeamFormer object
- *
- */
- BeamFormer() :
- Antenna()
- {
- m_local_phase_reference_position = transform_to_local_position(m_phase_reference_position);
- }
+ /**
+ * @brief Construct a new BeamFormer object given a coordinate system.
+ *
+ * @param coordinate_system
+ */
+ BeamFormer(const CoordinateSystem &coordinate_system)
+ : Antenna(coordinate_system) {
+ m_local_phase_reference_position =
+ transform_to_local_position(m_phase_reference_position);
+ }
- /**
- * @brief Construct a new BeamFormer object given a coordinate system.
- *
- * @param coordinate_system
- */
- BeamFormer(const CoordinateSystem &coordinate_system) :
- Antenna(coordinate_system)
- {
- m_local_phase_reference_position = transform_to_local_position(m_phase_reference_position);
- }
+ /**
+ * @brief Construct a new BeamFormer object given a coordinate system and a
+ * phase reference position
+ *
+ * @param coordinate_system
+ * @param phase_reference_position
+ */
+ BeamFormer(CoordinateSystem coordinate_system,
+ vector3r_t phase_reference_position)
+ : Antenna(coordinate_system, phase_reference_position) {
+ m_local_phase_reference_position =
+ transform_to_local_position(m_phase_reference_position);
+ }
- /**
- * @brief Construct a new BeamFormer object given a coordinate system and a phase reference position
- *
- * @param coordinate_system
- * @param phase_reference_position
- */
- BeamFormer(CoordinateSystem coordinate_system, vector3r_t phase_reference_position) :
- Antenna(coordinate_system, phase_reference_position)
- {
- m_local_phase_reference_position = transform_to_local_position(m_phase_reference_position);
- }
+ /**
+ * @brief Add an antenna to the m_antenna array.
+ *
+ * @param antenna
+ */
+ void add_antenna(Antenna::Ptr antenna) { m_antennas.push_back(antenna); }
- /**
- * @brief Add an antenna to the m_antenna array.
- *
- * @param antenna
- */
- void add_antenna(Antenna::Ptr antenna) {m_antennas.push_back(antenna);}
+ private:
+ vector3r_t
+ m_local_phase_reference_position; // in coordinate system of Antenna
-private:
+ // Transform position vector into a local position vector
+ vector3r_t transform_to_local_position(const vector3r_t &position);
- vector3r_t m_local_phase_reference_position; // in coordinate system of Antenna
+ // Compute the BeamFormer response in certain direction of arrival (ITRF, m)
+ // and return (Jones) matrix of response
+ virtual matrix22c_t local_response(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const Options &options) const override;
- // Transform position vector into a local position vector
- vector3r_t transform_to_local_position(const vector3r_t &position);
+ // Compute the local arrayFactor, with arrayFactor a vectorial
+ // "representation" of Jones matrix
+ virtual diag22c_t local_arrayFactor(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const Options &options) const override {
+ return {1.0, 1.0};
+ }
- // Compute the BeamFormer response in certain direction of arrival (ITRF, m)
- // and return (Jones) matrix of response
- virtual matrix22c_t local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const override;
+ // Compute the geometric response for all the antennas in the BeamFormer based
+ // on the probing frequency and a specified direction (either pointing dir or
+ // dir of interest).
+ std::vector<std::complex<double>> compute_geometric_response(
+ const double freq, const vector3r_t &direction) const;
- // Compute the local arrayFactor, with arrayFactor a vectorial "representation"
- // of Jones matrix
- virtual diag22c_t local_arrayFactor(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const override
- {
- return {1.0, 1.0};
- }
-
- // Compute the geometric response for all the antennas in the BeamFormer based on
- // the probing frequency and a specified direction (either pointing dir or dir of interest).
- std::vector<std::complex<double>> compute_geometric_response(const double freq, const vector3r_t &direction) const;
-
- // Compute the weights based on the pointing direction of the beam and the beam reference frequence.
- std::vector<std::pair<std::complex<double>,std::complex<double>>> compute_weights(const vector3r_t &direction, double freq) const;
+ // Compute the weights based on the pointing direction of the beam and the
+ // beam reference frequence.
+ std::vector<std::pair<std::complex<double>, std::complex<double>>>
+ compute_weights(const vector3r_t &direction, double freq) const;
- // List of antennas in BeamFormer
- // TODO: Maybe refactor to _m_antennas to indicate m_antennas is a private attribute
- std::vector<Antenna::Ptr> m_antennas;
+ // List of antennas in BeamFormer
+ // TODO: Maybe refactor to _m_antennas to indicate m_antennas is a private
+ // attribute
+ std::vector<Antenna::Ptr> m_antennas;
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/Constants.h b/Constants.h
index 53bf3eab7b4fb5c4c315f270dfe2c7043c69bc04..e335df7f13478782d70d6325fd1e31e9cd103819 100644
--- a/Constants.h
+++ b/Constants.h
@@ -31,11 +31,10 @@
namespace everybeam {
/** %Constants used in this library. */
-namespace constants
-{
+namespace constants {
/** Speed of light (m/s) */
const real_t c = 2.99792458e+08;
-} // namespace constants
-} // namespace everybeam
+} // namespace constants
+} // namespace everybeam
#endif
diff --git a/Element.h b/Element.h
index 83ba0ae0d843ca4a451824b70111cf25ab6b5192..ad5acf35ea8cce7e84978051291db94f4cfeb8f8 100644
--- a/Element.h
+++ b/Element.h
@@ -13,55 +13,48 @@ namespace everybeam {
/**
* @brief Elementary antenna, for which a response can be computed,
* but without any substructure like a beamformer
- *
+ *
*/
-class Element : public Antenna
-{
-public:
-
- typedef std::shared_ptr<Element> Ptr;
-
- /**
- * @brief Construct a new Element object
- *
- * @param coordinate_system (antenna) CoordinateSystem
- * @param element_response ElementResponseModel
- * @param id
- */
- Element(const CoordinateSystem &coordinate_system, ElementResponse::Ptr element_response, int id) :
- Antenna(coordinate_system),
+class Element : public Antenna {
+ public:
+ typedef std::shared_ptr<Element> Ptr;
+
+ /**
+ * @brief Construct a new Element object
+ *
+ * @param coordinate_system (antenna) CoordinateSystem
+ * @param element_response ElementResponseModel
+ * @param id
+ */
+ Element(const CoordinateSystem &coordinate_system,
+ ElementResponse::Ptr element_response, int id)
+ : Antenna(coordinate_system),
m_id(id),
- m_element_response(element_response)
- {}
-
- /**
- * @brief Compute the local response of the element.
- *
- * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
- * @param freq Frequency of the plane wave (Hz).
- * @param direction Direction of arrival (ITRF, m).
- * @param id ID of element
- * @param options
- * @return matrix22c_t
- */
- matrix22c_t local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- size_t id,
- const Options &options) const;
-
-private:
- virtual matrix22c_t local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const final override;
-
- int m_id;
- ElementResponse::Ptr m_element_response;
+ m_element_response(element_response) {}
+
+ /**
+ * @brief Compute the local response of the element.
+ *
+ * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * @param freq Frequency of the plane wave (Hz).
+ * @param direction Direction of arrival (ITRF, m).
+ * @param id ID of element
+ * @param options
+ * @return matrix22c_t
+ */
+ matrix22c_t local_response(real_t time, real_t freq,
+ const vector3r_t &direction, size_t id,
+ const Options &options) const;
+
+ private:
+ virtual matrix22c_t local_response(
+ real_t time, real_t freq, const vector3r_t &direction,
+ const Options &options) const final override;
+
+ int m_id;
+ ElementResponse::Ptr m_element_response;
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/ElementResponse.h b/ElementResponse.h
index eca9f106711237e4d8c38d46a42b4840b2cd2220..8ac85f1557f0efe82e527cb294d646ed2eec9ddd 100644
--- a/ElementResponse.h
+++ b/ElementResponse.h
@@ -9,60 +9,49 @@
namespace everybeam {
enum ElementResponseModel {
- Unknown,
- Hamaker,
- LOBES,
- OSKARDipole,
- OSKARSphericalWave
+ Unknown,
+ Hamaker,
+ LOBES,
+ OSKARDipole,
+ OSKARSphericalWave
};
std::ostream& operator<<(std::ostream& os, ElementResponseModel model);
/**
- * @brief Virtual class for the element response model. All the (antenna/element)
- * response models inherit from this class.
- *
+ * @brief Virtual class for the element response model. All the
+ * (antenna/element) response models inherit from this class.
+ *
*/
-class ElementResponse
-{
-public:
-
- typedef MutablePtr<ElementResponse> Ptr; //!< Pointer to ElementResponse object
-
- /**
- * @brief Virtual implementation of response method
- *
- * @param freq Frequency of the plane wave (Hz).
- * @param theta Angle wrt. z-axis (rad)
- * @param phi Angle in the xy-plane wrt. x-axis (rad)
- * @param result Pointer to 2x2 array of Jones matrix
- */
- virtual void response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&result)[2][2]) const = 0;
-
-
- /**
- * @brief Virtual implementation of response method
- *
- * @param element_id ID of element
- * @param freq Frequency of the plane wave (Hz).
- * @param theta Angle wrt. z-axis (rad)
- * @param phi Angle in the xy-plane wrt. x-axis (rad)
- * @param result Pointer to 2x2 array of Jones matrix
- */
- virtual void response(
- int element_id,
- double freq,
- double theta,
- double phi,
- std::complex<double> (&result)[2][2]) const
- {
- response(freq, theta, phi, result);
- }
-
+class ElementResponse {
+ public:
+ typedef MutablePtr<ElementResponse>
+ Ptr; //!< Pointer to ElementResponse object
+
+ /**
+ * @brief Virtual implementation of response method
+ *
+ * @param freq Frequency of the plane wave (Hz).
+ * @param theta Angle wrt. z-axis (rad)
+ * @param phi Angle in the xy-plane wrt. x-axis (rad)
+ * @param result Pointer to 2x2 array of Jones matrix
+ */
+ virtual void response(double freq, double theta, double phi,
+ std::complex<double> (&result)[2][2]) const = 0;
+
+ /**
+ * @brief Virtual implementation of response method
+ *
+ * @param element_id ID of element
+ * @param freq Frequency of the plane wave (Hz).
+ * @param theta Angle wrt. z-axis (rad)
+ * @param phi Angle in the xy-plane wrt. x-axis (rad)
+ * @param result Pointer to 2x2 array of Jones matrix
+ */
+ virtual void response(int element_id, double freq, double theta, double phi,
+ std::complex<double> (&result)[2][2]) const {
+ response(freq, theta, phi, result);
+ }
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/ITRFConverter.h b/ITRFConverter.h
index ff830b5e264a05d1c6cc980c5b4afbaa75007ba6..706d6ead515b2402417dab01b65d9fe50d3cbe64 100644
--- a/ITRFConverter.h
+++ b/ITRFConverter.h
@@ -17,29 +17,28 @@
namespace everybeam {
/**
- * @brief Class providing utilities for coordinate transformations
- * to and from ITRF (International Terrestrial Reference Frame)
- *
+ * @brief Class providing utilities for coordinate transformations
+ * to and from ITRF (International Terrestrial Reference Frame)
+ *
*/
-class ITRFConverter
-{
-public:
- typedef std::unique_ptr<ITRFDirection> Ptr;
- typedef std::unique_ptr<const ITRFDirection> ConstPtr;
-
- ITRFConverter(real_t time);
-
- void setTime(real_t time);
- vector3r_t j2000ToITRF(const vector2r_t &j2000Direction) const;
- vector3r_t j2000ToITRF(const vector3r_t &j2000Direction) const;
- vector3r_t toITRF(const casacore::MDirection &direction) const;
- casacore::MDirection toDirection(const vector2r_t &j2000Direction) const;
- casacore::MDirection toDirection(const vector3r_t &j2000Direction) const;
- casacore::MDirection toDirection(const casacore::MDirection &direction) const;
-
-private:
- casacore::MeasFrame itsFrame;
- mutable casacore::MDirection::Convert itsConverter;
+class ITRFConverter {
+ public:
+ typedef std::unique_ptr<ITRFDirection> Ptr;
+ typedef std::unique_ptr<const ITRFDirection> ConstPtr;
+
+ ITRFConverter(real_t time);
+
+ void setTime(real_t time);
+ vector3r_t j2000ToITRF(const vector2r_t &j2000Direction) const;
+ vector3r_t j2000ToITRF(const vector3r_t &j2000Direction) const;
+ vector3r_t toITRF(const casacore::MDirection &direction) const;
+ casacore::MDirection toDirection(const vector2r_t &j2000Direction) const;
+ casacore::MDirection toDirection(const vector3r_t &j2000Direction) const;
+ casacore::MDirection toDirection(const casacore::MDirection &direction) const;
+
+ private:
+ casacore::MeasFrame itsFrame;
+ mutable casacore::MDirection::Convert itsConverter;
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/ITRFDirection.h b/ITRFDirection.h
index 1a2e1428edded7976f846c6e29ca547f681ac680..4dc8e43b6942751d1bc49bbd977cbe8483c0822d 100644
--- a/ITRFDirection.h
+++ b/ITRFDirection.h
@@ -37,30 +37,29 @@
namespace everybeam {
-class ITRFDirection
-{
-public:
- typedef std::shared_ptr<ITRFDirection> Ptr;
- typedef std::shared_ptr<const ITRFDirection> ConstPtr;
+class ITRFDirection {
+ public:
+ typedef std::shared_ptr<ITRFDirection> Ptr;
+ typedef std::shared_ptr<const ITRFDirection> ConstPtr;
- ITRFDirection(const vector3r_t &position, const vector2r_t &direction);
- ITRFDirection(const vector3r_t &position, const vector3r_t &direction);
- ITRFDirection(const vector2r_t &direction);
- ITRFDirection(const vector3r_t &direction);
+ ITRFDirection(const vector3r_t &position, const vector2r_t &direction);
+ ITRFDirection(const vector3r_t &position, const vector3r_t &direction);
+ ITRFDirection(const vector2r_t &direction);
+ ITRFDirection(const vector3r_t &direction);
- vector3r_t at(real_t time) const;
+ vector3r_t at(real_t time) const;
- const static vector3r_t& LOFARPosition() { return itsLOFARPosition; }
+ const static vector3r_t &LOFARPosition() { return itsLOFARPosition; }
-private:
- //ITRF position of CS002LBA, just to use a fixed reference
- const static vector3r_t itsLOFARPosition;
+ private:
+ // ITRF position of CS002LBA, just to use a fixed reference
+ const static vector3r_t itsLOFARPosition;
- mutable casacore::MeasFrame itsFrame;
- mutable casacore::MDirection::Convert itsConverter;
- mutable std::mutex itsMutex;
+ mutable casacore::MeasFrame itsFrame;
+ mutable casacore::MDirection::Convert itsConverter;
+ mutable std::mutex itsMutex;
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/LofarMetaDataUtil.h b/LofarMetaDataUtil.h
index 8d042fc33e2c2c3ff9a0216476a5e4c91666f530..3975e90d06de0743b4aeb98b3da353bc598e0d39 100644
--- a/LofarMetaDataUtil.h
+++ b/LofarMetaDataUtil.h
@@ -37,47 +37,44 @@
namespace everybeam {
-const ElementResponseModel defaultElementResponseModel = ElementResponseModel::Hamaker;
+const ElementResponseModel defaultElementResponseModel =
+ ElementResponseModel::Hamaker;
/**
* @brief Read single station from MeasurementSet
- *
+ *
* @param ms Measurement set
* @param id Station id
* @param model Element response model
- * @return Station::Ptr
+ * @return Station::Ptr
*/
Station::Ptr readStation(
- const casacore::MeasurementSet &ms,
- unsigned int id,
+ const casacore::MeasurementSet &ms, unsigned int id,
const ElementResponseModel model = defaultElementResponseModel);
/**
* @brief Read multiple stations from measurment set into buffer out_it
* Loops over readStation for all the antennas in MeasurementSet
- *
+ *
* @tparam T Template type
* @param ms Measurement set
- * @param out_it Out buffer
+ * @param out_it Out buffer
* @param model Element Response buffer
*/
template <typename T>
void readStations(
- const casacore::MeasurementSet &ms,
- T out_it,
- const ElementResponseModel model = defaultElementResponseModel)
-{
- casacore::ROMSAntennaColumns antenna(ms.antenna());
- for(unsigned int i = 0; i < antenna.nrow(); ++i)
- {
- *out_it++ = readStation(ms, i, model);
- }
+ const casacore::MeasurementSet &ms, T out_it,
+ const ElementResponseModel model = defaultElementResponseModel) {
+ casacore::ROMSAntennaColumns antenna(ms.antenna());
+ for (unsigned int i = 0; i < antenna.nrow(); ++i) {
+ *out_it++ = readStation(ms, i, model);
+ }
}
// 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 everybeam
+} // namespace everybeam
#endif
diff --git a/MathUtil.h b/MathUtil.h
index 0a35cd0b0bc208758da98869a8a7c5cd518feb5a..077e01b466b23d7a762b18a713d811cdec36fa75 100644
--- a/MathUtil.h
+++ b/MathUtil.h
@@ -30,132 +30,112 @@
namespace everybeam {
-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 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 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 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 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 = {{M_PI_2 - atan2(cart[2], r), atan2(cart[1],
- cart[0])}};
- return thetaphi;
+inline vector2r_t cart2thetaphi(const vector3r_t &cart) {
+ real_t r = sqrt(cart[0] * cart[0] + cart[1] * cart[1]);
+ vector2r_t thetaphi = {{M_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;
+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;
+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 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 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 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 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 = {{ {{}}, {{}} }};
+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[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];
- }
+ 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;
+ return response;
}
-inline diag22c_t normalize(const raw_array_factor_t &raw)
-{
- diag22c_t af = {{}};
+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[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];
- }
+ if (raw.weight[1] != 0.0) {
+ af[1] = raw.factor[1] / raw.weight[1];
+ }
- return af;
+ return af;
}
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/MutablePtr.h b/MutablePtr.h
index b4ad46166c2b195ad0583d1aa4296381525767af..c8affd893f3d33d11fffc61b391875fec1fc687a 100644
--- a/MutablePtr.h
+++ b/MutablePtr.h
@@ -73,14 +73,15 @@ namespace everybeam {
*
* \endcode
*/
-template<typename T>
+template <typename T>
class MutablePtr : public std::shared_ptr<std::shared_ptr<T>> {
-public:
- MutablePtr(std::shared_ptr<T> ptr) : std::shared_ptr<std::shared_ptr<T>>(new std::shared_ptr<T>(ptr)) {}
- T& operator*() const { return **(this->get()); }
- std::shared_ptr<T> operator->() const { return *(this->get()); }
- void set(std::shared_ptr<T> ptr) { *(this->get()) = ptr;}
- explicit operator bool() const noexcept {return **this;}
+ public:
+ MutablePtr(std::shared_ptr<T> ptr)
+ : std::shared_ptr<std::shared_ptr<T>>(new std::shared_ptr<T>(ptr)) {}
+ T& operator*() const { return **(this->get()); }
+ std::shared_ptr<T> operator->() const { return *(this->get()); }
+ void set(std::shared_ptr<T> ptr) { *(this->get()) = ptr; }
+ explicit operator bool() const noexcept { return **this; }
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/Singleton.h b/Singleton.h
index 5d8006290dbd9878f280eeeb7645e2cdabe208bb..e0b80efaad65945fac3bce708fd2befa572f7af8 100644
--- a/Singleton.h
+++ b/Singleton.h
@@ -1,21 +1,19 @@
namespace everybeam {
-template<typename T>
-class Singleton
-{
- public:
- static std::shared_ptr<T> getInstance()
- {
- static std::shared_ptr<T> instance(new T()); // Guaranteed to be destroyed.
- // Instantiated on first use.
- return instance;
- }
- private:
- Singleton() {} // Constructor? (the {} brackets) are needed here.
+template <typename T>
+class Singleton {
+ public:
+ static std::shared_ptr<T> getInstance() {
+ static std::shared_ptr<T> instance(new T()); // Guaranteed to be destroyed.
+ // Instantiated on first use.
+ return instance;
+ }
- public:
- Singleton(Singleton const&) = delete;
- void operator=(Singleton const&) = delete;
+ private:
+ Singleton() {} // Constructor? (the {} brackets) are needed here.
+ public:
+ Singleton(Singleton const&) = delete;
+ void operator=(Singleton const&) = delete;
};
-} // namespace everybeam
+} // namespace everybeam
diff --git a/Station.h b/Station.h
index 9ea54f7ef78b916d32ac9ef555be12cb52bd6517..c0b9bb26ec9a6f484688bd51b42d44d124552fb2 100644
--- a/Station.h
+++ b/Station.h
@@ -36,385 +36,379 @@
#include <vector>
namespace everybeam {
-class Station
-{
-public:
- typedef std::shared_ptr<Station> Ptr;
- typedef std::shared_ptr<const Station> ConstPtr;
-// typedef std::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 std::string &name,
- const vector3r_t &position,
- const ElementResponseModel model);
-
- void setModel(const ElementResponseModel model);
-
- //! Return the name of the station.
- const std::string &name() const;
-
- //! 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);
-
- //! 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);
-
- //! 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).
- * \param rotate Boolean deciding if paralactic rotation should be applied.
- * \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 bool rotate = true)
- 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).
- * \param rotate Boolean deciding if paralactic rotation should be applied.
- * \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 rotate Boolean deciding if paralactic rotation should be applied.
- * \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 bool rotate=true) 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 rotate Boolean deciding if paralactic rotation should be applied.
- * \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 rotate Boolean deciding if paralactic rotation should be applied.
- * \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 bool rotate=true) 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 rotate Boolean deciding if paralactic rotation should be applied.
- * \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;
-
- // @}
-
- // ===================================================================
- // New methods introduced in refactor
- // ==================================================================
-
- const ElementResponse::Ptr get_element_response() {return itsElementResponse;}
-
- /**
- * @brief Compute the Jones matrix for the element response
- *
- * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
- * @param freq Frequency of the plane wave (Hz).
- * @param direction Direction of arrival (ITRF, m).
- * @param id Element id
- * @param rotate Boolean deciding if paralactic rotation should be applied.
- * @return matrix22c_t Jones matrix of element response
- */
- matrix22c_t elementResponse(real_t time, real_t freq,
- const vector3r_t &direction, size_t id, const bool rotate) const;
-
- /**
- * @brief Compute the Jones matrix for the element response
- *
- * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
- * @param freq Frequency of the plane wave (Hz).
- * @param direction Direction of arrival (ITRF, m).
- * @param rotate Boolean deciding if paralactic rotation should be applied.
- * @return matrix22c_t Jones matrix of element response
- */
- matrix22c_t elementResponse(real_t time, real_t freq,
- const vector3r_t &direction, const bool rotate = true) const;
-
- //! Specialized implementation of response function.
- matrix22c_t response(
- real_t time,
- real_t freq,
- const vector3r_t &direction) const
- {
- return itsAntenna->response(time, freq, direction);
- }
-
- //! Set antenna attribute, usually a BeamFormer, but can also be an Element
- void set_antenna(Antenna::Ptr antenna) {itsAntenna = antenna;}
-
- //! Set Element attribute
- void set_element(Element::Ptr element) {itsElement = element;}
-
-private:
-
- vector3r_t ncp(real_t time) const;
- vector3r_t ncppol0(real_t time) const;
- //! Compute the parallactic rotation.
- matrix22r_t rotation(real_t time, const vector3r_t &direction) const;
-
- std::string itsName;
- vector3r_t itsPosition;
- vector3r_t itsPhaseReference;
- ElementResponseModel itsElementResponseModel = ElementResponseModel::Unknown;
- ElementResponse::Ptr itsElementResponse;
- Element::Ptr itsElement;
-
- Antenna::Ptr itsAntenna;
-
- ITRFDirection::Ptr itsNCP;
- /** Reference direction for NCP observations.
- *
- * NCP pol0 is the direction used as reference in the coordinate system
- * when the target direction is close to/at the NCP. The regular coordinate
- * system rotates local east to that defined with respect to the NCP,
- * which is undefined at the NCP.
- * It is currently defined as ITRF position (1.0, 0.0, 0.0).
- *
- * Added by Maaijke Mevius, December 2018.
- */
- ITRFDirection::Ptr itsNCPPol0;
-
-
+class Station {
+ public:
+ typedef std::shared_ptr<Station> Ptr;
+ typedef std::shared_ptr<const Station> ConstPtr;
+ // typedef std::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 std::string &name, const vector3r_t &position,
+ const ElementResponseModel model);
+
+ void setModel(const ElementResponseModel model);
+
+ //! Return the name of the station.
+ const std::string &name() const;
+
+ //! 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);
+
+ //! 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);
+
+ //! 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).
+ * \param rotate Boolean deciding if paralactic rotation should be applied.
+ * \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 bool rotate = true) 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).
+ * \param rotate Boolean deciding if paralactic rotation should be applied.
+ * \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 rotate Boolean deciding if paralactic rotation should be applied.
+ * \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 bool rotate = true) 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 rotate Boolean deciding if paralactic rotation should be applied.
+ * \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 rotate Boolean deciding if paralactic rotation should be applied.
+ * \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 bool rotate = true) 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 rotate Boolean deciding if paralactic rotation should be applied.
+ * \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;
+
+ // @}
+
+ // ===================================================================
+ // New methods introduced in refactor
+ // ==================================================================
+
+ const ElementResponse::Ptr get_element_response() {
+ return itsElementResponse;
+ }
+
+ /**
+ * @brief Compute the Jones matrix for the element response
+ *
+ * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * @param freq Frequency of the plane wave (Hz).
+ * @param direction Direction of arrival (ITRF, m).
+ * @param id Element id
+ * @param rotate Boolean deciding if paralactic rotation should be applied.
+ * @return matrix22c_t Jones matrix of element response
+ */
+ matrix22c_t elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction, size_t id,
+ const bool rotate) const;
+
+ /**
+ * @brief Compute the Jones matrix for the element response
+ *
+ * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
+ * @param freq Frequency of the plane wave (Hz).
+ * @param direction Direction of arrival (ITRF, m).
+ * @param rotate Boolean deciding if paralactic rotation should be applied.
+ * @return matrix22c_t Jones matrix of element response
+ */
+ matrix22c_t elementResponse(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const bool rotate = true) const;
+
+ //! Specialized implementation of response function.
+ matrix22c_t response(real_t time, real_t freq,
+ const vector3r_t &direction) const {
+ return itsAntenna->response(time, freq, direction);
+ }
+
+ //! Set antenna attribute, usually a BeamFormer, but can also be an Element
+ void set_antenna(Antenna::Ptr antenna) { itsAntenna = antenna; }
+
+ //! Set Element attribute
+ void set_element(Element::Ptr element) { itsElement = element; }
+
+ private:
+ vector3r_t ncp(real_t time) const;
+ vector3r_t ncppol0(real_t time) const;
+ //! Compute the parallactic rotation.
+ matrix22r_t rotation(real_t time, const vector3r_t &direction) const;
+
+ std::string itsName;
+ vector3r_t itsPosition;
+ vector3r_t itsPhaseReference;
+ ElementResponseModel itsElementResponseModel = ElementResponseModel::Unknown;
+ ElementResponse::Ptr itsElementResponse;
+ Element::Ptr itsElement;
+
+ Antenna::Ptr itsAntenna;
+
+ ITRFDirection::Ptr itsNCP;
+ /** Reference direction for NCP observations.
+ *
+ * NCP pol0 is the direction used as reference in the coordinate system
+ * when the target direction is close to/at the NCP. The regular coordinate
+ * system rotates local east to that defined with respect to the NCP,
+ * which is undefined at the NCP.
+ * It is currently defined as ITRF position (1.0, 0.0, 0.0).
+ *
+ * Added by Maaijke Mevius, December 2018.
+ */
+ ITRFDirection::Ptr itsNCPPol0;
};
-
// ------------------------------------------------------------------------- //
// - 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 bool rotate) const
-{
- for(unsigned int i = 0; i < count; ++i)
- {
- *buffer++ = response(time, *freq++, direction, freq0, station0, tile0,
- rotate);
- }
+ const vector3r_t &direction, real_t freq0,
+ const vector3r_t &station0, const vector3r_t &tile0,
+ U buffer, const bool rotate) const {
+ for (unsigned int i = 0; i < count; ++i) {
+ *buffer++ =
+ response(time, *freq++, direction, freq0, station0, tile0, rotate);
+ }
}
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);
- }
+ 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 bool rotate) const
-{
- for(unsigned int i = 0; i < count; ++i)
- {
- *buffer++ = response(time, *freq++, direction, *freq0++, station0,
- tile0, rotate);
- }
+ const vector3r_t &direction, T freq0,
+ const vector3r_t &station0, const vector3r_t &tile0,
+ U buffer, const bool rotate) const {
+ for (unsigned int i = 0; i < count; ++i) {
+ *buffer++ =
+ response(time, *freq++, direction, *freq0++, station0, tile0, rotate);
+ }
}
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);
- }
+ 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 everybeam
+} // namespace everybeam
#endif
diff --git a/Types.h b/Types.h
index 21fd299577918716648dad5b7ccdee5a2cbd86c9..0385ca1e16a3d58a692dbbcde5c3b35b524b0b18 100644
--- a/Types.h
+++ b/Types.h
@@ -38,39 +38,38 @@ template <typename T, size_t N>
std::ostream &operator<<(std::ostream &out, const std::array<T, N> &obj);
/** Type used for real scalars. */
-typedef double real_t;
+typedef double real_t;
/** Type used for complex scalars. */
-typedef std::complex<double> complex_t;
+typedef std::complex<double> complex_t;
/** Type used for 2-dimensional real vectors. */
-typedef std::array<real_t, 2> vector2r_t;
+typedef std::array<real_t, 2> vector2r_t;
/** Type used for 3-dimensional real vectors. */
-typedef std::array<real_t, 3> vector3r_t;
+typedef std::array<real_t, 3> vector3r_t;
/** Type used for 2x2 real diagonal matrices. */
-typedef std::array<real_t, 2> diag22r_t;
+typedef std::array<real_t, 2> diag22r_t;
/** Type used for 2x2 complex diagonal matrices. */
-typedef std::array<complex_t, 2> diag22c_t;
+typedef std::array<complex_t, 2> diag22c_t;
/** Type used for 2x2 real matrices. */
-typedef std::array<std::array<real_t, 2>, 2> matrix22r_t;
+typedef std::array<std::array<real_t, 2>, 2> matrix22r_t;
/** Type used for 2x2 complex matrices. */
typedef std::array<std::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;
+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
@@ -83,24 +82,22 @@ struct raw_response_t
* 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;
+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;
};
template <typename T, size_t N>
-std::ostream &operator<<(std::ostream &out, const std::array<T, N> &obj)
-{
+std::ostream &operator<<(std::ostream &out, const std::array<T, N> &obj) {
print(out, obj.begin(), obj.end());
return out;
}
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/hamaker/HamakerCoeff.h b/hamaker/HamakerCoeff.h
index 1f46e1985916f0bd104e5ec6a78de6ffd0f17d82..0f15aaefa7d8c6df66fafa5308533285682bbe71 100644
--- a/hamaker/HamakerCoeff.h
+++ b/hamaker/HamakerCoeff.h
@@ -12,92 +12,78 @@
//! Hamaker coefficients
class HamakerCoefficients {
- public:
- //! Default constructor
- HamakerCoefficients();
-
- //! Constructor for reading coeff from file
- HamakerCoefficients(
- std::string& filename);
-
- //! Constructor for writing coeff to file
- HamakerCoefficients(
- const double freq_center,
- const double freq_range,
- const unsigned int nHarmonics,
- const unsigned int nPowerTheta,
- const unsigned int nPowerFreq);
-
- /**
- * @brief Set Hamaker coefficients
- *
- * @param n
- * @param t
- * @param f
- * @param value
- */
- void set_coeff(
- const unsigned int n,
- const unsigned int t,
- const unsigned int f,
- std::pair<std::complex<double>, std::complex<double>> value);
-
- void set_coeffs(
- const std::complex<double>* coeff);
-
- void set_coeffs(
- const std::vector<std::complex<double>> coeff);
-
- // Get
- size_t get_nr_coeffs() const;
-
- double get_freq_center() const { return m_freq_center; }
-
- double get_freq_range() const { return m_freq_range; }
-
- unsigned int get_nHarmonics() const { return m_nHarmonics; }
-
- unsigned int get_nPowerTheta() const { return m_nPowerTheta; }
-
- unsigned int get_nPowerFreq() const { return m_nPowerFreq; }
-
-
- std::pair<std::complex<double>, std::complex<double>> get_coeff(
- const unsigned int n,
- const unsigned int t,
- const unsigned int f);
-
- // HDF5 I/O
- void read_coeffs(
- std::string& filename);
-
- void write_coeffs(
- std::string& filename);
-
- // Debugging
- void print_coeffs();
-
- private:
- // Methods
- size_t get_index(
- const unsigned int n,
- const unsigned int t,
- const unsigned int f);
-
- // Parameters
- double m_freq_center;
- double m_freq_range;
- unsigned int m_nHarmonics;
- unsigned int m_nPowerTheta;
- unsigned int m_nPowerFreq;
- const unsigned int m_nInner = 2;
-
- // Data
- std::vector<std::complex<double>> m_coeff;
-
- // HDF5
- std::string m_dataset_name = "coeff";
- const unsigned int m_dataset_rank = 4;
+ public:
+ //! Default constructor
+ HamakerCoefficients();
+
+ //! Constructor for reading coeff from file
+ HamakerCoefficients(std::string& filename);
+
+ //! Constructor for writing coeff to file
+ HamakerCoefficients(const double freq_center, const double freq_range,
+ const unsigned int nHarmonics,
+ const unsigned int nPowerTheta,
+ const unsigned int nPowerFreq);
+
+ /**
+ * @brief Set Hamaker coefficients
+ *
+ * @param n
+ * @param t
+ * @param f
+ * @param value
+ */
+ void set_coeff(const unsigned int n, const unsigned int t,
+ const unsigned int f,
+ std::pair<std::complex<double>, std::complex<double>> value);
+
+ void set_coeffs(const std::complex<double>* coeff);
+
+ void set_coeffs(const std::vector<std::complex<double>> coeff);
+
+ // Get
+ size_t get_nr_coeffs() const;
+
+ double get_freq_center() const { return m_freq_center; }
+
+ double get_freq_range() const { return m_freq_range; }
+
+ unsigned int get_nHarmonics() const { return m_nHarmonics; }
+
+ unsigned int get_nPowerTheta() const { return m_nPowerTheta; }
+
+ unsigned int get_nPowerFreq() const { return m_nPowerFreq; }
+
+ std::pair<std::complex<double>, std::complex<double>> get_coeff(
+ const unsigned int n, const unsigned int t, const unsigned int f);
+
+ // HDF5 I/O
+ void read_coeffs(std::string& filename);
+
+ void write_coeffs(std::string& filename);
+
+ // Debugging
+ void print_coeffs();
+
+ private:
+ // Methods
+ size_t get_index(const unsigned int n, const unsigned int t,
+ const unsigned int f);
+
+ // Parameters
+ double m_freq_center;
+ double m_freq_range;
+ unsigned int m_nHarmonics;
+ unsigned int m_nPowerTheta;
+ unsigned int m_nPowerFreq;
+ const unsigned int m_nInner = 2;
+
+ // Data
+ std::vector<std::complex<double>> m_coeff;
+
+ // HDF5
+ std::string m_dataset_name = "coeff";
+ const unsigned int m_dataset_rank = 4;
};
#endif
diff --git a/hamaker/HamakerElementResponse.h b/hamaker/HamakerElementResponse.h
index 2a22c4eed721722509b657f58cf1c0437ce855f6..97edaab5d610552aca05f6c92979940a04330028 100644
--- a/hamaker/HamakerElementResponse.h
+++ b/hamaker/HamakerElementResponse.h
@@ -9,35 +9,31 @@
namespace everybeam {
//! Implementation of the Hamaker response model
-class HamakerElementResponse : public ElementResponse
-{
-public:
- virtual void response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const final override;
+class HamakerElementResponse : public ElementResponse {
+ public:
+ virtual void response(
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const final override;
- static std::shared_ptr<HamakerElementResponse> getInstance(const std::string &name);
+ static std::shared_ptr<HamakerElementResponse> getInstance(
+ const std::string &name);
-protected:
- std::string get_path(const char*) const;
+ protected:
+ std::string get_path(const char *) const;
- std::unique_ptr<HamakerCoefficients> m_coeffs;
+ std::unique_ptr<HamakerCoefficients> m_coeffs;
};
-class HamakerElementResponseHBA : public HamakerElementResponse
-{
-public:
- HamakerElementResponseHBA();
+class HamakerElementResponseHBA : public HamakerElementResponse {
+ public:
+ HamakerElementResponseHBA();
};
-class HamakerElementResponseLBA : public HamakerElementResponse
-{
-public:
- HamakerElementResponseLBA();
+class HamakerElementResponseLBA : public HamakerElementResponse {
+ public:
+ HamakerElementResponseLBA();
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/lobes/ElementResponse.h b/lobes/ElementResponse.h
index 90e65c18304557a864f258e99dbb935b42ddcb1a..1a9a0ae1e5103712382bce7620fb4e64de1a34b3 100644
--- a/lobes/ElementResponse.h
+++ b/lobes/ElementResponse.h
@@ -30,8 +30,7 @@
#include <complex>
-namespace everybeam
-{
+namespace everybeam {
// \addtogroup ElementResponse
// @{
@@ -49,7 +48,7 @@ namespace everybeam
// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
//
void element_response_lba(double freq, double theta, double phi,
- std::complex<double> (&response)[2][2]);
+ std::complex<double> (&response)[2][2]);
// Compute the response of an idealized LOFAR HBA dual dipole antenna to
// radiation at frequency freq (Hz) arriving from the direction given by theta,
@@ -64,7 +63,7 @@ void element_response_lba(double freq, double theta, double phi,
// phi: Azimuth in rad in the range [0.0, 2.0 * PI].
//
void element_response_hba(double freq, double theta, double phi,
- std::complex<double> (&response)[2][2]);
+ std::complex<double> (&response)[2][2]);
// Compute the response of an idealized LOFAR dual dipole antenna to radiation
// at frequency freq (Hz) arriving from the direction given by theta, phi (rad).
@@ -90,12 +89,13 @@ void element_response_hba(double freq, double theta, double phi,
// coeff_shape: Shape of the coefficient array, all dimensions should be > 0.
//
void element_response(double freq, double theta, double phi,
- std::complex<double> (&response)[2][2], double freq_center,
- double freq_range, const unsigned int (&coeff_shape)[3],
- const std::complex<double> coeff[]);
+ std::complex<double> (&response)[2][2],
+ double freq_center, double freq_range,
+ const unsigned int (&coeff_shape)[3],
+ const std::complex<double> coeff[]);
// @}
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/lobes/LOBESElementResponse.h b/lobes/LOBESElementResponse.h
index 738d6c21a2552b266133562b3392bb8f1603d9e0..39bae8884786b1b6086d46e4da1f925d8d28d484 100644
--- a/lobes/LOBESElementResponse.h
+++ b/lobes/LOBESElementResponse.h
@@ -8,21 +8,17 @@
namespace everybeam {
//! Implementation of the Lobes response model
-class LOBESElementResponse : public ElementResponse
-{
-public:
+class LOBESElementResponse : public ElementResponse {
+ public:
+ LOBESElementResponse(std::string name);
- LOBESElementResponse(std::string name);
+ virtual void response(
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const final override;
- virtual void response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const final override;
-
- static std::shared_ptr<LOBESElementResponse> getInstance(std::string name);
+ static std::shared_ptr<LOBESElementResponse> getInstance(std::string name);
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/lobes/lobes.h b/lobes/lobes.h
index 00112ee246ea5cd42badbd47e707b365e664e74a..d16241526dd7db27ce3f2282893ceba5848c9109 100644
--- a/lobes/lobes.h
+++ b/lobes/lobes.h
@@ -4,39 +4,39 @@
namespace py = pybind11;
-//! (Virtual) Beam model class
+//! (Virtual) Beam model class
class BeamModel {
+ public:
+ BeamModel() {}
-public:
-
- BeamModel() {}
-
- virtual py::array_t<std::complex<double>> eval(py::EigenDRef<const Eigen::ArrayXd> theta, py::EigenDRef<const Eigen::ArrayXd> phi)=0;
+ virtual py::array_t<std::complex<double>> eval(
+ py::EigenDRef<const Eigen::ArrayXd> theta,
+ py::EigenDRef<const Eigen::ArrayXd> phi) = 0;
};
//! Lobes beam model, wrapped with pybind11
class LobesBeamModel : public BeamModel {
+ public:
+ LobesBeamModel(const std::string &data_file_name);
-public:
-
- LobesBeamModel(const std::string &data_file_name);
+ py::array_t<std::complex<double>> eval(
+ py::EigenDRef<const Eigen::ArrayXd> theta,
+ py::EigenDRef<const Eigen::ArrayXd> phi) override;
- py::array_t<std::complex<double>> eval(py::EigenDRef<const Eigen::ArrayXd> theta, py::EigenDRef<const Eigen::ArrayXd> phi) override;
+ private:
+ Eigen::ArrayXXcd m_coefficients;
+ std::vector<size_t> m_coefficients_shape;
-private:
- Eigen::ArrayXXcd m_coefficients;
- std::vector<size_t> m_coefficients_shape;
-
- std::vector<double> m_frequencies;
- py::array_t<int> m_nms;
+ std::vector<double> m_frequencies;
+ py::array_t<int> m_nms;
};
//! Lobes beam model, not implemented!
class HamakerBeamModel : public BeamModel {
+ public:
+ HamakerBeamModel() {}
-public:
-
- HamakerBeamModel() {}
-
- py::array_t<std::complex<double>> eval(py::EigenDRef<const Eigen::ArrayXd> theta, py::EigenDRef<const Eigen::ArrayXd> phi) override {}
+ py::array_t<std::complex<double>> eval(
+ py::EigenDRef<const Eigen::ArrayXd> theta,
+ py::EigenDRef<const Eigen::ArrayXd> phi) override {}
};
diff --git a/oskar/OSKARDatafile.h b/oskar/OSKARDatafile.h
index 46de070e86d43716c821eda0280585a54f4be65c..a969423117957b677ca3a7d678736877f3436a3c 100644
--- a/oskar/OSKARDatafile.h
+++ b/oskar/OSKARDatafile.h
@@ -12,21 +12,19 @@
//! Oskar datafile interface
class Datafile {
- public:
- Datafile(
- const std::string& filename);
+ public:
+ Datafile(const std::string& filename);
- std::shared_ptr<Dataset> get(
- const unsigned int freq);
+ std::shared_ptr<Dataset> get(const unsigned int freq);
- private:
- // Coeffs;
- std::map<unsigned int, std::shared_ptr<Dataset>> m_map;
+ private:
+ // Coeffs;
+ std::map<unsigned int, std::shared_ptr<Dataset>> m_map;
- // HDF5
- std::string m_filename;
- std::unique_ptr<H5::H5File> m_h5_file;
- mutable std::mutex m_mutex;
+ // HDF5
+ std::string m_filename;
+ std::unique_ptr<H5::H5File> m_h5_file;
+ mutable std::mutex m_mutex;
};
#endif
\ No newline at end of file
diff --git a/oskar/OSKARDataset.h b/oskar/OSKARDataset.h
index 614167f02d3298134bd393b72699e37326fc916c..45a8779621024022e9f914c74d9d8d360ddbfad0 100644
--- a/oskar/OSKARDataset.h
+++ b/oskar/OSKARDataset.h
@@ -8,39 +8,34 @@
//! OSKAR dataset
class Dataset {
- public:
-
- /**
- * @brief Construct a new Dataset object given a h5 file and a
- * frequency
- *
- * @param h5_file H5 file (.h5)
- * @param freq Frequency to look for (Hz)
- */
- Dataset(
- H5::H5File& h5_file,
- const unsigned int freq);
-
- // Get
- size_t get_nr_elements() const { return m_nr_elements; };
- size_t get_l_max() const { return m_l_max; };
-
- std::complex<double>* get_alpha_ptr(
- const unsigned int element);
-
- private:
- // Methods
- size_t get_index(
- const unsigned int element) const;
-
- // Constants
- const unsigned int m_dataset_rank = 3;
-
- // Members
- std::vector<std::complex<double>> m_data;
- unsigned int m_nr_elements;
- unsigned int m_nr_coeffs;
- unsigned int m_l_max;
+ public:
+ /**
+ * @brief Construct a new Dataset object given a h5 file and a
+ * frequency
+ *
+ * @param h5_file H5 file (.h5)
+ * @param freq Frequency to look for (Hz)
+ */
+ Dataset(H5::H5File& h5_file, const unsigned int freq);
+
+ // Get
+ size_t get_nr_elements() const { return m_nr_elements; };
+ size_t get_l_max() const { return m_l_max; };
+
+ std::complex<double>* get_alpha_ptr(const unsigned int element);
+
+ private:
+ // Methods
+ size_t get_index(const unsigned int element) const;
+
+ // Constants
+ const unsigned int m_dataset_rank = 3;
+
+ // Members
+ std::vector<std::complex<double>> m_data;
+ unsigned int m_nr_elements;
+ unsigned int m_nr_coeffs;
+ unsigned int m_l_max;
};
#endif
\ No newline at end of file
diff --git a/oskar/OSKARElementResponse.h b/oskar/OSKARElementResponse.h
index ee8b772dd8e7f921aaaa4c45cce2604dfeb42918..4ff3f20ce0f79ad0cd0aac2ba5fa612b93dfe113 100644
--- a/oskar/OSKARElementResponse.h
+++ b/oskar/OSKARElementResponse.h
@@ -11,51 +11,39 @@
namespace everybeam {
//! Implementation of the OSKAR dipole response model
-class OSKARElementResponseDipole : public ElementResponse
-{
-public:
- static std::shared_ptr<OSKARElementResponseDipole> getInstance()
- {
- return Singleton<OSKARElementResponseDipole>::getInstance();
- }
-
- virtual void response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const final override;
-
+class OSKARElementResponseDipole : public ElementResponse {
+ public:
+ static std::shared_ptr<OSKARElementResponseDipole> getInstance() {
+ return Singleton<OSKARElementResponseDipole>::getInstance();
+ }
+
+ virtual void response(
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const final override;
};
//! Implementation of the OSKAR spherical wave response model
-class OSKARElementResponseSphericalWave : public ElementResponse
-{
-public:
- static std::shared_ptr<OSKARElementResponseSphericalWave> getInstance()
- {
- return Singleton<OSKARElementResponseSphericalWave>::getInstance();
- }
-
- OSKARElementResponseSphericalWave();
-
- virtual void response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const final override;
-
- virtual void response(
- int element_id,
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const final override;
-
-protected:
- std::string get_path(const char*) const;
-
- std::unique_ptr<Datafile> m_datafile;
+class OSKARElementResponseSphericalWave : public ElementResponse {
+ public:
+ static std::shared_ptr<OSKARElementResponseSphericalWave> getInstance() {
+ return Singleton<OSKARElementResponseSphericalWave>::getInstance();
+ }
+
+ OSKARElementResponseSphericalWave();
+
+ virtual void response(
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const final override;
+
+ virtual void response(
+ int element_id, double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const final override;
+
+ protected:
+ std::string get_path(const char*) const;
+
+ std::unique_ptr<Datafile> m_datafile;
};
-} // namespace everybeam
+} // namespace everybeam
#endif
diff --git a/oskar/oskar.h b/oskar/oskar.h
index 5b359d0ae9f4f2b90533dbdbde7cbe52899fb217..4488412ecf729f9139395ade9255d9517dfb21af 100644
--- a/oskar/oskar.h
+++ b/oskar/oskar.h
@@ -1,20 +1,14 @@
#include "oskar_vector_types.h"
#include "oskar_helper.h"
-
-void oskar_evaluate_dipole_pattern_double(
- const int num_points,
- const double* theta,
- const double* phi,
- const double freq_hz,
- const double dipole_length_m,
- std::complex<double>* pattern);
+void oskar_evaluate_dipole_pattern_double(const int num_points,
+ const double* theta,
+ const double* phi,
+ const double freq_hz,
+ const double dipole_length_m,
+ std::complex<double>* pattern);
void oskar_evaluate_spherical_wave_sum_double(
- const int num_points,
- const double* theta,
- const double* phi_x,
- const double* phi_y,
- const int l_max,
- const std::complex<double>* alpha,
+ const int num_points, const double* theta, const double* phi_x,
+ const double* phi_y, const int l_max, const std::complex<double>* alpha,
std::complex<double>* pattern);
diff --git a/oskar/oskar_helper.h b/oskar/oskar_helper.h
index 52b663c4266f0380e1046bf4e1f4644f142e8585..739d2a736bd527e039d9321d2be0e303e442d71b 100644
--- a/oskar/oskar_helper.h
+++ b/oskar/oskar_helper.h
@@ -27,92 +27,75 @@
*/
/* OUT_S += A * B */
-template<typename FP>
-inline void oskar_mul_add_complex(FP& out_s, FP A, FP B)
-{
- out_s.x += A.x * B.x; out_s.x -= A.y * B.y;
- out_s.y += A.x * B.y; out_s.y += A.y * B.x;
+template <typename FP>
+inline void oskar_mul_add_complex(FP& out_s, FP A, FP B) {
+ out_s.x += A.x * B.x;
+ out_s.x -= A.y * B.y;
+ out_s.y += A.x * B.y;
+ out_s.y += A.y * B.x;
}
/* OUT_S -= A * B */
-template<typename FP>
-inline void oskar_mul_sub_complex(FP& out_s, FP A, FP B)
-{
- out_s.x -= A.x * B.x; out_s.x += A.y * B.y;
- out_s.y -= A.x * B.y; out_s.y -= A.y * B.x;
+template <typename FP>
+inline void oskar_mul_sub_complex(FP& out_s, FP A, FP B) {
+ out_s.x -= A.x * B.x;
+ out_s.x += A.y * B.y;
+ out_s.y -= A.x * B.y;
+ out_s.y -= A.y * B.x;
}
-template<typename FP>
-inline void make_zero2(FP& X)
-{
- X.x = X.y = 0;
+template <typename FP>
+inline void make_zero2(FP& X) {
+ X.x = X.y = 0;
};
-template<typename FP, typename FP2>
-inline void oskar_sph_wave(
- FP pds,
- FP dpms,
- FP sin_p,
- FP cos_p,
- int M,
- FP2 a_tm,
- FP2 a_te,
- FP2& c_theta,
- FP2& c_phi)
-{
- FP2 qq, dd;
- qq.x = -cos_p * dpms;
- qq.y = -sin_p * dpms;
- dd.x = -sin_p * pds * (M);
- dd.y = cos_p * pds * (M);
- oskar_mul_add_complex(c_phi, qq, a_tm);
- oskar_mul_sub_complex(c_phi, dd, a_te);
- oskar_mul_add_complex(c_theta, dd, a_tm);
- oskar_mul_add_complex(c_theta, qq, a_te);
+template <typename FP, typename FP2>
+inline void oskar_sph_wave(FP pds, FP dpms, FP sin_p, FP cos_p, int M, FP2 a_tm,
+ FP2 a_te, FP2& c_theta, FP2& c_phi) {
+ FP2 qq, dd;
+ qq.x = -cos_p * dpms;
+ qq.y = -sin_p * dpms;
+ dd.x = -sin_p * pds * (M);
+ dd.y = cos_p * pds * (M);
+ oskar_mul_add_complex(c_phi, qq, a_tm);
+ oskar_mul_sub_complex(c_phi, dd, a_te);
+ oskar_mul_add_complex(c_theta, dd, a_tm);
+ oskar_mul_add_complex(c_theta, qq, a_te);
}
/* OUT0 is P_l^m(cos_theta),
* OUT1 is P_l^m(cos_theta) / sin_theta,
* OUT2 is d/d(cos_theta){P_l^m(cos_theta)} * sin_theta. */
-template<typename FP>
-inline void oskar_legendre2(
- int l,
- int m,
- FP cos_theta,
- FP sin_theta,
- FP& out0,
- FP& out1,
- FP& out2)
-{
- FP p0_ = (FP) 1, p1_;
- if (m > 0) {
- FP fact = (FP) 1;
- for (int i_ = 1; i_ <= m; ++i_) {
- p0_ *= (-fact) * sin_theta;
- fact += (FP) 2;
- }
+template <typename FP>
+inline void oskar_legendre2(int l, int m, FP cos_theta, FP sin_theta, FP& out0,
+ FP& out1, FP& out2) {
+ FP p0_ = (FP)1, p1_;
+ if (m > 0) {
+ FP fact = (FP)1;
+ for (int i_ = 1; i_ <= m; ++i_) {
+ p0_ *= (-fact) * sin_theta;
+ fact += (FP)2;
}
- out0 = cos_theta * (2 * m + 1) * p0_;
- if (l == m) {
- p1_ = out0;
- out0 = p0_;
- }
- else {
- p1_ = out0;
- for (int i_ = m + 2; i_ <= l + 1; ++i_) {
- out0 = p1_;
- p1_ = ((2 * i_ - 1) * cos_theta * out0 - (i_ + m - 1) * p0_) / (i_ - m);
- p0_ = out0;
- }
- }
- if (sin_theta != (FP) 0) {
- /* BOTH of these are divides. */
- out1 = out0 / sin_theta;
- out2 = (cos_theta * out0 * (l + 1) - p1_ * (l - m + 1)) / sin_theta;
- }
- else {
- out1 = out2 = (FP) 0;
+ }
+ out0 = cos_theta * (2 * m + 1) * p0_;
+ if (l == m) {
+ p1_ = out0;
+ out0 = p0_;
+ } else {
+ p1_ = out0;
+ for (int i_ = m + 2; i_ <= l + 1; ++i_) {
+ out0 = p1_;
+ p1_ = ((2 * i_ - 1) * cos_theta * out0 - (i_ + m - 1) * p0_) / (i_ - m);
+ p0_ = out0;
}
+ }
+ if (sin_theta != (FP)0) {
+ /* BOTH of these are divides. */
+ out1 = out0 / sin_theta;
+ out2 = (cos_theta * out0 * (l + 1) - p1_ * (l - m + 1)) / sin_theta;
+ } else {
+ out1 = out2 = (FP)0;
+ }
}
inline void oskar_sincos(float x, float* s, float* c) { sincosf(x, s, c); }
diff --git a/oskar/oskar_vector_types.h b/oskar/oskar_vector_types.h
index 461777e6a16bddf48fbe4cd6ccf17ba0fc7abe3f..0dde1ba37fa19bd91cadce803e0a08735f8a90cd 100644
--- a/oskar/oskar_vector_types.h
+++ b/oskar/oskar_vector_types.h
@@ -35,21 +35,21 @@
#ifdef __CUDACC__
/* Include the CUDA vector types header first, if we're compiling with nvcc. */
-# include <vector_types.h>
+#include <vector_types.h>
#endif
/* Memory alignment macros mirroring those used by CUDA. */
#if !(defined(__VECTOR_TYPES_H__) || defined(__CUDACC__))
-# if defined(__GNUC__)
-# define __align__(n) __attribute__((aligned(n)))
-# elif defined(_MSC_VER)
-# define __align__(n) __declspec(align(n))
-# endif
-# if defined(__GNUC__) || defined(_WIN64)
-# define __builtin_align__(a) __align__(a)
-# else
-# define __builtin_align__(a)
-# endif
+#if defined(__GNUC__)
+#define __align__(n) __attribute__((aligned(n)))
+#elif defined(_MSC_VER)
+#define __align__(n) __declspec(align(n))
+#endif
+#if defined(__GNUC__) || defined(_WIN64)
+#define __builtin_align__(a) __align__(a)
+#else
+#define __builtin_align__(a)
+#endif
/**
* @brief Two-element structure (single precision).
@@ -58,7 +58,9 @@
* Structure used to hold data for a length-2 vector.
* This must be compatible with the CUDA float2 type.
*/
-struct __builtin_align__(8) float2 { float x, y; };
+struct __builtin_align__(8) float2 {
+ float x, y;
+};
typedef struct float2 float2;
/**
@@ -68,7 +70,9 @@ typedef struct float2 float2;
* Structure used to hold data for a length-2 vector.
* This must be compatible with the CUDA double2 type.
*/
-struct __builtin_align__(16) double2 { double x, y; };
+struct __builtin_align__(16) double2 {
+ double x, y;
+};
typedef struct double2 double2;
#endif
@@ -82,7 +86,9 @@ typedef struct double2 double2;
* ( a b )
* ( c d )
*/
-struct __align__(32) float4c { float2 a, b, c, d; };
+struct __align__(32) float4c {
+ float2 a, b, c, d;
+};
typedef struct float4c float4c;
/**
@@ -95,7 +101,9 @@ typedef struct float4c float4c;
* ( a b )
* ( c d )
*/
-struct __align__(64) double4c { double2 a, b, c, d; };
+struct __align__(64) double4c {
+ double2 a, b, c, d;
+};
typedef struct double4c double4c;
#endif /* include guard */
diff --git a/scripts/run-clang-format.sh b/scripts/run-clang-format.sh
index cc985fcbfb53db8f3aacfa7281012fc2cdde1631..fd7998fd6cb413012a2a736fbe84ef670af6599a 100755
--- a/scripts/run-clang-format.sh
+++ b/scripts/run-clang-format.sh
@@ -12,4 +12,4 @@ cd $SCRIPT_PATH
cd ..
# Find all cpp headers ("*.h") and code files ("*.cc") source tree and execute clang-format. Exclude ./external director
-find . -path ./external -prune -o -name '*.h' -prune -o -name '*.cc' -exec clang-format -i -style=file \{\} +
+find . -path ./external -prune -o -type f \( -name "*.h" -o -name "*.cc" \) -exec clang-format -i -style=file \{\} +
diff --git a/test/beam-helper.h b/test/beam-helper.h
index 8a611629ec6e012a10c653d46863fbc94544eb8d..5566fa0e25593c9a94d591cea1cc349393ed74c6 100644
--- a/test/beam-helper.h
+++ b/test/beam-helper.h
@@ -16,77 +16,44 @@
using namespace everybeam;
-void GetPhaseCentreInfo(
- casacore::MeasurementSet& ms,
- size_t fieldId,
- double& ra,
- double& dec);
+void GetPhaseCentreInfo(casacore::MeasurementSet& ms, size_t fieldId,
+ double& ra, double& dec);
-void GetThetaPhiDirectionsZenith(
- vector2r_t* raDecDirections,
- size_t subgrid_size);
+void GetThetaPhiDirectionsZenith(vector2r_t* raDecDirections,
+ size_t subgrid_size);
-void GetITRFDirections(
- vector3r_t* itrfDirections,
- size_t subgrid_size,
- float image_size,
- double time,
- double phaseCentreRA,
- double phaseCentreDec);
+void GetITRFDirections(vector3r_t* itrfDirections, size_t subgrid_size,
+ float image_size, double time, double phaseCentreRA,
+ double phaseCentreDec);
-void StoreBeam(
- const std::string& filename,
- const std::complex<float>* buffer,
- size_t nStations,
- size_t width,
- size_t height);
+void StoreBeam(const std::string& filename, const std::complex<float>* buffer,
+ size_t nStations, size_t width, size_t height);
-void setITRFVector(
- const casacore::MDirection& itrfDir,
- vector3r_t& itrf);
+void setITRFVector(const casacore::MDirection& itrfDir, vector3r_t& itrf);
-inline matrix22c_t operator*(
- const matrix22c_t &arg0,
- const matrix22c_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 matrix22c_t operator*(const matrix22c_t& arg0, const matrix22c_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;
}
-template<typename T>
-void XYToLM(
- size_t x,
- size_t y,
- T pixelSizeX,
- T pixelSizeY,
- size_t width,
- size_t height,
- T &l,
- T &m)
-{
- T midX = (T) width / 2.0, midY = (T) height / 2.0;
- l = (midX - (T) x) * pixelSizeX;
- m = ((T) y - midY) * pixelSizeY;
+template <typename T>
+void XYToLM(size_t x, size_t y, T pixelSizeX, T pixelSizeY, size_t width,
+ size_t height, T& l, T& m) {
+ T midX = (T)width / 2.0, midY = (T)height / 2.0;
+ l = (midX - (T)x) * pixelSizeX;
+ m = ((T)y - midY) * pixelSizeY;
}
-void GetRaDecZenith(
- vector3r_t position,
- double time,
- double& ra,
- double& dec);
+void GetRaDecZenith(vector3r_t position, double time, double& ra, double& dec);
-std::string GetFieldName(
- casacore::MeasurementSet& ms,
- unsigned int field_id = 0);
+std::string GetFieldName(casacore::MeasurementSet& ms,
+ unsigned int field_id = 0);
-std::string GetStationName(
- casacore::MeasurementSet& ms,
- unsigned int station_id);
+std::string GetStationName(casacore::MeasurementSet& ms,
+ unsigned int station_id);
-unsigned int GetNrAntennas(
- casacore::MeasurementSet& ms,
- unsigned int field_id);
\ No newline at end of file
+unsigned int GetNrAntennas(casacore::MeasurementSet& ms, unsigned int field_id);
\ No newline at end of file
diff --git a/test/tElementBeamCommon.h b/test/tElementBeamCommon.h
index 4d49106fdfde0f12109861d5ce5c2b710ff39f15..8b64c91eadd2eb6a80cdf90bd7922a6e40a0eeab 100644
--- a/test/tElementBeamCommon.h
+++ b/test/tElementBeamCommon.h
@@ -4,144 +4,137 @@
#include "beam-helper.h"
-void calculateElementBeams(
- everybeam::Station::Ptr& station,
- std::vector<vector2r_t>& thetaPhiDirections,
- size_t nr_antennas,
- unsigned int subgrid_size,
- double frequency,
- std::vector<std::complex<float>>& buffer)
-{
- typedef std::complex<float> Data[nr_antennas][subgrid_size][subgrid_size][4];
- Data* data_ptr = (Data *) buffer.data();
-
- auto elementResponse = station->get_element_response();
-
- #pragma omp parallel for
- for (size_t a = 0; a < nr_antennas; a++) {
- for (unsigned y = 0; y < subgrid_size; y++) {
- for (unsigned x = 0; x < subgrid_size; x++) {
- // Get theta, phi
- auto direction_thetaphi = thetaPhiDirections[y * subgrid_size + x];
- double theta = direction_thetaphi[0];
- double phi = direction_thetaphi[1];
-
- // Compute gain
- std::complex<double> gainMatrix[2][2] = {0.0};
- if (std::isfinite(theta) && std::isfinite(phi)) {
- elementResponse->response(a, frequency, theta, phi, gainMatrix);
- }
-
- // Store gain
- std::complex<float>* antBufferPtr = (*data_ptr)[a][y][x];
- antBufferPtr[0] = gainMatrix[0][0];
- antBufferPtr[1] = gainMatrix[0][1];
- antBufferPtr[2] = gainMatrix[1][0];
- antBufferPtr[3] = gainMatrix[1][1];
- }
+void calculateElementBeams(everybeam::Station::Ptr& station,
+ std::vector<vector2r_t>& thetaPhiDirections,
+ size_t nr_antennas, unsigned int subgrid_size,
+ double frequency,
+ std::vector<std::complex<float>>& buffer) {
+ typedef std::complex<float> Data[nr_antennas][subgrid_size][subgrid_size][4];
+ Data* data_ptr = (Data*)buffer.data();
+
+ auto elementResponse = station->get_element_response();
+
+#pragma omp parallel for
+ for (size_t a = 0; a < nr_antennas; a++) {
+ for (unsigned y = 0; y < subgrid_size; y++) {
+ for (unsigned x = 0; x < subgrid_size; x++) {
+ // Get theta, phi
+ auto direction_thetaphi = thetaPhiDirections[y * subgrid_size + x];
+ double theta = direction_thetaphi[0];
+ double phi = direction_thetaphi[1];
+
+ // Compute gain
+ std::complex<double> gainMatrix[2][2] = {0.0};
+ if (std::isfinite(theta) && std::isfinite(phi)) {
+ elementResponse->response(a, frequency, theta, phi, gainMatrix);
}
+
+ // Store gain
+ std::complex<float>* antBufferPtr = (*data_ptr)[a][y][x];
+ antBufferPtr[0] = gainMatrix[0][0];
+ antBufferPtr[1] = gainMatrix[0][1];
+ antBufferPtr[2] = gainMatrix[1][0];
+ antBufferPtr[3] = gainMatrix[1][1];
+ }
}
+ }
}
-void calculateElementBeams(
- everybeam::Station::Ptr& station,
- std::vector<vector3r_t>& itrfDirections,
- size_t nr_antennas,
- unsigned int subgrid_size,
- double time,
- double frequency,
- std::vector<std::complex<float>>& buffer)
-{
- typedef std::complex<float> Data[nr_antennas][subgrid_size][subgrid_size][4];
- Data* data_ptr = (Data *) buffer.data();
-
- #pragma omp parallel for
- for (size_t a = 0; a < nr_antennas; a++) {
- for (unsigned y = 0; y < subgrid_size; y++) {
- for (unsigned x = 0; x < subgrid_size; x++) {
- // Get direction
- auto direction = itrfDirections[y * subgrid_size + x];
-
- // Compute gain
- matrix22c_t gainMatrix = { 0.0 };
- if (std::isfinite(direction[0])) {
- gainMatrix = station->elementResponse(time, frequency, direction, a, true);
- }
-
- // Store gain
- std::complex<float>* antBufferPtr = (*data_ptr)[a][y][x];
- antBufferPtr[0] = gainMatrix[0][0];
- antBufferPtr[1] = gainMatrix[0][1];
- antBufferPtr[2] = gainMatrix[1][0];
- antBufferPtr[3] = gainMatrix[1][1];
- }
+void calculateElementBeams(everybeam::Station::Ptr& station,
+ std::vector<vector3r_t>& itrfDirections,
+ size_t nr_antennas, unsigned int subgrid_size,
+ double time, double frequency,
+ std::vector<std::complex<float>>& buffer) {
+ typedef std::complex<float> Data[nr_antennas][subgrid_size][subgrid_size][4];
+ Data* data_ptr = (Data*)buffer.data();
+
+#pragma omp parallel for
+ for (size_t a = 0; a < nr_antennas; a++) {
+ for (unsigned y = 0; y < subgrid_size; y++) {
+ for (unsigned x = 0; x < subgrid_size; x++) {
+ // Get direction
+ auto direction = itrfDirections[y * subgrid_size + x];
+
+ // Compute gain
+ matrix22c_t gainMatrix = {0.0};
+ if (std::isfinite(direction[0])) {
+ gainMatrix =
+ station->elementResponse(time, frequency, direction, a, true);
}
+
+ // Store gain
+ std::complex<float>* antBufferPtr = (*data_ptr)[a][y][x];
+ antBufferPtr[0] = gainMatrix[0][0];
+ antBufferPtr[1] = gainMatrix[0][1];
+ antBufferPtr[2] = gainMatrix[1][0];
+ antBufferPtr[3] = gainMatrix[1][1];
+ }
}
+ }
}
-void run(
- everybeam::ElementResponseModel elementResponseModel,
- double frequency,
- std::string& input_filename,
- std::string& output_filename)
-{
- // Open measurement set
- std::cout << ">> Opening measurementset: " << input_filename << std::endl;
- casacore::MeasurementSet ms(input_filename);
-
- // Print frequency
- std::clog << "Frequency: " << frequency * 1e-6 << " Mhz" << std::endl;
-
- // Set number of stations to 1
- size_t nr_stations = 1;
-
- // Read number of timesteps
- size_t nr_timesteps = ms.nrow();
- std::clog << "Number of timesteps: " << nr_timesteps << std::endl;
-
-
- // Read observation time
- casacore::ScalarColumn<double> timeColumn(ms, ms.columnName(casacore::MSMainEnums::TIME));
- double currentTime = timeColumn(nr_timesteps / 2);
-
- // Read station
- size_t field_id = 0;
- size_t station_id = 0;
- auto station = readStation(ms, station_id, elementResponseModel);
- auto field_name = GetFieldName(ms, field_id);
- auto station_name = GetStationName(ms, station_id);
- auto nr_antennas = GetNrAntennas(ms, field_id);
- std::cout << "field: " << field_name << std::endl;
- std::cout << "station: " << station_name << std::endl;
- std::cout << "nr_antennas: " << nr_antennas << std::endl;
-
- // Compute RA and DEC of zenith at currentTime
- double zenithRA, zenithDec;
- GetRaDecZenith(station->position(), currentTime, zenithRA, zenithDec);
- std::clog << "RA: " << zenithRA << std::endl;
- std::clog << "DEC: " << zenithDec << std::endl;
-
- // Imaging parameters
- float image_size = M_PI; // in radians
- size_t subgrid_size = 32; // in pixels
-
- // Compute element beams from theta, phi
- std::cout << ">>> Computing element beams theta, phi" << std::endl;
- std::vector<vector2r_t> thetaPhiDirections(subgrid_size * subgrid_size);
- GetThetaPhiDirectionsZenith(thetaPhiDirections.data(), subgrid_size);
- std::vector<std::complex<float>> beam_thetaphi(subgrid_size*subgrid_size*4*nr_antennas);
- calculateElementBeams(station, thetaPhiDirections, nr_antennas, subgrid_size, frequency, beam_thetaphi);
-
- // Compute element beams from itrf coordinates
- // TODO: the Station::elementResponse method does not work properly
- #if 0
+void run(everybeam::ElementResponseModel elementResponseModel, double frequency,
+ std::string& input_filename, std::string& output_filename) {
+ // Open measurement set
+ std::cout << ">> Opening measurementset: " << input_filename << std::endl;
+ casacore::MeasurementSet ms(input_filename);
+
+ // Print frequency
+ std::clog << "Frequency: " << frequency * 1e-6 << " Mhz" << std::endl;
+
+ // Set number of stations to 1
+ size_t nr_stations = 1;
+
+ // Read number of timesteps
+ size_t nr_timesteps = ms.nrow();
+ std::clog << "Number of timesteps: " << nr_timesteps << std::endl;
+
+ // Read observation time
+ casacore::ScalarColumn<double> timeColumn(
+ ms, ms.columnName(casacore::MSMainEnums::TIME));
+ double currentTime = timeColumn(nr_timesteps / 2);
+
+ // Read station
+ size_t field_id = 0;
+ size_t station_id = 0;
+ auto station = readStation(ms, station_id, elementResponseModel);
+ auto field_name = GetFieldName(ms, field_id);
+ auto station_name = GetStationName(ms, station_id);
+ auto nr_antennas = GetNrAntennas(ms, field_id);
+ std::cout << "field: " << field_name << std::endl;
+ std::cout << "station: " << station_name << std::endl;
+ std::cout << "nr_antennas: " << nr_antennas << std::endl;
+
+ // Compute RA and DEC of zenith at currentTime
+ double zenithRA, zenithDec;
+ GetRaDecZenith(station->position(), currentTime, zenithRA, zenithDec);
+ std::clog << "RA: " << zenithRA << std::endl;
+ std::clog << "DEC: " << zenithDec << std::endl;
+
+ // Imaging parameters
+ float image_size = M_PI; // in radians
+ size_t subgrid_size = 32; // in pixels
+
+ // Compute element beams from theta, phi
+ std::cout << ">>> Computing element beams theta, phi" << std::endl;
+ std::vector<vector2r_t> thetaPhiDirections(subgrid_size * subgrid_size);
+ GetThetaPhiDirectionsZenith(thetaPhiDirections.data(), subgrid_size);
+ std::vector<std::complex<float>> beam_thetaphi(subgrid_size * subgrid_size *
+ 4 * nr_antennas);
+ calculateElementBeams(station, thetaPhiDirections, nr_antennas, subgrid_size,
+ frequency, beam_thetaphi);
+
+// Compute element beams from itrf coordinates
+// TODO: the Station::elementResponse method does not work properly
+#if 0
std::cout << ">>> Computing element beams itrfs" << std::endl;
std::vector<vector3r_t> itrfDirections(subgrid_size * subgrid_size);
GetITRFDirections(itrfDirections.data(), subgrid_size, image_size, currentTime, zenithRA, zenithDec);
std::vector<std::complex<float>> beam_itrf(subgrid_size*subgrid_size*4*nr_antennas);
calculateElementBeams(station, itrfDirections, nr_antennas, subgrid_size, currentTime, frequency, beam_itrf);
- #endif
+#endif
- // Store aterm
- StoreBeam(output_filename, beam_thetaphi.data(), nr_antennas, subgrid_size, subgrid_size);
+ // Store aterm
+ StoreBeam(output_filename, beam_thetaphi.data(), nr_antennas, subgrid_size,
+ subgrid_size);
}
diff --git a/test/tStationBeamCommon.h b/test/tStationBeamCommon.h
index 300599621279fc010d0ab24473621d58d509e7d8..c3d8e41c403f55a7a6ad137cc00fbd1c22639e7f 100644
--- a/test/tStationBeamCommon.h
+++ b/test/tStationBeamCommon.h
@@ -4,104 +4,104 @@
#include "beam-helper.h"
-void calculateStationBeams(
- std::vector<everybeam::Station::Ptr>& stations,
- std::vector<vector3r_t>& itrfDirections,
- vector3r_t stationDirection,
- vector3r_t tileDirection,
- unsigned int subgrid_size,
- std::vector<std::complex<float>>& buffer,
- double time, double frequency)
-{
- typedef std::complex<float> Data[stations.size()][subgrid_size][subgrid_size][4];
- Data* data_ptr = (Data *) buffer.data();
-
- #pragma omp parallel for
- for (size_t s = 0; s < stations.size(); s++) {
- for (unsigned y = 0; y < subgrid_size; y++) {
- for (unsigned x = 0; x < subgrid_size; x++) {
- auto direction = itrfDirections[y * subgrid_size + x];
- auto freq_beamformer = frequency;
- matrix22c_t gainMatrix =
- stations[s]->response(
- time, frequency, direction, freq_beamformer,
- stationDirection, tileDirection);
-
- std::complex<float>* antBufferPtr = (*data_ptr)[s][y][x];
-
- matrix22c_t stationGain = gainMatrix;
- antBufferPtr[0] = stationGain[0][0];
- antBufferPtr[1] = stationGain[0][1];
- antBufferPtr[2] = stationGain[1][0];
- antBufferPtr[3] = stationGain[1][1];
- }
- }
+void calculateStationBeams(std::vector<everybeam::Station::Ptr>& stations,
+ std::vector<vector3r_t>& itrfDirections,
+ vector3r_t stationDirection,
+ vector3r_t tileDirection, unsigned int subgrid_size,
+ std::vector<std::complex<float>>& buffer,
+ double time, double frequency) {
+ typedef std::complex<float> Data[stations.size()][subgrid_size][subgrid_size]
+ [4];
+ Data* data_ptr = (Data*)buffer.data();
+
+#pragma omp parallel for
+ for (size_t s = 0; s < stations.size(); s++) {
+ for (unsigned y = 0; y < subgrid_size; y++) {
+ for (unsigned x = 0; x < subgrid_size; x++) {
+ auto direction = itrfDirections[y * subgrid_size + x];
+ auto freq_beamformer = frequency;
+ matrix22c_t gainMatrix =
+ stations[s]->response(time, frequency, direction, freq_beamformer,
+ stationDirection, tileDirection);
+
+ std::complex<float>* antBufferPtr = (*data_ptr)[s][y][x];
+
+ matrix22c_t stationGain = gainMatrix;
+ antBufferPtr[0] = stationGain[0][0];
+ antBufferPtr[1] = stationGain[0][1];
+ antBufferPtr[2] = stationGain[1][0];
+ antBufferPtr[3] = stationGain[1][1];
+ }
}
+ }
}
-void run(
- everybeam::ElementResponseModel elementResponseModel,
- double frequency,
- std::string& input_filename,
- std::string& output_filename)
-{
- // Open measurement set
- std::cout << ">> Opening measurementset: " << input_filename << std::endl;
- casacore::MeasurementSet ms(input_filename);
-
- // Print frequency
- std::clog << "Frequency: " << frequency * 1e-6 << " Mhz" << std::endl;
-
- // Read number of stations
- size_t nr_stations = ms.antenna().nrow();
- std::clog << "Number of stations: " << nr_stations << std::endl;
-
- // Read number of timesteps
- size_t nr_timesteps = ms.nrow();
- std::clog << "Number of timesteps: " << nr_timesteps << std::endl;
-
- // Read field id
- casacore::ROScalarColumn<int> fieldIdColumn(ms, casacore::MS::columnName(casacore::MSMainEnums::FIELD_ID));
- size_t fieldId = fieldIdColumn.getColumn()[0];
- std::clog << "Field ID: " << fieldId << std::endl;
-
- // Read phase centre info
- double phaseCentreRA, phaseCentreDec;
- GetPhaseCentreInfo(ms, fieldId, phaseCentreRA, phaseCentreDec);
- std::clog << "RA: " << phaseCentreRA << std::endl;
- std::clog << "DEC: " << phaseCentreDec << std::endl;
-
- // Read observation time
- casacore::ScalarColumn<double> timeColumn(ms, ms.columnName(casacore::MSMainEnums::TIME));
- double currentTime = timeColumn(nr_timesteps / 2);
-
- // Read stations
- std::vector<everybeam::Station::Ptr> stations;
- stations.resize(nr_stations);
- readStations(ms, stations.begin(), elementResponseModel);
-
- // Imaging parameters
- float image_size = 0.5; // in radians
- size_t subgrid_size = 32; // in pixels
-
- // Evaluate beam directions in ITRF coordinates
- std::cout << ">>> Computing directions to evaluate beam" << std::endl;
- std::vector<vector3r_t> itrfDirections(subgrid_size * subgrid_size);
- GetITRFDirections(itrfDirections.data(), subgrid_size, image_size, currentTime, phaseCentreRA, phaseCentreDec);
-
- // Set station beam direction to centre of field
- vector3r_t stationDirection = itrfDirections[(subgrid_size/2) * subgrid_size + (subgrid_size/2)];
-
- // Set tile beam direction equal to station direction
- vector3r_t tileDirection = stationDirection;
-
- // Compute station beams
- std::cout << ">>> Computing station beams" << std::endl;
- std::vector<std::complex<float>> beams;
- beams.resize(subgrid_size*subgrid_size*4*nr_stations);
- calculateStationBeams(stations, itrfDirections, stationDirection, tileDirection, subgrid_size, beams, currentTime, frequency);
-
- // Store aterm
- std::cout << ">>> Writing beam images to: " << output_filename << std::endl;
- StoreBeam(output_filename, beams.data(), nr_stations, subgrid_size, subgrid_size);
+void run(everybeam::ElementResponseModel elementResponseModel, double frequency,
+ std::string& input_filename, std::string& output_filename) {
+ // Open measurement set
+ std::cout << ">> Opening measurementset: " << input_filename << std::endl;
+ casacore::MeasurementSet ms(input_filename);
+
+ // Print frequency
+ std::clog << "Frequency: " << frequency * 1e-6 << " Mhz" << std::endl;
+
+ // Read number of stations
+ size_t nr_stations = ms.antenna().nrow();
+ std::clog << "Number of stations: " << nr_stations << std::endl;
+
+ // Read number of timesteps
+ size_t nr_timesteps = ms.nrow();
+ std::clog << "Number of timesteps: " << nr_timesteps << std::endl;
+
+ // Read field id
+ casacore::ROScalarColumn<int> fieldIdColumn(
+ ms, casacore::MS::columnName(casacore::MSMainEnums::FIELD_ID));
+ size_t fieldId = fieldIdColumn.getColumn()[0];
+ std::clog << "Field ID: " << fieldId << std::endl;
+
+ // Read phase centre info
+ double phaseCentreRA, phaseCentreDec;
+ GetPhaseCentreInfo(ms, fieldId, phaseCentreRA, phaseCentreDec);
+ std::clog << "RA: " << phaseCentreRA << std::endl;
+ std::clog << "DEC: " << phaseCentreDec << std::endl;
+
+ // Read observation time
+ casacore::ScalarColumn<double> timeColumn(
+ ms, ms.columnName(casacore::MSMainEnums::TIME));
+ double currentTime = timeColumn(nr_timesteps / 2);
+
+ // Read stations
+ std::vector<everybeam::Station::Ptr> stations;
+ stations.resize(nr_stations);
+ readStations(ms, stations.begin(), elementResponseModel);
+
+ // Imaging parameters
+ float image_size = 0.5; // in radians
+ size_t subgrid_size = 32; // in pixels
+
+ // Evaluate beam directions in ITRF coordinates
+ std::cout << ">>> Computing directions to evaluate beam" << std::endl;
+ std::vector<vector3r_t> itrfDirections(subgrid_size * subgrid_size);
+ GetITRFDirections(itrfDirections.data(), subgrid_size, image_size,
+ currentTime, phaseCentreRA, phaseCentreDec);
+
+ // Set station beam direction to centre of field
+ vector3r_t stationDirection =
+ itrfDirections[(subgrid_size / 2) * subgrid_size + (subgrid_size / 2)];
+
+ // Set tile beam direction equal to station direction
+ vector3r_t tileDirection = stationDirection;
+
+ // Compute station beams
+ std::cout << ">>> Computing station beams" << std::endl;
+ std::vector<std::complex<float>> beams;
+ beams.resize(subgrid_size * subgrid_size * 4 * nr_stations);
+ calculateStationBeams(stations, itrfDirections, stationDirection,
+ tileDirection, subgrid_size, beams, currentTime,
+ frequency);
+
+ // Store aterm
+ std::cout << ">>> Writing beam images to: " << output_filename << std::endl;
+ StoreBeam(output_filename, beams.data(), nr_stations, subgrid_size,
+ subgrid_size);
}