diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000000000000000000000000000000000000..a0a17c17ab6507017f353333e987fd5c810a43df
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,48 @@
+---
+# BasedOnStyle: Google
+AccessModifierOffset: -1
+ConstructorInitializerIndentWidth: 4
+AlignEscapedNewlinesLeft: true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortIfStatementsOnASingleLine: true
+AllowShortLoopsOnASingleLine: true
+AlwaysBreakTemplateDeclarations: true
+AlwaysBreakBeforeMultilineStrings: true
+BreakBeforeBinaryOperators: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BinPackParameters: true
+ColumnLimit: 80
+ConstructorInitializerAllOnOneLineOrOnePerLine: true
+DerivePointerBinding: true
+ExperimentalAutoDetectBinPacking: false
+IndentCaseLabels: true
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCSpaceBeforeProtocolList: false
+PenaltyBreakBeforeFirstCallParameter: 1
+PenaltyBreakComment: 60
+PenaltyBreakString: 1000
+PenaltyBreakFirstLessLess: 120
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 200
+PointerBindsToType: true
+SpacesBeforeTrailingComments: 2
+Cpp11BracedListStyle: true
+Standard: Auto
+IndentWidth: 2
+TabWidth: 8
+UseTab: Never
+BreakBeforeBraces: Attach
+IndentFunctionDeclarationAfterType: true
+SortIncludes: false
+SpacesInParentheses: false
+SpacesInAngles: false
+SpaceInEmptyParentheses: false
+SpacesInCStyleCastParentheses: false
+SpaceAfterControlStatementKeyword: true
+SpaceBeforeAssignmentOperators: true
+ContinuationIndentWidth: 4
+...
+
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 8fcd2c9b697691b092100f5a3eacbc555d4c45e1..a35760fc9058ceba5a3d43fb5bbed792d16ead5f 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,5 +1,6 @@
stages:
- build-base
+ - clang-format
- build-lofarbeam
- build-doc
- deploy-doc
@@ -18,6 +19,13 @@ build-base:
-f ./test/Dockerfile-base .
- docker push $CI_REGISTRY_IMAGE:base
+clang-format:
+ stage: clang-format
+ image: $CI_REGISTRY_IMAGE:base
+ script:
+ - apt-get update && apt-get -y install clang-format
+ - ./scripts/clang-format-check.sh
+
build-lofarbeam:
stage: build-lofarbeam
script:
diff --git a/Antenna.cc b/Antenna.cc
index 1685705dc9734e2f11fd80bed91b48144d9fab81..c5f4b5bb702384d2dbc4fee0844d50a69f20f5c8 100644
--- a/Antenna.cc
+++ b/Antenna.cc
@@ -6,14 +6,14 @@ namespace LOFAR {
namespace StationResponse {
vector3r_t Antenna::transform_to_local_direction(const vector3r_t &direction) {
- vector3r_t local_direction {
- dot(m_coordinate_system.axes.p, direction),
- dot(m_coordinate_system.axes.q, direction),
- dot(m_coordinate_system.axes.r, direction),
- };
+ vector3r_t local_direction{
+ dot(m_coordinate_system.axes.p, direction),
+ dot(m_coordinate_system.axes.q, direction),
+ dot(m_coordinate_system.axes.r, direction),
+ };
- return local_direction;
+ return local_direction;
}
-} // namespace StationResponse
-} // namespace LOFAR
\ No newline at end of file
+} // namespace StationResponse
+} // namespace LOFAR
\ No newline at end of file
diff --git a/BeamFormer.cc b/BeamFormer.cc
index f9280888f20040e55133f861dc75e8869904e8cc..c0efd2fe3a4ad0010bf811cb3993e1de8800c8a5 100644
--- a/BeamFormer.cc
+++ b/BeamFormer.cc
@@ -5,109 +5,110 @@
#include <cmath>
-
namespace LOFAR {
namespace StationResponse {
vector3r_t BeamFormer::transform_to_local_position(const vector3r_t &position) {
- vector3r_t dposition {
- position[0] - m_coordinate_system.origin[0],
- position[1] - m_coordinate_system.origin[1],
- position[2] - m_coordinate_system.origin[2]
- };
- vector3r_t local_position {
- dot(m_coordinate_system.axes.p, dposition),
- dot(m_coordinate_system.axes.q, dposition),
- dot(m_coordinate_system.axes.r, dposition),
- };
-
- return local_position;
+ vector3r_t dposition{position[0] - m_coordinate_system.origin[0],
+ position[1] - m_coordinate_system.origin[1],
+ position[2] - m_coordinate_system.origin[2]};
+ vector3r_t local_position{
+ dot(m_coordinate_system.axes.p, dposition),
+ dot(m_coordinate_system.axes.q, dposition),
+ dot(m_coordinate_system.axes.r, dposition),
+ };
+
+ return local_position;
}
-
-std::vector<std::complex<double>> BeamFormer::compute_geometric_response(double freq, const vector3r_t &direction) const
-{
- std::vector<std::complex<double>> result;
- result.reserve(m_antennas.size());
- for (auto &antenna : m_antennas)
- {
-// std::cout << "(" << antenna->m_phase_reference_position[0] << ", " <<
-// antenna->m_phase_reference_position[1] << ", " <<
-// antenna->m_phase_reference_position[2] << ")" << std::endl;
-//
-// std::cout << "(" << m_local_phase_reference_position[0] << ", " <<
-// m_local_phase_reference_position[1] << ", " <<
-// m_local_phase_reference_position[2] << ")" << std::endl;
-//
-// std::cout << "(" << (antenna->m_phase_reference_position[0] - m_local_phase_reference_position[0]) << ", " <<
-// (antenna->m_phase_reference_position[1] - m_local_phase_reference_position[1]) << ", " <<
-// (antenna->m_phase_reference_position[2] - m_local_phase_reference_position[2]) << ")" << std::endl;
-//
-// std::cout << "======" << std::endl;
-
- double dl = direction[0] * (antenna->m_phase_reference_position[0] - m_local_phase_reference_position[0]) +
- direction[1] * (antenna->m_phase_reference_position[1] - m_local_phase_reference_position[1]) +
- direction[2] * (antenna->m_phase_reference_position[2] - m_local_phase_reference_position[2]);
-
- double phase = -2 * M_PI * dl / (Constants::c / freq);
- result.push_back({std::sin(phase), std::cos(phase)});
- }
- return result;
+std::vector<std::complex<double>> BeamFormer::compute_geometric_response(
+ double freq, const vector3r_t &direction) const {
+ std::vector<std::complex<double>> result;
+ result.reserve(m_antennas.size());
+ for (auto &antenna : m_antennas) {
+ // std::cout << "(" << antenna->m_phase_reference_position[0] << ",
+ // " <<
+ // antenna->m_phase_reference_position[1] << ", " <<
+ // antenna->m_phase_reference_position[2] << ")" << std::endl;
+ //
+ // std::cout << "(" << m_local_phase_reference_position[0] << ", "
+ // <<
+ // m_local_phase_reference_position[1] << ", " <<
+ // m_local_phase_reference_position[2] << ")" << std::endl;
+ //
+ // std::cout << "(" << (antenna->m_phase_reference_position[0] -
+ // m_local_phase_reference_position[0]) << ", " <<
+ // (antenna->m_phase_reference_position[1] -
+ // m_local_phase_reference_position[1]) << ", " <<
+ // (antenna->m_phase_reference_position[2] -
+ // m_local_phase_reference_position[2]) << ")" << std::endl;
+ //
+ // std::cout << "======" << std::endl;
+
+ double dl = direction[0] * (antenna->m_phase_reference_position[0] -
+ m_local_phase_reference_position[0]) +
+ direction[1] * (antenna->m_phase_reference_position[1] -
+ m_local_phase_reference_position[1]) +
+ direction[2] * (antenna->m_phase_reference_position[2] -
+ m_local_phase_reference_position[2]);
+
+ double phase = -2 * M_PI * dl / (Constants::c / freq);
+ result.push_back({std::sin(phase), std::cos(phase)});
+ }
+ return result;
}
-std::vector<std::pair<std::complex<double>,std::complex<double>>> BeamFormer::compute_weights(const vector3r_t &pointing, double freq) const
-{
- std::vector<std::pair<std::complex<double>,std::complex<double>>> result;
- double weight_sum[2] = {0.0, 0.0};
- auto geometric_response = compute_geometric_response(freq, pointing);
- result.reserve(geometric_response.size());
- for (unsigned int antenna_idx = 0; antenna_idx < m_antennas.size(); ++antenna_idx)
- {
- auto phasor = geometric_response[antenna_idx];
- result.push_back({
- std::conj(phasor) * (1.0 * m_antennas[antenna_idx]->m_enabled[0]),
- std::conj(phasor) * (1.0 * m_antennas[antenna_idx]->m_enabled[1])
- });
- weight_sum[0] += (1.0 * m_antennas[antenna_idx]->m_enabled[0]);
- weight_sum[1] += (1.0 * m_antennas[antenna_idx]->m_enabled[1]);
- }
- for (unsigned int antenna_idx = 0; antenna_idx < m_antennas.size(); ++antenna_idx)
- {
- result[antenna_idx].first /= weight_sum[0];
- result[antenna_idx].second /= weight_sum[1];
- }
-
- return result;
+std::vector<std::pair<std::complex<double>, std::complex<double>>>
+ BeamFormer::compute_weights(const vector3r_t &pointing, double freq) const {
+ std::vector<std::pair<std::complex<double>, std::complex<double>>> result;
+ double weight_sum[2] = {0.0, 0.0};
+ auto geometric_response = compute_geometric_response(freq, pointing);
+ result.reserve(geometric_response.size());
+ for (unsigned int antenna_idx = 0; antenna_idx < m_antennas.size();
+ ++antenna_idx) {
+ auto phasor = geometric_response[antenna_idx];
+ result.push_back(
+ {std::conj(phasor) * (1.0 * m_antennas[antenna_idx]->m_enabled[0]),
+ std::conj(phasor) * (1.0 * m_antennas[antenna_idx]->m_enabled[1])});
+ weight_sum[0] += (1.0 * m_antennas[antenna_idx]->m_enabled[0]);
+ weight_sum[1] += (1.0 * m_antennas[antenna_idx]->m_enabled[1]);
+ }
+ for (unsigned int antenna_idx = 0; antenna_idx < m_antennas.size();
+ ++antenna_idx) {
+ result[antenna_idx].first /= weight_sum[0];
+ result[antenna_idx].second /= weight_sum[1];
+ }
+
+ return result;
}
-
-matrix22c_t BeamFormer::local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const
-{
- auto weights = compute_weights(options.station0, options.freq0);
- auto geometric_response = compute_geometric_response(freq, direction);
-
- matrix22c_t result = {0};
- for (unsigned int antenna_idx = 0; antenna_idx < m_antennas.size(); ++antenna_idx) {
- auto antenna = m_antennas[antenna_idx];
- auto antenna_weight = weights[antenna_idx];
- auto antenna_geometric_reponse = geometric_response[antenna_idx];
-
- matrix22c_t antenna_response = antenna->response(time, freq, direction, options);
-
- result[0][0] += antenna_weight.first * antenna_geometric_reponse * antenna_response[0][0];
- result[0][1] += antenna_weight.first * antenna_geometric_reponse * antenna_response[0][1];
- result[1][0] += antenna_weight.second * antenna_geometric_reponse * antenna_response[1][0];
- result[1][1] += antenna_weight.second * antenna_geometric_reponse * antenna_response[1][1];
-
- }
- return result;
+matrix22c_t BeamFormer::local_response(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const Options &options) const {
+ auto weights = compute_weights(options.station0, options.freq0);
+ auto geometric_response = compute_geometric_response(freq, direction);
+
+ matrix22c_t result = {0};
+ for (unsigned int antenna_idx = 0; antenna_idx < m_antennas.size();
+ ++antenna_idx) {
+ auto antenna = m_antennas[antenna_idx];
+ auto antenna_weight = weights[antenna_idx];
+ auto antenna_geometric_reponse = geometric_response[antenna_idx];
+
+ matrix22c_t antenna_response =
+ antenna->response(time, freq, direction, options);
+
+ result[0][0] += antenna_weight.first * antenna_geometric_reponse *
+ antenna_response[0][0];
+ result[0][1] += antenna_weight.first * antenna_geometric_reponse *
+ antenna_response[0][1];
+ result[1][0] += antenna_weight.second * antenna_geometric_reponse *
+ antenna_response[1][0];
+ result[1][1] += antenna_weight.second * antenna_geometric_reponse *
+ antenna_response[1][1];
+ }
+ return result;
}
-
-} // namespace StationResponse
-} // namespace LOFAR
-
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/Element.cc b/Element.cc
index 91a6f3fc4cabf6eb5357f1dbb8d18c4f1749a475..608330836302e194a433554bab592db2ff9f3a97 100644
--- a/Element.cc
+++ b/Element.cc
@@ -4,39 +4,34 @@
namespace LOFAR {
namespace StationResponse {
-matrix22c_t Element::local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- size_t id,
- const Options &options) const
-{
- vector2r_t thetaphi = cart2thetaphi(direction);
+matrix22c_t Element::local_response(real_t time, real_t freq,
+ const vector3r_t &direction, size_t id,
+ const Options &options) const {
+ vector2r_t thetaphi = cart2thetaphi(direction);
- matrix22c_t result;
- static_assert(sizeof(std::complex<double>[2][2]) == sizeof(matrix22c_t));
- m_element_response->response(id, freq, thetaphi[0], thetaphi[1], reinterpret_cast<std::complex<double> (&)[2][2]>(result));
+ matrix22c_t result;
+ static_assert(sizeof(std::complex<double>[2][2]) == sizeof(matrix22c_t));
+ m_element_response->response(
+ id, freq, thetaphi[0], thetaphi[1],
+ reinterpret_cast<std::complex<double>(&)[2][2]>(result));
- if (options.rotate) {
- vector3r_t up = {0.0, 0.0, 1.0};
- vector3r_t e_phi = normalize(cross(up, direction));
- vector3r_t e_theta = cross(e_phi, direction);
- matrix22r_t rotation;
- rotation[0] = {dot(e_theta, options.north), dot(e_theta, options.east)};
- rotation[1] = {dot(e_phi, options.north), dot(e_phi, options.east)};
- result = result * rotation;
- }
- return result;
+ if (options.rotate) {
+ vector3r_t up = {0.0, 0.0, 1.0};
+ vector3r_t e_phi = normalize(cross(up, direction));
+ vector3r_t e_theta = cross(e_phi, direction);
+ matrix22r_t rotation;
+ rotation[0] = {dot(e_theta, options.north), dot(e_theta, options.east)};
+ rotation[1] = {dot(e_phi, options.north), dot(e_phi, options.east)};
+ result = result * rotation;
+ }
+ return result;
}
-matrix22c_t Element::local_response(
- real_t time,
- real_t freq,
- const vector3r_t &direction,
- const Options &options) const
-{
- return local_response(time, freq, direction, m_id, options);
+matrix22c_t Element::local_response(real_t time, real_t freq,
+ const vector3r_t &direction,
+ const Options &options) const {
+ return local_response(time, freq, direction, m_id, options);
}
-} // namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/ElementResponse.cc b/ElementResponse.cc
index 79809f269477497f56bde6c6888135b8df56cd7a..6b913d9aab4da243d40ce20a44860d5bf2cb2b07 100644
--- a/ElementResponse.cc
+++ b/ElementResponse.cc
@@ -3,19 +3,28 @@
namespace LOFAR {
namespace StationResponse {
-std::ostream& operator<<(std::ostream& os, ElementResponseModel model )
-{
- switch (model)
- {
- case Unknown: os << "Unknown"; break;
- case Hamaker: os << "Hamaker"; break;
- case LOBES : os << "LOBES"; break;
- case OSKARDipole : os << "OSKARDipole"; break;
- case OSKARSphericalWave : os << "OSKARSphericalWave"; break;
- default : os.setstate(std::ios_base::failbit);
- }
- return os;
+std::ostream& operator<<(std::ostream& os, ElementResponseModel model) {
+ switch (model) {
+ case Unknown:
+ os << "Unknown";
+ break;
+ case Hamaker:
+ os << "Hamaker";
+ break;
+ case LOBES:
+ os << "LOBES";
+ break;
+ case OSKARDipole:
+ os << "OSKARDipole";
+ break;
+ case OSKARSphericalWave:
+ os << "OSKARSphericalWave";
+ break;
+ default:
+ os.setstate(std::ios_base::failbit);
+ }
+ return os;
}
-} // namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/ITRFConverter.cc b/ITRFConverter.cc
index 16eb6d1bef5c44642e7cb8c7fd47c0e693e4793d..e55204f85386e55a6ddeacf6edfc2215f3aa8f5f 100644
--- a/ITRFConverter.cc
+++ b/ITRFConverter.cc
@@ -10,76 +10,75 @@
namespace LOFAR {
namespace StationResponse {
- //TODO: Initialize converter with a time (and fixed position) and convert specific directions.
- // Needed for wslean as well as for the makestationresponse executable.
-
-
-ITRFConverter::ITRFConverter(real_t time)
-{
- //create ITRF Direction from fixed stationposition
- casacore::MVPosition mvPosition(ITRFDirection::LOFARPosition()[0], ITRFDirection::LOFARPosition()[1], ITRFDirection::LOFARPosition()[2]);
- casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
- casacore::MEpoch timeEpoch(casacore::Quantity(time, "s"));
- itsFrame = casacore::MeasFrame(timeEpoch, mPosition);
-
- // Order of angles seems to be longitude (along the equator), lattitude
- // (towards the pole).
- itsConverter = casacore::MDirection::Convert(
- casacore::MDirection::J2000,
- casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
+// TODO: Initialize converter with a time (and fixed position) and convert
+// specific directions.
+// Needed for wslean as well as for the makestationresponse executable.
+
+ITRFConverter::ITRFConverter(real_t time) {
+ // create ITRF Direction from fixed stationposition
+ casacore::MVPosition mvPosition(ITRFDirection::LOFARPosition()[0],
+ ITRFDirection::LOFARPosition()[1],
+ ITRFDirection::LOFARPosition()[2]);
+ casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
+ casacore::MEpoch timeEpoch(casacore::Quantity(time, "s"));
+ itsFrame = casacore::MeasFrame(timeEpoch, mPosition);
+
+ // Order of angles seems to be longitude (along the equator), lattitude
+ // (towards the pole).
+ itsConverter = casacore::MDirection::Convert(
+ casacore::MDirection::J2000,
+ casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
}
-void ITRFConverter::setTime(real_t time)
-{
- // Cannot use MeasFrame::resetEpoch(Double), because that assumes the
- // argument is UTC in (fractional) days (MJD).
- itsFrame.resetEpoch(casacore::Quantity(time, "s"));
+void ITRFConverter::setTime(real_t time) {
+ // Cannot use MeasFrame::resetEpoch(Double), because that assumes the
+ // argument is UTC in (fractional) days (MJD).
+ itsFrame.resetEpoch(casacore::Quantity(time, "s"));
}
-vector3r_t ITRFConverter::j2000ToITRF(const vector2r_t &j2000Direction) const
-{
- casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1]);
- casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
- const casacore::MVDirection mvITRF = itsConverter(mDirection).getValue();
+vector3r_t ITRFConverter::j2000ToITRF(const vector2r_t &j2000Direction) const {
+ casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+ const casacore::MVDirection mvITRF = itsConverter(mDirection).getValue();
- return vector3r_t{{mvITRF(0), mvITRF(1), mvITRF(2)}};
+ return vector3r_t{{mvITRF(0), mvITRF(1), mvITRF(2)}};
}
-vector3r_t ITRFConverter::j2000ToITRF(const vector3r_t &j2000Direction) const
-{
- casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1], j2000Direction[2]);
- casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+vector3r_t ITRFConverter::j2000ToITRF(const vector3r_t &j2000Direction) const {
+ casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1],
+ j2000Direction[2]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
- const casacore::MVDirection mvITRF = itsConverter(mDirection).getValue();
+ const casacore::MVDirection mvITRF = itsConverter(mDirection).getValue();
- return vector3r_t{{mvITRF(0), mvITRF(1), mvITRF(2)}};
+ return vector3r_t{{mvITRF(0), mvITRF(1), mvITRF(2)}};
}
-vector3r_t ITRFConverter::toITRF(const casacore::MDirection &direction) const{
- const casacore::MVDirection mvITRF = itsConverter(direction).getValue();
- return vector3r_t{{mvITRF(0), mvITRF(1), mvITRF(2)}};
+vector3r_t ITRFConverter::toITRF(const casacore::MDirection &direction) const {
+ const casacore::MVDirection mvITRF = itsConverter(direction).getValue();
+ return vector3r_t{{mvITRF(0), mvITRF(1), mvITRF(2)}};
}
-casacore::MDirection ITRFConverter::toDirection(const vector2r_t &j2000Direction) const
-{
- casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1]);
- casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
- return itsConverter(mDirection);
+casacore::MDirection ITRFConverter::toDirection(
+ const vector2r_t &j2000Direction) const {
+ casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+ return itsConverter(mDirection);
}
-casacore::MDirection ITRFConverter::toDirection(const vector3r_t &j2000Direction) const
-{
- casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1], j2000Direction[2]);
- casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+casacore::MDirection ITRFConverter::toDirection(
+ const vector3r_t &j2000Direction) const {
+ casacore::MVDirection mvDirection(j2000Direction[0], j2000Direction[1],
+ j2000Direction[2]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
- return itsConverter(mDirection);
+ return itsConverter(mDirection);
}
-casacore::MDirection ITRFConverter::toDirection(const casacore::MDirection &direction) const
-{
- return itsConverter(direction);
+casacore::MDirection ITRFConverter::toDirection(
+ const casacore::MDirection &direction) const {
+ return itsConverter(direction);
}
-
-} //# namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/ITRFDirection.cc b/ITRFDirection.cc
index 1f484926320826e469c952930f194e00e344cefd..7a0b264c1e575073ca4f63ea57295b32afb29348 100644
--- a/ITRFDirection.cc
+++ b/ITRFDirection.cc
@@ -29,70 +29,67 @@
namespace LOFAR {
namespace StationResponse {
+// ITRF position of CS002LBA, just to use a fixed reference
+const vector3r_t ITRFDirection::itsLOFARPosition = {
+ {826577.022720000, 461022.995082000, 5064892.814}};
- //ITRF position of CS002LBA, just to use a fixed reference
- const vector3r_t ITRFDirection::itsLOFARPosition = {{ 826577.022720000,461022.995082000,5064892.814 }};
-
- //TODO: initialize converter with a time (and fixed position) and convert specific directions. Needed for wslean as well as for the makestationresponse executable.
-
+// TODO: initialize converter with a time (and fixed position) and convert
+// specific directions. Needed for wslean as well as for the makestationresponse
+// executable.
ITRFDirection::ITRFDirection(const vector3r_t &position,
- const vector2r_t &direction)
-{
- casacore::MVPosition mvPosition(position[0], position[1], position[2]);
- casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
- itsFrame = casacore::MeasFrame(casacore::MEpoch(), mPosition);
-
- // Order of angles seems to be longitude (along the equator), lattitude
- // (towards the pole).
- casacore::MVDirection mvDirection(direction[0], direction[1]);
- casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
- itsConverter = casacore::MDirection::Convert(mDirection,
- casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
+ const vector2r_t &direction) {
+ casacore::MVPosition mvPosition(position[0], position[1], position[2]);
+ casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
+ itsFrame = casacore::MeasFrame(casacore::MEpoch(), mPosition);
+
+ // Order of angles seems to be longitude (along the equator), lattitude
+ // (towards the pole).
+ casacore::MVDirection mvDirection(direction[0], direction[1]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+ itsConverter = casacore::MDirection::Convert(
+ mDirection,
+ casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
}
-ITRFDirection::ITRFDirection(const vector2r_t &direction):
- ITRFDirection(itsLOFARPosition, direction)
-{
- //create ITRF Direction from fixed stationposition
+ITRFDirection::ITRFDirection(const vector2r_t &direction)
+ : ITRFDirection(itsLOFARPosition, direction) {
+ // create ITRF Direction from fixed stationposition
}
ITRFDirection::ITRFDirection(const vector3r_t &position,
- const vector3r_t &direction)
-{
- casacore::MVPosition mvPosition(position[0], position[1], position[2]);
- casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
- itsFrame = casacore::MeasFrame(casacore::MEpoch(), mPosition);
-
- casacore::MVDirection mvDirection(direction[0], direction[1], direction[2]);
- casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
- itsConverter = casacore::MDirection::Convert(mDirection,
- casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
+ const vector3r_t &direction) {
+ casacore::MVPosition mvPosition(position[0], position[1], position[2]);
+ casacore::MPosition mPosition(mvPosition, casacore::MPosition::ITRF);
+ itsFrame = casacore::MeasFrame(casacore::MEpoch(), mPosition);
+
+ casacore::MVDirection mvDirection(direction[0], direction[1], direction[2]);
+ casacore::MDirection mDirection(mvDirection, casacore::MDirection::J2000);
+ itsConverter = casacore::MDirection::Convert(
+ mDirection,
+ casacore::MDirection::Ref(casacore::MDirection::ITRF, itsFrame));
}
-ITRFDirection::ITRFDirection(const vector3r_t &direction):
- ITRFDirection(itsLOFARPosition, direction)
+ITRFDirection::ITRFDirection(const vector3r_t &direction)
+ : ITRFDirection(itsLOFARPosition, direction)
{
- //create ITRF Direction from fixed stationposition
+ // create ITRF Direction from fixed stationposition
}
+vector3r_t ITRFDirection::at(real_t time) const {
+ std::lock_guard<std::mutex> lock(itsMutex);
-vector3r_t ITRFDirection::at(real_t time) const
-{
- std::lock_guard<std::mutex> lock(itsMutex);
-
- // Cannot use MeasFrame::resetEpoch(Double), because that assumes the
- // argument is UTC in (fractional) days (MJD).
- itsFrame.resetEpoch(casacore::Quantity(time, "s"));
+ // Cannot use MeasFrame::resetEpoch(Double), because that assumes the
+ // argument is UTC in (fractional) days (MJD).
+ itsFrame.resetEpoch(casacore::Quantity(time, "s"));
- const casacore::MDirection &mITRF = itsConverter();
- const casacore::MVDirection &mvITRF = mITRF.getValue();
+ const casacore::MDirection &mITRF = itsConverter();
+ const casacore::MVDirection &mvITRF = mITRF.getValue();
- vector3r_t itrf = {{mvITRF(0), mvITRF(1), mvITRF(2)}};
- return itrf;
+ vector3r_t itrf = {{mvITRF(0), mvITRF(1), mvITRF(2)}};
+ return itrf;
}
-
-} //# namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/LofarMetaDataUtil.cc b/LofarMetaDataUtil.cc
index 3572d4353a77b2914e5054a489c60e21709e7705..6a7f6f9d3d7829215496bbb4c584b7b9233bcd01 100644
--- a/LofarMetaDataUtil.cc
+++ b/LofarMetaDataUtil.cc
@@ -58,8 +58,16 @@ namespace LOFAR {
namespace StationResponse {
constexpr Antenna::CoordinateSystem::Axes lofar_antenna_orientation = {
- {-std::sqrt(.5), -std::sqrt(.5), 0.0,},
- { std::sqrt(.5), -std::sqrt(.5), 0.0,},
+ {
+ -std::sqrt(.5),
+ -std::sqrt(.5),
+ 0.0,
+ },
+ {
+ std::sqrt(.5),
+ -std::sqrt(.5),
+ 0.0,
+ },
{0.0, 0.0, 1.0},
};
@@ -69,64 +77,53 @@ constexpr Antenna::CoordinateSystem::Axes lofar_antenna_orientation = {
// {0.0, 0.0, 1.0},
// };
-
using namespace casacore;
-typedef std::array<vector3r_t, 16> TileConfig;
+typedef std::array<vector3r_t, 16> TileConfig;
-bool hasColumn(const Table &table, const string &column)
-{
- return table.tableDesc().isColumn(column);
+bool hasColumn(const Table &table, const string &column) {
+ return table.tableDesc().isColumn(column);
}
-bool hasSubTable(const Table &table, const string &name)
-{
- return table.keywordSet().isDefined(name);
+bool hasSubTable(const Table &table, const string &name) {
+ return table.keywordSet().isDefined(name);
}
-Table getSubTable(const Table &table, const string &name)
-{
- return table.keywordSet().asTable(name);
+Table getSubTable(const Table &table, const string &name) {
+ return table.keywordSet().asTable(name);
}
-TileConfig readTileConfig(const Table &table, unsigned int row)
-{
- ROArrayQuantColumn<Double> c_tile_offset(table, "TILE_ELEMENT_OFFSET", "m");
+TileConfig readTileConfig(const Table &table, unsigned int row) {
+ ROArrayQuantColumn<Double> c_tile_offset(table, "TILE_ELEMENT_OFFSET", "m");
- // Read tile configuration for HBA antenna fields.
- Matrix<Quantity> aips_offset = c_tile_offset(row);
- assert(aips_offset.ncolumn() == TileAntenna::TileConfig::size());
+ // Read tile configuration for HBA antenna fields.
+ Matrix<Quantity> aips_offset = c_tile_offset(row);
+ assert(aips_offset.ncolumn() == TileAntenna::TileConfig::size());
- TileConfig config;
- for(unsigned int i = 0; i < config.size(); ++i)
- {
- config[i][0] = aips_offset(0, i).getValue();
- config[i][1] = aips_offset(1, i).getValue();
- config[i][2] = aips_offset(2, i).getValue();
- }
- return config;
+ TileConfig config;
+ for (unsigned int i = 0; i < config.size(); ++i) {
+ config[i][0] = aips_offset(0, i).getValue();
+ config[i][1] = aips_offset(1, i).getValue();
+ config[i][2] = aips_offset(2, i).getValue();
+ }
+ return config;
}
void transformToFieldCoordinates(TileConfig &config,
- const Antenna::CoordinateSystem::Axes &axes)
-{
- for(unsigned int i = 0; i < config.size(); ++i)
- {
- const vector3r_t position = config[i];
- config[i][0] = dot(position, axes.p);
- config[i][1] = dot(position, axes.q);
- config[i][2] = dot(position, axes.r);
- }
+ const Antenna::CoordinateSystem::Axes &axes) {
+ for (unsigned int i = 0; i < config.size(); ++i) {
+ const vector3r_t position = config[i];
+ config[i][0] = dot(position, axes.p);
+ config[i][1] = dot(position, axes.q);
+ config[i][2] = dot(position, axes.r);
+ }
}
-vector3r_t transformToFieldCoordinates(const vector3r_t &position,
- const Antenna::CoordinateSystem::Axes &axes)
-{
- const vector3r_t result{
- dot(position, axes.p),
- dot(position, axes.q),
- dot(position, axes.r)};
- return result;
+vector3r_t transformToFieldCoordinates(
+ const vector3r_t &position, const Antenna::CoordinateSystem::Axes &axes) {
+ const vector3r_t result{dot(position, axes.p), dot(position, axes.q),
+ dot(position, axes.r)};
+ return result;
}
// AntennaField::CoordinateSystem readCoordinateSystemAartfaac(
@@ -156,237 +153,233 @@ vector3r_t transformToFieldCoordinates(const vector3r_t &position,
// }
Antenna::CoordinateSystem readCoordinateSystem(const Table &table,
- unsigned int id)
-{
- ROArrayQuantColumn<Double> c_position(table, "POSITION", "m");
- ROArrayQuantColumn<Double> c_axes(table, "COORDINATE_AXES", "m");
-
- // Read antenna field center (ITRF).
- Vector<Quantity> aips_position = c_position(id);
- assert(aips_position.size() == 3);
-
- vector3r_t position = {{aips_position(0).getValue(),
- aips_position(1).getValue(), aips_position(2).getValue()}};
-
- // Read antenna field coordinate axes (ITRF).
- Matrix<Quantity> aips_axes = c_axes(id);
- assert(aips_axes.shape().isEqual(IPosition(2, 3, 3)));
-
- vector3r_t p = {{aips_axes(0, 0).getValue(), aips_axes(1, 0).getValue(),
- aips_axes(2, 0).getValue()}};
- vector3r_t q = {{aips_axes(0, 1).getValue(), aips_axes(1, 1).getValue(),
- aips_axes(2, 1).getValue()}};
- vector3r_t r = {{aips_axes(0, 2).getValue(), aips_axes(1, 2).getValue(),
- aips_axes(2, 2).getValue()}};
-
- Antenna::CoordinateSystem coordinate_system = {position, {p, q, r}};
- return coordinate_system;
+ unsigned int id) {
+ ROArrayQuantColumn<Double> c_position(table, "POSITION", "m");
+ ROArrayQuantColumn<Double> c_axes(table, "COORDINATE_AXES", "m");
+
+ // Read antenna field center (ITRF).
+ Vector<Quantity> aips_position = c_position(id);
+ assert(aips_position.size() == 3);
+
+ vector3r_t position = {{aips_position(0).getValue(),
+ aips_position(1).getValue(),
+ aips_position(2).getValue()}};
+
+ // Read antenna field coordinate axes (ITRF).
+ Matrix<Quantity> aips_axes = c_axes(id);
+ assert(aips_axes.shape().isEqual(IPosition(2, 3, 3)));
+
+ vector3r_t p = {{aips_axes(0, 0).getValue(), aips_axes(1, 0).getValue(),
+ aips_axes(2, 0).getValue()}};
+ vector3r_t q = {{aips_axes(0, 1).getValue(), aips_axes(1, 1).getValue(),
+ aips_axes(2, 1).getValue()}};
+ vector3r_t r = {{aips_axes(0, 2).getValue(), aips_axes(1, 2).getValue(),
+ aips_axes(2, 2).getValue()}};
+
+ Antenna::CoordinateSystem coordinate_system = {position, {p, q, r}};
+ return coordinate_system;
}
-BeamFormer::Ptr make_tile(std::string name, Antenna::CoordinateSystem coordinate_system,
- TileConfig tile_config, ElementResponse::Ptr)
-{
-}
-
-BeamFormer::Ptr readAntennaField(const Table &table, unsigned int id, ElementResponse::Ptr element_response)
-{
- Antenna::CoordinateSystem coordinate_system = readCoordinateSystem(table, id);
-// std::cout << "coordinate_system: " << std::endl;
-// std::cout << " axes.p: " << coordinate_system.axes.p[0] << ", " << coordinate_system.axes.p[1] << ", " << coordinate_system.axes.p[2] << std::endl;
-// std::cout << " axes.q: " << coordinate_system.axes.q[0] << ", " << coordinate_system.axes.q[1] << ", " << coordinate_system.axes.q[2] << std::endl;
-// std::cout << " axes.r: " << coordinate_system.axes.r[0] << ", " << coordinate_system.axes.r[1] << ", " << coordinate_system.axes.r[2] << std::endl;
- BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system));
-
- ROScalarColumn<String> c_name(table, "NAME");
- ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
- ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG");
-
- const string &name = c_name(id);
-
- // Read element offsets and flags.
- Matrix<Quantity> aips_offset = c_offset(id);
- assert(aips_offset.shape().isEqual(IPosition(2, 3, aips_offset.ncolumn())));
-
- Matrix<Bool> aips_flag = c_flag(id);
- assert(aips_flag.shape().isEqual(IPosition(2, 2, aips_offset.ncolumn())));
-
-// TileConfig tile_config;
-// if(name != "LBA") readTileConfig(table, id);
-// transformToFieldCoordinates(tile_config, coordinate_system.axes);
-
- for(size_t i = 0; i < aips_offset.ncolumn(); ++i)
- {
- vector3r_t antenna_position = {
- aips_offset(0, i).getValue(),
- aips_offset(1, i).getValue(),
- aips_offset(2, i).getValue()
- };
- antenna_position = transformToFieldCoordinates(antenna_position, coordinate_system.axes);
-// std::cout << "antenna_position: " << antenna_position[0] << ", " << antenna_position[1] << ", " << antenna_position[2] << std::endl;
- Antenna::Ptr antenna;
- Antenna::CoordinateSystem antenna_coordinate_system{
- antenna_position,
- lofar_antenna_orientation
- };
- if (name == "LBA") {
- antenna = Element::Ptr(new Element(antenna_coordinate_system, element_response, id));
- } else {
- antenna = Element::Ptr(new Element(antenna_coordinate_system, element_response, id));
- // TODO
- // HBA, HBA0, HBA1
- // antenna = make_tile(name, coordinate_system, tile_config, element_response);
- }
-
- antenna->m_enabled[0] = !aips_flag(0, i);
- antenna->m_enabled[1] = !aips_flag(1, i);
- beam_former->add_antenna(antenna);
+BeamFormer::Ptr make_tile(std::string name,
+ Antenna::CoordinateSystem coordinate_system,
+ TileConfig tile_config, ElementResponse::Ptr) {}
+
+BeamFormer::Ptr readAntennaField(const Table &table, unsigned int id,
+ ElementResponse::Ptr element_response) {
+ Antenna::CoordinateSystem coordinate_system = readCoordinateSystem(table, id);
+ // std::cout << "coordinate_system: " << std::endl;
+ // std::cout << " axes.p: " << coordinate_system.axes.p[0] << ", " <<
+ // coordinate_system.axes.p[1] << ", " << coordinate_system.axes.p[2] <<
+ // std::endl; std::cout << " axes.q: " << coordinate_system.axes.q[0] <<
+ // ", " << coordinate_system.axes.q[1] << ", " <<
+ // coordinate_system.axes.q[2] << std::endl; std::cout << " axes.r: " <<
+ // coordinate_system.axes.r[0] << ", " << coordinate_system.axes.r[1] <<
+ // ", " << coordinate_system.axes.r[2] << std::endl;
+ BeamFormer::Ptr beam_former(new BeamFormer(coordinate_system));
+
+ ROScalarColumn<String> c_name(table, "NAME");
+ ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");
+ ROArrayColumn<Bool> c_flag(table, "ELEMENT_FLAG");
+
+ const string &name = c_name(id);
+
+ // Read element offsets and flags.
+ Matrix<Quantity> aips_offset = c_offset(id);
+ assert(aips_offset.shape().isEqual(IPosition(2, 3, aips_offset.ncolumn())));
+
+ Matrix<Bool> aips_flag = c_flag(id);
+ assert(aips_flag.shape().isEqual(IPosition(2, 2, aips_offset.ncolumn())));
+
+ // TileConfig tile_config;
+ // if(name != "LBA") readTileConfig(table, id);
+ // transformToFieldCoordinates(tile_config, coordinate_system.axes);
+
+ for (size_t i = 0; i < aips_offset.ncolumn(); ++i) {
+ vector3r_t antenna_position = {aips_offset(0, i).getValue(),
+ aips_offset(1, i).getValue(),
+ aips_offset(2, i).getValue()};
+ antenna_position =
+ transformToFieldCoordinates(antenna_position, coordinate_system.axes);
+ // std::cout << "antenna_position: " << antenna_position[0] << ", "
+ // << antenna_position[1] << ", " << antenna_position[2] <<
+ // std::endl;
+ Antenna::Ptr antenna;
+ Antenna::CoordinateSystem antenna_coordinate_system{
+ antenna_position, lofar_antenna_orientation};
+ if (name == "LBA") {
+ antenna = Element::Ptr(
+ new Element(antenna_coordinate_system, element_response, id));
+ } else {
+ antenna = Element::Ptr(
+ new Element(antenna_coordinate_system, element_response, id));
+ // TODO
+ // HBA, HBA0, HBA1
+ // antenna = make_tile(name, coordinate_system, tile_config,
+ // element_response);
}
- return beam_former;
-}
-
-
-BeamFormer::Ptr readAntennaFieldAartfaac(const Table &table, const string &ant_type,
- unsigned int id)
-{
- BeamFormer::Ptr field;
-// AntennaField::CoordinateSystem system = readCoordinateSystemAartfaac(table, id);
-//
-// if (ant_type == "LBA")
-// {
-// DualDipoleAntenna::Ptr model(new DualDipoleAntenna());
-// field = AntennaField::Ptr(new AntennaFieldLBA(ant_type, system, model));
-// }
-// else // HBA
-// {
-// // TODO: implement this
-// throw std::runtime_error("HBA for Aartfaac is not implemented yet.");
-// }
-//
-// // Add only one antenna to the field (no offset, always enabled)
-// AntennaField::Antenna antenna;
-// antenna.position[0] = 0.;
-// antenna.position[1] = 0.;
-// antenna.position[2] = 0.;
-// antenna.enabled[0] = true;
-// antenna.enabled[1] = true;
-//
-// field->addAntenna(antenna);
+ antenna->m_enabled[0] = !aips_flag(0, i);
+ antenna->m_enabled[1] = !aips_flag(1, i);
+ beam_former->add_antenna(antenna);
+ }
+ return beam_former;
+}
- return field;
+BeamFormer::Ptr readAntennaFieldAartfaac(const Table &table,
+ const string &ant_type,
+ unsigned int id) {
+ BeamFormer::Ptr field;
+ // AntennaField::CoordinateSystem system =
+ // readCoordinateSystemAartfaac(table, id);
+ //
+ // if (ant_type == "LBA")
+ // {
+ // DualDipoleAntenna::Ptr model(new DualDipoleAntenna());
+ // field = AntennaField::Ptr(new AntennaFieldLBA(ant_type, system,
+ // model));
+ // }
+ // else // HBA
+ // {
+ // // TODO: implement this
+ // throw std::runtime_error("HBA for Aartfaac is not implemented
+ // yet.");
+ // }
+ //
+ // // Add only one antenna to the field (no offset, always enabled)
+ // AntennaField::Antenna antenna;
+ // antenna.position[0] = 0.;
+ // antenna.position[1] = 0.;
+ // antenna.position[2] = 0.;
+ // antenna.enabled[0] = true;
+ // antenna.enabled[1] = true;
+ //
+ // field->addAntenna(antenna);
+
+ return field;
}
-vector3r_t readStationPhaseReference(const Table &table, unsigned int id)
-{
- vector3r_t phase_reference = {0.0, 0.0, 0.0};
- const string columnName("LOFAR_PHASE_REFERENCE");
- if(hasColumn(table, columnName))
- {
- ROScalarMeasColumn<MPosition> c_reference(table, columnName);
- MPosition mReference = MPosition::Convert(c_reference(id),
- MPosition::ITRF)();
- MVPosition mvReference = mReference.getValue();
- phase_reference = {mvReference(0), mvReference(1), mvReference(2)};
- }
- return phase_reference;
+vector3r_t readStationPhaseReference(const Table &table, unsigned int id) {
+ vector3r_t phase_reference = {0.0, 0.0, 0.0};
+ const string columnName("LOFAR_PHASE_REFERENCE");
+ if (hasColumn(table, columnName)) {
+ ROScalarMeasColumn<MPosition> c_reference(table, columnName);
+ MPosition mReference =
+ MPosition::Convert(c_reference(id), MPosition::ITRF)();
+ MVPosition mvReference = mReference.getValue();
+ phase_reference = {mvReference(0), mvReference(1), mvReference(2)};
+ }
+ return phase_reference;
}
-Station::Ptr readStation(
- const MeasurementSet &ms,
- unsigned int id,
- const ElementResponseModel model)
-{
- ROMSAntennaColumns antenna(ms.antenna());
- assert(antenna.nrow() > id && !antenna.flagRow()(id));
+Station::Ptr readStation(const MeasurementSet &ms, unsigned int id,
+ const ElementResponseModel model) {
+ ROMSAntennaColumns antenna(ms.antenna());
+ assert(antenna.nrow() > id && !antenna.flagRow()(id));
- // Get station name.
- const string name(antenna.name()(id));
+ // Get station name.
+ const string name(antenna.name()(id));
- // Get station position (ITRF).
- MPosition mPosition = MPosition::Convert(antenna.positionMeas()(id),
- MPosition::ITRF)();
- MVPosition mvPosition = mPosition.getValue();
- const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
+ // Get station position (ITRF).
+ MPosition mPosition =
+ MPosition::Convert(antenna.positionMeas()(id), MPosition::ITRF)();
+ MVPosition mvPosition = mPosition.getValue();
+ const vector3r_t position = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
- // Create station.
- Station::Ptr station(new Station(name, position, model));
+ // Create station.
+ Station::Ptr station(new Station(name, position, model));
- // Read phase reference position (if available).
- station->setPhaseReference(readStationPhaseReference(ms.antenna(), id));
+ // Read phase reference position (if available).
+ station->setPhaseReference(readStationPhaseReference(ms.antenna(), id));
- // Read antenna field information.
- ROScalarColumn<String> telescope_name_col(getSubTable(ms, "OBSERVATION"),
- "TELESCOPE_NAME");
- string telescope_name = telescope_name_col(0);
+ // Read antenna field information.
+ ROScalarColumn<String> telescope_name_col(getSubTable(ms, "OBSERVATION"),
+ "TELESCOPE_NAME");
+ string telescope_name = telescope_name_col(0);
- if (telescope_name == "LOFAR")
- {
- Table tab_field = getSubTable(ms, "LOFAR_ANTENNA_FIELD");
- tab_field = tab_field(tab_field.col("ANTENNA_ID") == static_cast<Int>(id));
-
- // The Station will consist of a BeamFormer that combines the fields
- // coordinate system is ITRF
- // phase reference is station position
- auto beam_former = BeamFormer::Ptr(new BeamFormer(
- Antenna::identity_coordinate_system,
- station->phaseReference()));
-
- for(size_t i = 0; i < tab_field.nrow(); ++i)
- {
- beam_former->add_antenna(readAntennaField(tab_field, i, station->get_element_response()));
- }
-
- // TODO
- // If There is only one field, the top level beamformer is not needed
- // and the station antenna can be set the the beamformer of the field
-
- station->set_antenna(beam_former);
-
- size_t field_id = 0;
- size_t element_id = 0;
- Antenna::CoordinateSystem coordinate_system = readCoordinateSystem(tab_field, field_id);
- auto model = station->get_element_response();
- // TODO: rotate coordinate system for antenna
- auto element = Element::Ptr(new Element(coordinate_system, model, element_id));
- station->set_element(element);
- }
- else if (telescope_name == "AARTFAAC")
- {
- ROScalarColumn<String> ant_type_col(getSubTable(ms, "OBSERVATION"),
- "AARTFAAC_ANTENNA_TYPE");
- string ant_type = ant_type_col(0);
+ if (telescope_name == "LOFAR") {
+ Table tab_field = getSubTable(ms, "LOFAR_ANTENNA_FIELD");
+ tab_field = tab_field(tab_field.col("ANTENNA_ID") == static_cast<Int>(id));
+
+ // The Station will consist of a BeamFormer that combines the fields
+ // coordinate system is ITRF
+ // phase reference is station position
+ auto beam_former = BeamFormer::Ptr(new BeamFormer(
+ Antenna::identity_coordinate_system, station->phaseReference()));
- Table tab_field = getSubTable(ms, "ANTENNA");
- station->set_antenna(readAntennaFieldAartfaac(tab_field, ant_type, id));
+ for (size_t i = 0; i < tab_field.nrow(); ++i) {
+ beam_former->add_antenna(
+ readAntennaField(tab_field, i, station->get_element_response()));
}
- return station;
+ // TODO
+ // If There is only one field, the top level beamformer is not needed
+ // and the station antenna can be set the the beamformer of the field
+
+ station->set_antenna(beam_former);
+
+ size_t field_id = 0;
+ size_t element_id = 0;
+ Antenna::CoordinateSystem coordinate_system =
+ readCoordinateSystem(tab_field, field_id);
+ auto model = station->get_element_response();
+ // TODO: rotate coordinate system for antenna
+ auto element =
+ Element::Ptr(new Element(coordinate_system, model, element_id));
+ station->set_element(element);
+ } else if (telescope_name == "AARTFAAC") {
+ ROScalarColumn<String> ant_type_col(getSubTable(ms, "OBSERVATION"),
+ "AARTFAAC_ANTENNA_TYPE");
+ string ant_type = ant_type_col(0);
+
+ Table tab_field = getSubTable(ms, "ANTENNA");
+ station->set_antenna(readAntennaFieldAartfaac(tab_field, ant_type, id));
+ }
+
+ return station;
}
MDirection readTileBeamDirection(const casacore::MeasurementSet &ms) {
- MDirection tileBeamDir;
-
- Table fieldTable = getSubTable(ms, "FIELD");
-
- if (fieldTable.nrow() != 1) {
- throw std::runtime_error("MS has multiple fields, this does not work with the LOFAR beam library.");
- }
-
- if (hasColumn(fieldTable, "LOFAR_TILE_BEAM_DIR"))
- {
- ROArrayMeasColumn<MDirection> tileBeamCol(fieldTable,
- "LOFAR_TILE_BEAM_DIR");
- tileBeamDir = *(tileBeamCol(0).data());
- }
- else
- {
- ROArrayMeasColumn<MDirection> tileBeamCol(fieldTable,
- "DELAY_DIR");
- tileBeamDir = *(tileBeamCol(0).data());
- }
-
- return tileBeamDir;
+ MDirection tileBeamDir;
+
+ Table fieldTable = getSubTable(ms, "FIELD");
+
+ if (fieldTable.nrow() != 1) {
+ throw std::runtime_error(
+ "MS has multiple fields, this does not work with the LOFAR beam "
+ "library.");
+ }
+
+ if (hasColumn(fieldTable, "LOFAR_TILE_BEAM_DIR")) {
+ ROArrayMeasColumn<MDirection> tileBeamCol(fieldTable,
+ "LOFAR_TILE_BEAM_DIR");
+ tileBeamDir = *(tileBeamCol(0).data());
+ } else {
+ ROArrayMeasColumn<MDirection> tileBeamCol(fieldTable, "DELAY_DIR");
+ tileBeamDir = *(tileBeamCol(0).data());
+ }
+
+ return tileBeamDir;
}
-} //# namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/MathUtil.cc b/MathUtil.cc
index bdff62b70438ac091a1975a7f73c1e238f4f92ad..4a73b10cc02277e0e0729f3aad53c3e4069f85ed 100644
--- a/MathUtil.cc
+++ b/MathUtil.cc
@@ -23,8 +23,5 @@
#include "MathUtil.h"
namespace LOFAR {
-namespace StationResponse
-{
-
-} //# namespace StationResponse
-} // namespace LOFAR
+namespace StationResponse {} // namespace StationResponse
+} // namespace LOFAR
diff --git a/Station.cc b/Station.cc
index f2fd43ff34be9507319a4651ba6a772414a82fcd..5723bdb6ff872bc9db4b4900e2b5073f9243aee8 100644
--- a/Station.cc
+++ b/Station.cc
@@ -30,240 +30,212 @@
// #include "TileAntenna.h"
namespace LOFAR {
-namespace StationResponse
-{
-
-Station::Station(
- const std::string &name,
- const vector3r_t &position,
- const ElementResponseModel model)
- : itsName(name),
- itsPosition(position),
- itsPhaseReference(position),
- itsElementResponse(nullptr)
-{
- setModel(model);
- vector3r_t ncp = {{0.0, 0.0, 1.0}};
- itsNCP.reset(new ITRFDirection(ncp));
- vector3r_t ncppol0 = {{1.0, 0.0, 0.0}};
- itsNCPPol0.reset(new ITRFDirection(ncppol0));
+namespace StationResponse {
+
+Station::Station(const std::string &name, const vector3r_t &position,
+ const ElementResponseModel model)
+ : itsName(name),
+ itsPosition(position),
+ itsPhaseReference(position),
+ itsElementResponse(nullptr) {
+ setModel(model);
+ vector3r_t ncp = {{0.0, 0.0, 1.0}};
+ itsNCP.reset(new ITRFDirection(ncp));
+ vector3r_t ncppol0 = {{1.0, 0.0, 0.0}};
+ itsNCPPol0.reset(new ITRFDirection(ncppol0));
}
-void Station::setModel(const ElementResponseModel model)
-{
- switch (model)
- {
- case Hamaker:
- itsElementResponse.set(HamakerElementResponse::getInstance(itsName));
- break;
- case OSKARDipole:
- itsElementResponse.set(OSKARElementResponseDipole::getInstance());
- break;
- case OSKARSphericalWave:
- itsElementResponse.set(OSKARElementResponseSphericalWave::getInstance());
- break;
- case LOBES:
- itsElementResponse.set(LOBESElementResponse::getInstance(itsName));
- break;
- default:
- std::stringstream message;
- message << "The requested element response model '"
- << model << "' is not implemented.";
- throw std::runtime_error(message.str());
- }
+void Station::setModel(const ElementResponseModel model) {
+ switch (model) {
+ case Hamaker:
+ itsElementResponse.set(HamakerElementResponse::getInstance(itsName));
+ break;
+ case OSKARDipole:
+ itsElementResponse.set(OSKARElementResponseDipole::getInstance());
+ break;
+ case OSKARSphericalWave:
+ itsElementResponse.set(OSKARElementResponseSphericalWave::getInstance());
+ break;
+ case LOBES:
+ itsElementResponse.set(LOBESElementResponse::getInstance(itsName));
+ break;
+ default:
+ std::stringstream message;
+ message << "The requested element response model '" << model
+ << "' is not implemented.";
+ throw std::runtime_error(message.str());
+ }
}
-const std::string &Station::name() const
-{
- return itsName;
-}
+const std::string &Station::name() const { return itsName; }
-const vector3r_t &Station::position() const
-{
- return itsPosition;
-}
+const vector3r_t &Station::position() const { return itsPosition; }
-void Station::setPhaseReference(const vector3r_t &reference)
-{
- itsPhaseReference = reference;
+void Station::setPhaseReference(const vector3r_t &reference) {
+ itsPhaseReference = reference;
}
-const vector3r_t &Station::phaseReference() const
-{
- return itsPhaseReference;
-}
+const vector3r_t &Station::phaseReference() const { return itsPhaseReference; }
// ========================================================
matrix22c_t Station::elementResponse(real_t time, real_t freq,
- const vector3r_t &direction, size_t id, const bool rotate) const
-{
- Antenna::Options options;
- options.rotate = rotate;
-
- if (rotate) {
- vector3r_t ncp_ = ncp(time);
- vector3r_t east = normalize(cross(ncp_, direction));
- vector3r_t north = cross(direction, east);
- options.east = east;
- options.north = north;
- }
+ const vector3r_t &direction, size_t id,
+ const bool rotate) const {
+ Antenna::Options options;
+ options.rotate = rotate;
+
+ if (rotate) {
+ vector3r_t ncp_ = ncp(time);
+ vector3r_t east = normalize(cross(ncp_, direction));
+ vector3r_t north = cross(direction, east);
+ options.east = east;
+ options.north = north;
+ }
return itsElement->local_response(time, freq, direction, id, options);
}
matrix22c_t Station::elementResponse(real_t time, real_t freq,
- const vector3r_t &direction, const bool rotate) const
-{
-// if (rotate)
-// return itsElement->response(time, freq, direction)
-// * rotation(time, direction);
-// else
-// return itsElement->response(time, freq, direction);
+ const vector3r_t &direction,
+ const bool rotate) const {
+ // if (rotate)
+ // return itsElement->response(time, freq, direction)
+ // * rotation(time, direction);
+ // else
+ // return itsElement->response(time, freq, direction);
return itsElement->response(time, freq, direction);
}
matrix22c_t Station::response(real_t time, real_t freq,
- const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
- const vector3r_t &tile0, const bool rotate) const
-{
- Antenna::Options options = {
- .freq0 = freq0,
- .station0 = station0,
- .tile0 = tile0,
- .rotate = rotate
- };
-
- if (rotate) {
- vector3r_t ncp_ = ncp(time);
- vector3r_t east = normalize(cross(ncp_, direction));
- vector3r_t north = cross(direction, east);
- options.east = east;
- options.north = north;
- }
-
- matrix22c_t response = itsAntenna->response(time, freq, direction, options);
-
-// if (rotate) {
-// std::cout << "rotate" << std::endl;
-// auto r = rotation(time, direction);
-// std::cout << r[0][0] << ", " << r[0][1] << std::endl;
-// std::cout << r[1][0] << ", " << r[1][1] << std::endl;
-// response = response * r;
-// std::cout << response[0][0] << std::endl;
-// }
- return response;
+ const vector3r_t &direction, real_t freq0,
+ const vector3r_t &station0,
+ const vector3r_t &tile0,
+ const bool rotate) const {
+ Antenna::Options options = {
+ .freq0 = freq0, .station0 = station0, .tile0 = tile0, .rotate = rotate};
+
+ if (rotate) {
+ vector3r_t ncp_ = ncp(time);
+ vector3r_t east = normalize(cross(ncp_, direction));
+ vector3r_t north = cross(direction, east);
+ options.east = east;
+ options.north = north;
+ }
+
+ matrix22c_t response = itsAntenna->response(time, freq, direction, options);
+
+ // if (rotate) {
+ // std::cout << "rotate" << std::endl;
+ // auto r = rotation(time, direction);
+ // std::cout << r[0][0] << ", " << r[0][1] << std::endl;
+ // std::cout << r[1][0] << ", " << r[1][1] << std::endl;
+ // response = response * r;
+ // std::cout << response[0][0] << std::endl;
+ // }
+ return response;
}
diag22c_t Station::arrayFactor(real_t time, real_t freq,
- const vector3r_t &direction, real_t freq0, const vector3r_t &station0,
- const vector3r_t &tile0) const
-{
- Antenna::Options options = {
- .freq0 = freq0,
- .station0 = station0,
- .tile0 = tile0};
- return itsAntenna->arrayFactor(time, freq, direction, options);
-}
-
-matrix22r_t Station::rotation(real_t time, const vector3r_t &direction)
- const
-{
- //rotation needs to be optional, normally you only want to rotate your coordinatesytem for the center of your (mosaiced) image
- // Compute the cross product of the NCP and the target direction. This
- // yields a vector tangent to the celestial sphere at the target
- // direction, pointing towards the East (the direction of +Y in the IAU
- // definition, or positive right ascension).
- // Test if the direction is equal to the NCP. If it is, take a random
- // vector orthogonal to v1 (the east is not defined here).
- vector3r_t v1;
- if (std::abs(ncp(time)[0]-direction[0])<1e-9 &&
- std::abs(ncp(time)[1]-direction[1])<1e-9 &&
- std::abs(ncp(time)[2]-direction[2])<1e-9) {
- // Make sure v1 is orthogonal to ncp(time). In the direction of the meridian
- v1 = normalize(ncppol0(time));
- } else {
- v1 = normalize(cross(ncp(time), direction));
- }
-
- // Compute the cross product of the antenna field normal (R) and the
- // target direction. This yields a vector tangent to the topocentric
- // spherical coordinate system at the target direction, pointing towards
- // the direction of positive phi (which runs East over North around the
- // pseudo zenith).
- // Test if the normal is equal to the target direction. If it is, take
- // a random vector orthogonal to the normal.
- vector3r_t v2;
- if (std::abs(itsAntenna->m_coordinate_system.axes.r[0]-direction[0])<1e-9 &&
- std::abs(itsAntenna->m_coordinate_system.axes.r[1]-direction[1])<1e-9 &&
- std::abs(itsAntenna->m_coordinate_system.axes.r[2]-direction[2])<1e-9)
- {
- // Nothing to be rotated if the direction is equal to zenith
- v2 = v1;
- } else {
- v2 = normalize(cross(itsAntenna->m_coordinate_system.axes.r, direction));
- }
-
- // Compute the cosine and sine of the parallactic angle, i.e. the angle
- // between v1 and v2, both tangent to a latitude circle of their
- // respective spherical coordinate systems.
- real_t coschi = dot(v1, v2);
- real_t sinchi;
- if (coschi==1.0)
- sinchi = 0.0;
- else
- sinchi = dot(cross(v1, v2), direction);
-
- // The input coordinate system is a right handed system with its third
- // axis along the direction of propagation (IAU +Z). The output
- // coordinate system is right handed as well, but its third axis points
- // in the direction of arrival (i.e. exactly opposite).
- //
- // Because the electromagnetic field is always perpendicular to the
- // direction of propagation, we only need to relate the (X, Y) axes of
- // the input system to the corresponding (theta, phi) axes of the output
- // system.
- //
- // To this end, we first rotate the input system around its third axis
- // to align the Y axis with the phi axis. The X and theta axis are
- // parallel after this rotation, but point in opposite directions. To
- // align the X axis with the theta axis, we flip it.
- //
- // The Jones matrix to align the Y axis with the phi axis when these are
- // separated by an angle phi (measured counter-clockwise around the
- // direction of propagation, looking towards the origin), is given by:
- //
- // [ cos(phi) sin(phi)]
- // [-sin(phi) cos(phi)]
- //
- // Here, cos(phi) and sin(phi) can be computed directly, without having
- // to compute phi first (see the computation of coschi and sinchi
- // above).
- //
- // Now, sinchi as computed above is opposite to sin(phi), because the
- // direction used in the computation is the direction of arrival instead
- // of the direction of propagation. Therefore, the sign of sinchi needs
- // to be reversed. Furthermore, as explained above, the X axis has to be
- // flipped to align with the theta axis. The Jones matrix returned from
- // this function is therefore given by:
- //
- // [-coschi sinchi]
- // [ sinchi coschi]
- matrix22r_t rotation = {{{{-coschi, sinchi}}, {{sinchi, coschi}}}};
- return rotation;
-}
-
-vector3r_t Station::ncp(real_t time) const
-{
- return itsNCP->at(time);
+ const vector3r_t &direction, real_t freq0,
+ const vector3r_t &station0,
+ const vector3r_t &tile0) const {
+ Antenna::Options options = {
+ .freq0 = freq0, .station0 = station0, .tile0 = tile0};
+ return itsAntenna->arrayFactor(time, freq, direction, options);
}
-vector3r_t Station::ncppol0(real_t time) const
-{
- return itsNCPPol0->at(time);
+matrix22r_t Station::rotation(real_t time, const vector3r_t &direction) const {
+ // rotation needs to be optional, normally you only want to rotate your
+ // coordinatesytem for the center of your (mosaiced) image
+ // Compute the cross product of the NCP and the target direction. This
+ // yields a vector tangent to the celestial sphere at the target
+ // direction, pointing towards the East (the direction of +Y in the IAU
+ // definition, or positive right ascension).
+ // Test if the direction is equal to the NCP. If it is, take a random
+ // vector orthogonal to v1 (the east is not defined here).
+ vector3r_t v1;
+ if (std::abs(ncp(time)[0] - direction[0]) < 1e-9 &&
+ std::abs(ncp(time)[1] - direction[1]) < 1e-9 &&
+ std::abs(ncp(time)[2] - direction[2]) < 1e-9) {
+ // Make sure v1 is orthogonal to ncp(time). In the direction of the meridian
+ v1 = normalize(ncppol0(time));
+ } else {
+ v1 = normalize(cross(ncp(time), direction));
+ }
+
+ // Compute the cross product of the antenna field normal (R) and the
+ // target direction. This yields a vector tangent to the topocentric
+ // spherical coordinate system at the target direction, pointing towards
+ // the direction of positive phi (which runs East over North around the
+ // pseudo zenith).
+ // Test if the normal is equal to the target direction. If it is, take
+ // a random vector orthogonal to the normal.
+ vector3r_t v2;
+ if (std::abs(itsAntenna->m_coordinate_system.axes.r[0] - direction[0]) <
+ 1e-9 &&
+ std::abs(itsAntenna->m_coordinate_system.axes.r[1] - direction[1]) <
+ 1e-9 &&
+ std::abs(itsAntenna->m_coordinate_system.axes.r[2] - direction[2]) <
+ 1e-9) {
+ // Nothing to be rotated if the direction is equal to zenith
+ v2 = v1;
+ } else {
+ v2 = normalize(cross(itsAntenna->m_coordinate_system.axes.r, direction));
+ }
+
+ // Compute the cosine and sine of the parallactic angle, i.e. the angle
+ // between v1 and v2, both tangent to a latitude circle of their
+ // respective spherical coordinate systems.
+ real_t coschi = dot(v1, v2);
+ real_t sinchi;
+ if (coschi == 1.0)
+ sinchi = 0.0;
+ else
+ sinchi = dot(cross(v1, v2), direction);
+
+ // The input coordinate system is a right handed system with its third
+ // axis along the direction of propagation (IAU +Z). The output
+ // coordinate system is right handed as well, but its third axis points
+ // in the direction of arrival (i.e. exactly opposite).
+ //
+ // Because the electromagnetic field is always perpendicular to the
+ // direction of propagation, we only need to relate the (X, Y) axes of
+ // the input system to the corresponding (theta, phi) axes of the output
+ // system.
+ //
+ // To this end, we first rotate the input system around its third axis
+ // to align the Y axis with the phi axis. The X and theta axis are
+ // parallel after this rotation, but point in opposite directions. To
+ // align the X axis with the theta axis, we flip it.
+ //
+ // The Jones matrix to align the Y axis with the phi axis when these are
+ // separated by an angle phi (measured counter-clockwise around the
+ // direction of propagation, looking towards the origin), is given by:
+ //
+ // [ cos(phi) sin(phi)]
+ // [-sin(phi) cos(phi)]
+ //
+ // Here, cos(phi) and sin(phi) can be computed directly, without having
+ // to compute phi first (see the computation of coschi and sinchi
+ // above).
+ //
+ // Now, sinchi as computed above is opposite to sin(phi), because the
+ // direction used in the computation is the direction of arrival instead
+ // of the direction of propagation. Therefore, the sign of sinchi needs
+ // to be reversed. Furthermore, as explained above, the X axis has to be
+ // flipped to align with the theta axis. The Jones matrix returned from
+ // this function is therefore given by:
+ //
+ // [-coschi sinchi]
+ // [ sinchi coschi]
+ matrix22r_t rotation = {{{{-coschi, sinchi}}, {{sinchi, coschi}}}};
+ return rotation;
}
+vector3r_t Station::ncp(real_t time) const { return itsNCP->at(time); }
+vector3r_t Station::ncppol0(real_t time) const { return itsNCPPol0->at(time); }
-} //# namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/Types.cc b/Types.cc
index 43c926234172556444a484929cba395109f9f604..ccfb3c1ba933e64c155eafee8be58c9373150632 100644
--- a/Types.cc
+++ b/Types.cc
@@ -23,8 +23,5 @@
#include "Types.h"
namespace LOFAR {
-namespace StationResponse
-{
-
-} //# namespace StationResponse
-} // namespace LOFAR
+namespace StationResponse {} // namespace StationResponse
+} // namespace LOFAR
diff --git a/hamaker/HamakerCoeff.cc b/hamaker/HamakerCoeff.cc
index 1becc96ab8eba8b2c1cd5feaeee0670e5f37d08d..66ba45ddce4b74db1f75f0eb647631fe8ce0398c 100644
--- a/hamaker/HamakerCoeff.cc
+++ b/hamaker/HamakerCoeff.cc
@@ -1,163 +1,145 @@
#include "HamakerCoeff.h"
-H5::CompType get_complex_double_type()
-{
- H5::CompType complex_type(sizeof(std::complex<double>));
- complex_type.insertMember("r", 0, H5::PredType::NATIVE_DOUBLE);
- complex_type.insertMember("i", sizeof(double), H5::PredType::NATIVE_DOUBLE);
- return complex_type;
+H5::CompType get_complex_double_type() {
+ H5::CompType complex_type(sizeof(std::complex<double>));
+ complex_type.insertMember("r", 0, H5::PredType::NATIVE_DOUBLE);
+ complex_type.insertMember("i", sizeof(double), H5::PredType::NATIVE_DOUBLE);
+ return complex_type;
}
-size_t HamakerCoefficients::get_index(
- const unsigned int n,
- const unsigned int t,
- const unsigned int f)
-{
- return n * m_nPowerTheta * m_nPowerFreq * m_nInner +
- t * m_nPowerFreq * m_nInner +
- f * m_nInner;
+size_t HamakerCoefficients::get_index(const unsigned int n,
+ const unsigned int t,
+ const unsigned int f) {
+ return n * m_nPowerTheta * m_nPowerFreq * m_nInner +
+ t * m_nPowerFreq * m_nInner + f * m_nInner;
}
// Constructor for reading coeff from file
-HamakerCoefficients::HamakerCoefficients(
- std::string& filename)
-{
- read_coeffs(filename);
+HamakerCoefficients::HamakerCoefficients(std::string& filename) {
+ read_coeffs(filename);
};
// Constructor for writing coeff to file
-HamakerCoefficients::HamakerCoefficients(
- const double freq_center,
- const double freq_range,
- const unsigned int nHarmonics,
- const unsigned int nPowerTheta,
- const unsigned int nPowerFreq) :
- m_freq_center(freq_center),
- m_freq_range(freq_range),
- m_nHarmonics(nHarmonics),
- m_nPowerTheta(nPowerTheta),
- m_nPowerFreq(nPowerFreq),
- m_coeff(get_nr_coeffs())
-{}
-
-size_t HamakerCoefficients::get_nr_coeffs() const
-{
- return m_nHarmonics * m_nPowerTheta * m_nPowerFreq * m_nInner;
+HamakerCoefficients::HamakerCoefficients(const double freq_center,
+ const double freq_range,
+ const unsigned int nHarmonics,
+ const unsigned int nPowerTheta,
+ const unsigned int nPowerFreq)
+ : m_freq_center(freq_center),
+ m_freq_range(freq_range),
+ m_nHarmonics(nHarmonics),
+ m_nPowerTheta(nPowerTheta),
+ m_nPowerFreq(nPowerFreq),
+ m_coeff(get_nr_coeffs()) {}
+
+size_t HamakerCoefficients::get_nr_coeffs() const {
+ return m_nHarmonics * m_nPowerTheta * m_nPowerFreq * m_nInner;
}
void HamakerCoefficients::set_coeff(
- const unsigned int n,
- const unsigned int t,
- const unsigned int f,
- std::pair<std::complex<double>, std::complex<double>> value)
-{
- size_t index = get_index(n, t, f);
- m_coeff[index] = value.first;
- m_coeff[index + 1] = value.second;
+ const unsigned int n, const unsigned int t, const unsigned int f,
+ std::pair<std::complex<double>, std::complex<double>> value) {
+ size_t index = get_index(n, t, f);
+ m_coeff[index] = value.first;
+ m_coeff[index + 1] = value.second;
}
-void HamakerCoefficients::set_coeffs(
- const std::complex<double>* coeff) {
- memcpy(m_coeff.data(), coeff, m_coeff.size() * sizeof(std::complex<double>));
+void HamakerCoefficients::set_coeffs(const std::complex<double>* coeff) {
+ memcpy(m_coeff.data(), coeff, m_coeff.size() * sizeof(std::complex<double>));
}
void HamakerCoefficients::set_coeffs(
- const std::vector<std::complex<double>> coeff)
-{
- assert(coeff.size() == m_coeff.size());
- std::copy(coeff.begin(), coeff.end(), m_coeff.begin());
+ const std::vector<std::complex<double>> coeff) {
+ assert(coeff.size() == m_coeff.size());
+ std::copy(coeff.begin(), coeff.end(), m_coeff.begin());
}
-std::pair<std::complex<double>, std::complex<double>> HamakerCoefficients::get_coeff(
- const unsigned int n,
- const unsigned int t,
- const unsigned int f)
-{
- size_t index = get_index(n, t, f);
- std::pair<std::complex<double>, std::complex<double>> value;
- value.first = m_coeff[index];
- value.second = m_coeff[index + 1];
- return value;
+std::pair<std::complex<double>, std::complex<double>>
+ HamakerCoefficients::get_coeff(const unsigned int n, const unsigned int t,
+ const unsigned int f) {
+ size_t index = get_index(n, t, f);
+ std::pair<std::complex<double>, std::complex<double>> value;
+ value.first = m_coeff[index];
+ value.second = m_coeff[index + 1];
+ return value;
}
-void HamakerCoefficients::read_coeffs(
- std::string& filename)
-{
- // Open file
- H5::H5File file(filename, H5F_ACC_RDONLY);
-
- // Read dataset
- H5::DataSet dataset = file.openDataSet(m_dataset_name);
-
- // Open attribute and read its contents
- H5::Attribute freq_center_attr = dataset.openAttribute("freq_center");
- H5::Attribute freq_range_attr = dataset.openAttribute("freq_range");
- freq_center_attr.read(H5::PredType::NATIVE_DOUBLE, &m_freq_center);
- freq_range_attr.read(H5::PredType::NATIVE_DOUBLE, &m_freq_range);
-
- // Read dataset dimensions
- H5::DataSpace dataspace = dataset.getSpace();
- unsigned int rank = dataspace.getSimpleExtentNdims();
- assert(rank == m_dataset_rank);
- hsize_t dims[rank];
- dataspace.getSimpleExtentDims(dims, NULL);
- m_nHarmonics = dims[0];
- m_nPowerTheta = dims[1];
- m_nPowerFreq = dims[2];
-
- // Read coeffs
- m_coeff.resize(get_nr_coeffs());
- H5::DataType data_type = dataset.getDataType();
- assert(data_type.getSize() == sizeof(std::complex<double>));
- dataset.read(m_coeff.data(), data_type, dataspace);
+void HamakerCoefficients::read_coeffs(std::string& filename) {
+ // Open file
+ H5::H5File file(filename, H5F_ACC_RDONLY);
+
+ // Read dataset
+ H5::DataSet dataset = file.openDataSet(m_dataset_name);
+
+ // Open attribute and read its contents
+ H5::Attribute freq_center_attr = dataset.openAttribute("freq_center");
+ H5::Attribute freq_range_attr = dataset.openAttribute("freq_range");
+ freq_center_attr.read(H5::PredType::NATIVE_DOUBLE, &m_freq_center);
+ freq_range_attr.read(H5::PredType::NATIVE_DOUBLE, &m_freq_range);
+
+ // Read dataset dimensions
+ H5::DataSpace dataspace = dataset.getSpace();
+ unsigned int rank = dataspace.getSimpleExtentNdims();
+ assert(rank == m_dataset_rank);
+ hsize_t dims[rank];
+ dataspace.getSimpleExtentDims(dims, NULL);
+ m_nHarmonics = dims[0];
+ m_nPowerTheta = dims[1];
+ m_nPowerFreq = dims[2];
+
+ // Read coeffs
+ m_coeff.resize(get_nr_coeffs());
+ H5::DataType data_type = dataset.getDataType();
+ assert(data_type.getSize() == sizeof(std::complex<double>));
+ dataset.read(m_coeff.data(), data_type, dataspace);
}
-void HamakerCoefficients::write_coeffs(
- std::string& filename)
-{
- // Open file
- H5::H5File file(filename, H5F_ACC_TRUNC);
+void HamakerCoefficients::write_coeffs(std::string& filename) {
+ // Open file
+ H5::H5File file(filename, H5F_ACC_TRUNC);
- // Create dataspace
- const unsigned int rank = 4;
- hsize_t dims[rank] = { m_nHarmonics, m_nPowerTheta, m_nPowerFreq, m_nInner };
- H5::DataSpace coeff_dataspace(rank, dims);
+ // Create dataspace
+ const unsigned int rank = 4;
+ hsize_t dims[rank] = {m_nHarmonics, m_nPowerTheta, m_nPowerFreq, m_nInner};
+ H5::DataSpace coeff_dataspace(rank, dims);
- // Create pointer to coeff
- typedef std::complex<double> coeff_type[m_nHarmonics][m_nPowerTheta][m_nPowerFreq][m_nInner];
- coeff_type *coeff_ptr = (coeff_type *) m_coeff.data();
+ // Create pointer to coeff
+ typedef std::complex<double> coeff_type[m_nHarmonics][m_nPowerTheta]
+ [m_nPowerFreq][m_nInner];
+ coeff_type* coeff_ptr = (coeff_type*)m_coeff.data();
- // Create complex type
- H5::CompType complex_type = get_complex_double_type();
+ // Create complex type
+ H5::CompType complex_type = get_complex_double_type();
- // Write dataset
- H5::DataSet dataset = file.createDataSet("coeff", complex_type, coeff_dataspace);
- dataset.write(coeff_ptr, complex_type);
+ // Write dataset
+ H5::DataSet dataset =
+ file.createDataSet("coeff", complex_type, coeff_dataspace);
+ dataset.write(coeff_ptr, complex_type);
- // Create new dataspace for attribute
- H5::DataSpace attr_dataspace(H5S_SCALAR);
+ // Create new dataspace for attribute
+ H5::DataSpace attr_dataspace(H5S_SCALAR);
- // Write frequency center attribute
- H5::Attribute freq_center_attr = dataset.createAttribute("freq_center", H5::PredType::NATIVE_DOUBLE, attr_dataspace);
- freq_center_attr.write(H5::PredType::NATIVE_DOUBLE, &m_freq_center);
+ // Write frequency center attribute
+ H5::Attribute freq_center_attr = dataset.createAttribute(
+ "freq_center", H5::PredType::NATIVE_DOUBLE, attr_dataspace);
+ freq_center_attr.write(H5::PredType::NATIVE_DOUBLE, &m_freq_center);
- // Write frequency range attribute
- H5::Attribute freq_range_attr = dataset.createAttribute("freq_range", H5::PredType::NATIVE_DOUBLE, attr_dataspace);
- freq_range_attr.write(H5::PredType::NATIVE_DOUBLE, &m_freq_range);
+ // Write frequency range attribute
+ H5::Attribute freq_range_attr = dataset.createAttribute(
+ "freq_range", H5::PredType::NATIVE_DOUBLE, attr_dataspace);
+ freq_range_attr.write(H5::PredType::NATIVE_DOUBLE, &m_freq_range);
- file.flush(H5F_SCOPE_LOCAL);
+ file.flush(H5F_SCOPE_LOCAL);
}
-void HamakerCoefficients::print_coeffs()
-{
- for (unsigned int h = 0; h < m_nHarmonics; h++) {
- for (unsigned int t = 0; t < m_nPowerTheta; t++) {
- for (unsigned int f = 0; f < m_nPowerFreq; f++) {
- auto coeff = get_coeff(h, t, f);
- std::cout << coeff.first << ", "
- << coeff.second << std::endl;
- }
- }
+void HamakerCoefficients::print_coeffs() {
+ for (unsigned int h = 0; h < m_nHarmonics; h++) {
+ for (unsigned int t = 0; t < m_nPowerTheta; t++) {
+ for (unsigned int f = 0; f < m_nPowerFreq; f++) {
+ auto coeff = get_coeff(h, t, f);
+ std::cout << coeff.first << ", " << coeff.second << std::endl;
+ }
}
- std::cout << std::endl;
+ }
+ std::cout << std::endl;
}
diff --git a/hamaker/HamakerElementResponse.cc b/hamaker/HamakerElementResponse.cc
index eb47fd59e9785c000f50b6e0fab20960de3e2043..34a8a5f23339a1e7fccfff149705a656b024641f 100644
--- a/hamaker/HamakerElementResponse.cc
+++ b/hamaker/HamakerElementResponse.cc
@@ -2,8 +2,7 @@
//# Functions to compute the (idealized) response of a LOFAR
//# LBA or HBA dual dipole antenna.
-
-#include<stdexcept>
+#include <stdexcept>
#include "config.h"
@@ -13,127 +12,117 @@
namespace LOFAR {
namespace StationResponse {
-
-std::shared_ptr<HamakerElementResponse> HamakerElementResponse::getInstance(const std::string &name)
-{
- if (name.find("LBA") != std::string::npos) {
- return Singleton<HamakerElementResponseLBA>::getInstance();
- }
- if (name.find("HBA") != std::string::npos) {
- return Singleton<HamakerElementResponseHBA>::getInstance();
- }
- throw std::invalid_argument("HamakerElementResponse::getInstance: name should end in either 'LBA' or 'HBA'");
+std::shared_ptr<HamakerElementResponse> HamakerElementResponse::getInstance(
+ const std::string& name) {
+ if (name.find("LBA") != std::string::npos) {
+ return Singleton<HamakerElementResponseLBA>::getInstance();
+ }
+ if (name.find("HBA") != std::string::npos) {
+ return Singleton<HamakerElementResponseHBA>::getInstance();
+ }
+ throw std::invalid_argument(
+ "HamakerElementResponse::getInstance: name should end in either 'LBA' or "
+ "'HBA'");
}
-
-std::string HamakerElementResponse::get_path(
- const char* filename) const
-{
- std::stringstream ss;
- ss << LOFARBEAM_DATA_DIR << "/";
- ss << filename;
- return ss.str();
+std::string HamakerElementResponse::get_path(const char* filename) const {
+ std::stringstream ss;
+ ss << LOFARBEAM_DATA_DIR << "/";
+ ss << filename;
+ return ss.str();
}
void HamakerElementResponse::response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const
-{
- // Initialize the response to zero.
- response[0][0] = 0.0;
- response[0][1] = 0.0;
- response[1][0] = 0.0;
- response[1][1] = 0.0;
-
- // Clip directions below the horizon.
- if(theta >= M_PI_2)
- {
- return;
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const {
+ // Initialize the response to zero.
+ response[0][0] = 0.0;
+ response[0][1] = 0.0;
+ response[1][0] = 0.0;
+ response[1][1] = 0.0;
+
+ // Clip directions below the horizon.
+ if (theta >= M_PI_2) {
+ return;
+ }
+
+ const double freq_center = m_coeffs->get_freq_center();
+ const double freq_range = m_coeffs->get_freq_range();
+ const unsigned int nHarmonics = m_coeffs->get_nHarmonics();
+ const unsigned int nPowerTheta = m_coeffs->get_nPowerTheta();
+ const unsigned int nPowerFreq = m_coeffs->get_nPowerFreq();
+ ;
+
+ // The model is parameterized in terms of a normalized frequency in the
+ // range [-1, 1]. The appropriate conversion is taken care of below.
+ freq = (freq - freq_center) / freq_range;
+
+ // The variables sign and kappa are used to compute the value of kappa
+ // mentioned in the description of the beam model [kappa = (-1)^k * (2 * k
+ //+ 1)] incrementally.
+ int sign = 1, kappa = 1;
+
+ std::pair<std::complex<double>, std::complex<double>> P;
+ std::pair<std::complex<double>, std::complex<double>> Pj;
+ for (unsigned int k = 0; k < nHarmonics; ++k) {
+ // Compute the (diagonal) projection matrix P for the current harmonic.
+ // This requires the evaluation of two polynomials in theta and freq (of
+ // degree nPowerTheta in theta and nPowerFreq in freq), one for each
+ // element of P. The polynomials are evaluated using Horner's rule.
+
+ // Horner's rule requires backward iteration of the coefficients, so
+ // start indexing the block of coefficients at the last element
+
+ // Evaluate the highest order term.
+ P = m_coeffs->get_coeff(k, nPowerTheta - 1, nPowerFreq - 1);
+
+ for (unsigned int i = 0; i < nPowerFreq - 1; ++i) {
+ auto Pk = m_coeffs->get_coeff(k, nPowerTheta - 1, nPowerFreq - i - 2);
+ P.first = P.first * freq + Pk.first;
+ P.second = P.second * freq + Pk.second;
}
- const double freq_center = m_coeffs->get_freq_center();
- const double freq_range = m_coeffs->get_freq_range();
- const unsigned int nHarmonics = m_coeffs->get_nHarmonics();
- const unsigned int nPowerTheta = m_coeffs->get_nPowerTheta();
- const unsigned int nPowerFreq = m_coeffs->get_nPowerFreq();;
-
- // The model is parameterized in terms of a normalized frequency in the
- // range [-1, 1]. The appropriate conversion is taken care of below.
- freq = (freq - freq_center) / freq_range;
-
- // The variables sign and kappa are used to compute the value of kappa
- // mentioned in the description of the beam model [kappa = (-1)^k * (2 * k
- //+ 1)] incrementally.
- int sign = 1, kappa = 1;
-
- std::pair<std::complex<double>, std::complex<double>> P;
- std::pair<std::complex<double>, std::complex<double>> Pj;
- for(unsigned int k = 0; k < nHarmonics; ++k)
- {
- // Compute the (diagonal) projection matrix P for the current harmonic.
- // This requires the evaluation of two polynomials in theta and freq (of
- // degree nPowerTheta in theta and nPowerFreq in freq), one for each
- // element of P. The polynomials are evaluated using Horner's rule.
-
- // Horner's rule requires backward iteration of the coefficients, so
- // start indexing the block of coefficients at the last element
-
- // Evaluate the highest order term.
- P = m_coeffs->get_coeff(k, nPowerTheta-1, nPowerFreq-1);
-
- for(unsigned int i = 0; i < nPowerFreq - 1; ++i)
- {
- auto Pk = m_coeffs->get_coeff(k, nPowerTheta-1, nPowerFreq-i-2);
- P.first = P.first * freq + Pk.first;
- P.second = P.second * freq + Pk.second;
- }
-
- // Evaluate the remaining terms.
- for(unsigned int j = 0; j < nPowerTheta - 1; ++j)
- {
- Pj = m_coeffs->get_coeff(k, nPowerTheta-j-2, nPowerFreq-1);
-
- for(unsigned int i = 0; i < nPowerFreq - 1; ++i)
- {
- auto Pk = m_coeffs->get_coeff(k, nPowerTheta-j-2, nPowerFreq-i-2);
- Pj.first = Pj.first * freq + Pk.first;
- Pj.second = Pj.second * freq + Pk.second;
- }
-
- P.first = P.first * theta + Pj.first;
- P.second = P.second * theta + Pj.second;
- }
-
- // Compute the Jones matrix for the current harmonic, by rotating P over
- // kappa * az, and add it to the result.
- const double angle = sign * kappa * phi;
- const double caz = std::cos(angle);
- const double saz = std::sin(angle);
-
- response[0][0] += caz * P.first;
- response[0][1] += -saz * P.second;
- response[1][0] += saz * P.first;
- response[1][1] += caz * P.second;
-
- // Update sign and kappa.
- sign = -sign;
- kappa += 2;
+ // Evaluate the remaining terms.
+ for (unsigned int j = 0; j < nPowerTheta - 1; ++j) {
+ Pj = m_coeffs->get_coeff(k, nPowerTheta - j - 2, nPowerFreq - 1);
+
+ for (unsigned int i = 0; i < nPowerFreq - 1; ++i) {
+ auto Pk =
+ m_coeffs->get_coeff(k, nPowerTheta - j - 2, nPowerFreq - i - 2);
+ Pj.first = Pj.first * freq + Pk.first;
+ Pj.second = Pj.second * freq + Pk.second;
+ }
+
+ P.first = P.first * theta + Pj.first;
+ P.second = P.second * theta + Pj.second;
}
+
+ // Compute the Jones matrix for the current harmonic, by rotating P over
+ // kappa * az, and add it to the result.
+ const double angle = sign * kappa * phi;
+ const double caz = std::cos(angle);
+ const double saz = std::sin(angle);
+
+ response[0][0] += caz * P.first;
+ response[0][1] += -saz * P.second;
+ response[1][0] += saz * P.first;
+ response[1][1] += caz * P.second;
+
+ // Update sign and kappa.
+ sign = -sign;
+ kappa += 2;
+ }
}
-HamakerElementResponseHBA::HamakerElementResponseHBA()
-{
- std::string path = get_path("HamakerHBACoeff.h5");
- m_coeffs.reset(new HamakerCoefficients(path));
+HamakerElementResponseHBA::HamakerElementResponseHBA() {
+ std::string path = get_path("HamakerHBACoeff.h5");
+ m_coeffs.reset(new HamakerCoefficients(path));
}
-HamakerElementResponseLBA::HamakerElementResponseLBA()
-{
- std::string path = get_path("HamakerLBACoeff.h5");
- m_coeffs.reset(new HamakerCoefficients(path));
+HamakerElementResponseLBA::HamakerElementResponseLBA() {
+ std::string path = get_path("HamakerLBACoeff.h5");
+ m_coeffs.reset(new HamakerCoefficients(path));
}
-} // namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/lobes/DefaultCoeffHBA.cc b/lobes/DefaultCoeffHBA.cc
index c37da7e5e27a609311dcc132013ba358ab9d3eb7..f57c8f720fa5cb6f16ea75f8b13fda4bacb893fb 100644
--- a/lobes/DefaultCoeffHBA.cc
+++ b/lobes/DefaultCoeffHBA.cc
@@ -1,5 +1,5 @@
// Beam model coefficients converted by convert_coeff.py.
-// Conversion performed on 2011/09/14/08:23:16 UTC using:
+// Conversion performed on 2011/09/14/08:23:16 UTC using:
// convert_coeff.py element_beam_HBA.coeff DefaultCoeffHBA.cc default_hba
#include <complex>
@@ -21,54 +21,103 @@ const unsigned int default_hba_coeff_shape[3] = {2, 5, 5};
// The array of coefficients in row-major order ("C"-order).
//
const std::complex<double> default_hba_coeff[100] = {
- std::complex<double>(0.9989322499459223, 0.0003305895124867), std::complex<double>(1.0030546028872600, 0.0002157249025076),
- std::complex<double>(0.0003002209532403, 0.0007909077657054), std::complex<double>(0.0022051270911392, 0.0003834815341981),
- std::complex<double>(-0.0003856663268042, 0.0008435910525861), std::complex<double>(0.0004887765294093, 0.0002777796480946),
- std::complex<double>(-0.0000699366665322, 0.0005136144371953), std::complex<double>(0.0001520602842105, 0.0001303481681886),
- std::complex<double>(0.0000512381993616, 0.0001550785137302), std::complex<double>(0.0000819244737818, 0.0000466470412396),
- std::complex<double>(0.0249658150445263, -0.0122024663463393), std::complex<double>(-0.0917825091832822, -0.0062606338208358),
- std::complex<double>(-0.0083709499453879, -0.0289759752488368), std::complex<double>(-0.0689260153643395, -0.0111348626546314),
- std::complex<double>(0.0116296166994115, -0.0307342946951178), std::complex<double>(-0.0171249717275797, -0.0080642275561593),
- std::complex<double>(0.0012408055399100, -0.0191295543986957), std::complex<double>(-0.0051031652662961, -0.0037143632875100),
- std::complex<double>(-0.0022414352263751, -0.0060474723525871), std::complex<double>(-0.0024377933436567, -0.0012852163337395),
- std::complex<double>(-0.6730977722052307, 0.0940030437973656), std::complex<double>(0.3711597596859299, 0.0557089394867947),
- std::complex<double>(0.2119250520015808, 0.2155514942677135), std::complex<double>(0.6727380529527980, 0.0989550572104158),
- std::complex<double>(-0.0419944347289523, 0.2355624543349744), std::complex<double>(0.1917656461134636, 0.0732470381581913),
- std::complex<double>(0.0048918921441903, 0.1588912409502319), std::complex<double>(0.0575369727210951, 0.0344677222786687),
- std::complex<double>(0.0241014578366618, 0.0547046570516960), std::complex<double>(0.0219986510834463, 0.0112189146988984),
- std::complex<double>(0.0665319393516388, -0.1418009730472832), std::complex<double>(-0.7576728614553603, -0.0472040122949963),
- std::complex<double>(-0.1017024786435272, -0.3302620837788515), std::complex<double>(-0.5600906156274197, -0.0797555201430585),
- std::complex<double>(0.0889729243872774, -0.3406964719938829), std::complex<double>(-0.1342560801672904, -0.0515926960946038),
- std::complex<double>(-0.0149335262655201, -0.2084962323582034), std::complex<double>(-0.0327252678958813, -0.0172371907472848),
- std::complex<double>(-0.0362395089905272, -0.0661322227928722), std::complex<double>(-0.0141568558526096, -0.0042676979206835),
- std::complex<double>(0.1121669548152054, 0.0504713119323919), std::complex<double>(0.1882531376700409, 0.0088411256350159),
- std::complex<double>(0.0066968933526899, 0.1181452711088882), std::complex<double>(0.0981630367567397, 0.0129921405004959),
- std::complex<double>(-0.0347327225501659, 0.1186585563636635), std::complex<double>(0.0102831315790362, 0.0046275244914932),
- std::complex<double>(0.0070209144233666, 0.0689639468490938), std::complex<double>(-0.0020239346031291, -0.0025499069613344),
- std::complex<double>(0.0132702874173192, 0.0207916487187541), std::complex<double>(0.0004387107229914, -0.0017223838914815),
- std::complex<double>(-0.0004916757488397, 0.0000266213616248), std::complex<double>(0.0006516553273188, -0.0000433166563288),
- std::complex<double>(-0.0004357897643121, 0.0000320567996700), std::complex<double>(0.0005818285824826, -0.0001021069650381),
- std::complex<double>(-0.0001047488648808, -0.0000302146563592), std::complex<double>(0.0001593350153828, -0.0000879125663990),
- std::complex<double>(-0.0000141882506567, -0.0000941521783975), std::complex<double>(-0.0000004226298134, -0.0000245060763932),
- std::complex<double>(-0.0000177429496833, -0.0000561890408003), std::complex<double>(-0.0000018388829279, 0.0000032387726477),
- std::complex<double>(0.0162495046881796, -0.0010736997976255), std::complex<double>(-0.0175635905033026, 0.0012997068962173),
- std::complex<double>(0.0138897851110661, -0.0014876219938565), std::complex<double>(-0.0150211436594772, 0.0029712291209158),
- std::complex<double>(0.0031705620225488, 0.0004838463688512), std::complex<double>(-0.0034418973689263, 0.0024603729467258),
- std::complex<double>(0.0003028387544878, 0.0026905629457281), std::complex<double>(0.0006768121359769, 0.0005901486396051),
- std::complex<double>(0.0004634797107989, 0.0016976603895716), std::complex<double>(0.0003344773954073, -0.0001499932789294),
- std::complex<double>(-0.1492097398080444, 0.0123735410547393), std::complex<double>(0.1393121453502456, -0.0121117146246749),
- std::complex<double>(-0.1217628319418324, 0.0222643129255504), std::complex<double>(0.1108579917761457, -0.0262986164183475),
- std::complex<double>(-0.0273147374272124, 0.0098595182007132), std::complex<double>(0.0208992817013466, -0.0205929453727953),
- std::complex<double>(-0.0002152227668601, -0.0089220757225133), std::complex<double>(-0.0074792188817697, -0.0043562231368076),
- std::complex<double>(-0.0012019994038721, -0.0079939660050373), std::complex<double>(-0.0035807498769946, 0.0014801422733613),
- std::complex<double>(0.1567990061437258, -0.0143275575385193), std::complex<double>(-0.1043118778001582, 0.0106756004832779),
- std::complex<double>(0.1151024257152241, -0.0225518489392044), std::complex<double>(-0.0593437249231851, 0.0216080058910987),
- std::complex<double>(0.0142781186223020, -0.0057037138045721), std::complex<double>(0.0151043140114779, 0.0141435752121475),
- std::complex<double>(-0.0057143555179676, 0.0141142700941743), std::complex<double>(0.0251435557201315, -0.0005753615445942),
- std::complex<double>(0.0004475745352473, 0.0102135659618127), std::complex<double>(0.0090474375150397, -0.0032177128650026),
- std::complex<double>(-0.0459124372023251, 0.0044990718645418), std::complex<double>(0.0135433541303599, -0.0021789296923529),
- std::complex<double>(-0.0306136798186735, 0.0064963361606382), std::complex<double>(-0.0046440676338940, -0.0037281688158807),
- std::complex<double>(-0.0006372791846825, 0.0008894047150233), std::complex<double>(-0.0181611528840412, -0.0011106177431486),
- std::complex<double>(0.0032325387394458, -0.0048123509184894), std::complex<double>(-0.0136340313457176, 0.0021185000810664),
- std::complex<double>(0.0001287985092565, -0.0032079544559908), std::complex<double>(-0.0045503800737417, 0.0015366231416036)
-};
+ std::complex<double>(0.9989322499459223, 0.0003305895124867),
+ std::complex<double>(1.0030546028872600, 0.0002157249025076),
+ std::complex<double>(0.0003002209532403, 0.0007909077657054),
+ std::complex<double>(0.0022051270911392, 0.0003834815341981),
+ std::complex<double>(-0.0003856663268042, 0.0008435910525861),
+ std::complex<double>(0.0004887765294093, 0.0002777796480946),
+ std::complex<double>(-0.0000699366665322, 0.0005136144371953),
+ std::complex<double>(0.0001520602842105, 0.0001303481681886),
+ std::complex<double>(0.0000512381993616, 0.0001550785137302),
+ std::complex<double>(0.0000819244737818, 0.0000466470412396),
+ std::complex<double>(0.0249658150445263, -0.0122024663463393),
+ std::complex<double>(-0.0917825091832822, -0.0062606338208358),
+ std::complex<double>(-0.0083709499453879, -0.0289759752488368),
+ std::complex<double>(-0.0689260153643395, -0.0111348626546314),
+ std::complex<double>(0.0116296166994115, -0.0307342946951178),
+ std::complex<double>(-0.0171249717275797, -0.0080642275561593),
+ std::complex<double>(0.0012408055399100, -0.0191295543986957),
+ std::complex<double>(-0.0051031652662961, -0.0037143632875100),
+ std::complex<double>(-0.0022414352263751, -0.0060474723525871),
+ std::complex<double>(-0.0024377933436567, -0.0012852163337395),
+ std::complex<double>(-0.6730977722052307, 0.0940030437973656),
+ std::complex<double>(0.3711597596859299, 0.0557089394867947),
+ std::complex<double>(0.2119250520015808, 0.2155514942677135),
+ std::complex<double>(0.6727380529527980, 0.0989550572104158),
+ std::complex<double>(-0.0419944347289523, 0.2355624543349744),
+ std::complex<double>(0.1917656461134636, 0.0732470381581913),
+ std::complex<double>(0.0048918921441903, 0.1588912409502319),
+ std::complex<double>(0.0575369727210951, 0.0344677222786687),
+ std::complex<double>(0.0241014578366618, 0.0547046570516960),
+ std::complex<double>(0.0219986510834463, 0.0112189146988984),
+ std::complex<double>(0.0665319393516388, -0.1418009730472832),
+ std::complex<double>(-0.7576728614553603, -0.0472040122949963),
+ std::complex<double>(-0.1017024786435272, -0.3302620837788515),
+ std::complex<double>(-0.5600906156274197, -0.0797555201430585),
+ std::complex<double>(0.0889729243872774, -0.3406964719938829),
+ std::complex<double>(-0.1342560801672904, -0.0515926960946038),
+ std::complex<double>(-0.0149335262655201, -0.2084962323582034),
+ std::complex<double>(-0.0327252678958813, -0.0172371907472848),
+ std::complex<double>(-0.0362395089905272, -0.0661322227928722),
+ std::complex<double>(-0.0141568558526096, -0.0042676979206835),
+ std::complex<double>(0.1121669548152054, 0.0504713119323919),
+ std::complex<double>(0.1882531376700409, 0.0088411256350159),
+ std::complex<double>(0.0066968933526899, 0.1181452711088882),
+ std::complex<double>(0.0981630367567397, 0.0129921405004959),
+ std::complex<double>(-0.0347327225501659, 0.1186585563636635),
+ std::complex<double>(0.0102831315790362, 0.0046275244914932),
+ std::complex<double>(0.0070209144233666, 0.0689639468490938),
+ std::complex<double>(-0.0020239346031291, -0.0025499069613344),
+ std::complex<double>(0.0132702874173192, 0.0207916487187541),
+ std::complex<double>(0.0004387107229914, -0.0017223838914815),
+ std::complex<double>(-0.0004916757488397, 0.0000266213616248),
+ std::complex<double>(0.0006516553273188, -0.0000433166563288),
+ std::complex<double>(-0.0004357897643121, 0.0000320567996700),
+ std::complex<double>(0.0005818285824826, -0.0001021069650381),
+ std::complex<double>(-0.0001047488648808, -0.0000302146563592),
+ std::complex<double>(0.0001593350153828, -0.0000879125663990),
+ std::complex<double>(-0.0000141882506567, -0.0000941521783975),
+ std::complex<double>(-0.0000004226298134, -0.0000245060763932),
+ std::complex<double>(-0.0000177429496833, -0.0000561890408003),
+ std::complex<double>(-0.0000018388829279, 0.0000032387726477),
+ std::complex<double>(0.0162495046881796, -0.0010736997976255),
+ std::complex<double>(-0.0175635905033026, 0.0012997068962173),
+ std::complex<double>(0.0138897851110661, -0.0014876219938565),
+ std::complex<double>(-0.0150211436594772, 0.0029712291209158),
+ std::complex<double>(0.0031705620225488, 0.0004838463688512),
+ std::complex<double>(-0.0034418973689263, 0.0024603729467258),
+ std::complex<double>(0.0003028387544878, 0.0026905629457281),
+ std::complex<double>(0.0006768121359769, 0.0005901486396051),
+ std::complex<double>(0.0004634797107989, 0.0016976603895716),
+ std::complex<double>(0.0003344773954073, -0.0001499932789294),
+ std::complex<double>(-0.1492097398080444, 0.0123735410547393),
+ std::complex<double>(0.1393121453502456, -0.0121117146246749),
+ std::complex<double>(-0.1217628319418324, 0.0222643129255504),
+ std::complex<double>(0.1108579917761457, -0.0262986164183475),
+ std::complex<double>(-0.0273147374272124, 0.0098595182007132),
+ std::complex<double>(0.0208992817013466, -0.0205929453727953),
+ std::complex<double>(-0.0002152227668601, -0.0089220757225133),
+ std::complex<double>(-0.0074792188817697, -0.0043562231368076),
+ std::complex<double>(-0.0012019994038721, -0.0079939660050373),
+ std::complex<double>(-0.0035807498769946, 0.0014801422733613),
+ std::complex<double>(0.1567990061437258, -0.0143275575385193),
+ std::complex<double>(-0.1043118778001582, 0.0106756004832779),
+ std::complex<double>(0.1151024257152241, -0.0225518489392044),
+ std::complex<double>(-0.0593437249231851, 0.0216080058910987),
+ std::complex<double>(0.0142781186223020, -0.0057037138045721),
+ std::complex<double>(0.0151043140114779, 0.0141435752121475),
+ std::complex<double>(-0.0057143555179676, 0.0141142700941743),
+ std::complex<double>(0.0251435557201315, -0.0005753615445942),
+ std::complex<double>(0.0004475745352473, 0.0102135659618127),
+ std::complex<double>(0.0090474375150397, -0.0032177128650026),
+ std::complex<double>(-0.0459124372023251, 0.0044990718645418),
+ std::complex<double>(0.0135433541303599, -0.0021789296923529),
+ std::complex<double>(-0.0306136798186735, 0.0064963361606382),
+ std::complex<double>(-0.0046440676338940, -0.0037281688158807),
+ std::complex<double>(-0.0006372791846825, 0.0008894047150233),
+ std::complex<double>(-0.0181611528840412, -0.0011106177431486),
+ std::complex<double>(0.0032325387394458, -0.0048123509184894),
+ std::complex<double>(-0.0136340313457176, 0.0021185000810664),
+ std::complex<double>(0.0001287985092565, -0.0032079544559908),
+ std::complex<double>(-0.0045503800737417, 0.0015366231416036)};
diff --git a/lobes/DefaultCoeffLBA.cc b/lobes/DefaultCoeffLBA.cc
index d196d2356d16a71b0f153b5364ac153bb5829f81..deb203e3a42c96ecc173ee024cc80515c229e948 100644
--- a/lobes/DefaultCoeffLBA.cc
+++ b/lobes/DefaultCoeffLBA.cc
@@ -1,5 +1,5 @@
// Beam model coefficients converted by convert_coeff.py.
-// Conversion performed on 2011/09/14/08:23:09 UTC using:
+// Conversion performed on 2011/09/14/08:23:09 UTC using:
// convert_coeff.py element_beam_LBA.coeff DefaultCoeffLBA.cc default_lba
#include <complex>
@@ -21,54 +21,103 @@ const unsigned int default_lba_coeff_shape[3] = {2, 5, 5};
// The array of coefficients in row-major order ("C"-order).
//
const std::complex<double> default_lba_coeff[100] = {
- std::complex<double>(0.9982445079290715, 0.0000650863154389), std::complex<double>(1.0006230902158257, -0.0022053287681416),
- std::complex<double>(0.0001002692200362, 0.0006838211278268), std::complex<double>(-0.0003660049052840, -0.0008418920419220),
- std::complex<double>(-0.0010581424498791, 0.0015237878543047), std::complex<double>(0.0007398729642721, 0.0028468649470433),
- std::complex<double>(-0.0039458389254656, -0.0007048354913730), std::complex<double>(0.0007040177887611, 0.0007856369612188),
- std::complex<double>(-0.0031701591723043, -0.0010521154166512), std::complex<double>(-0.0007213036752903, -0.0007227764008022),
- std::complex<double>(0.0550606068782634, 0.0011958385659938), std::complex<double>(-0.0160912944232080, 0.0703645376267940),
- std::complex<double>(0.0033849565901213, -0.0244636379385135), std::complex<double>(0.0234264238829944, 0.0084068836453700),
- std::complex<double>(0.0557107413978542, -0.0634701730653090), std::complex<double>(-0.0139549526991330, -0.1175401658864208),
- std::complex<double>(0.1336911750356096, 0.0202651327657687), std::complex<double>(-0.0113385668361727, -0.0339262369086247),
- std::complex<double>(0.0962263571740972, 0.0440074333288440), std::complex<double>(0.0313595045238824, 0.0230763038515351),
- std::complex<double>(-0.8624889445327827, -0.1522883072804402), std::complex<double>(-0.0386800869486029, -0.7569350701887934),
- std::complex<double>(0.0891332399420108, 0.1876527151756476), std::complex<double>(-0.1012363483900640, -0.1975118891151966),
- std::complex<double>(-0.6404795825927633, 0.7568775384981410), std::complex<double>(0.0767245154665722, 1.3441875993523555),
- std::complex<double>(-0.8758406699506004, 0.3350237639226141), std::complex<double>(0.2824832769101577, 0.6821307442669313),
- std::complex<double>(-0.3144282315609649, -0.2763869580286276), std::complex<double>(-0.1705959031354030, -0.0712085950559831),
- std::complex<double>(0.4039567648146965, 0.0810473144253429), std::complex<double>(-0.0350803390479135, 0.5214591717801087),
- std::complex<double>(0.2232030356124932, -0.2248154851829713), std::complex<double>(0.4704343293662089, -0.3552101485419532),
- std::complex<double>(0.9646419509627557, -0.8095088593139815), std::complex<double>(0.1635280638865702, -1.4854352979459096),
- std::complex<double>(1.0331569921006993, 0.0509705885336283), std::complex<double>(0.1501121326521990, -0.5193414816770609),
- std::complex<double>(0.4715775965513117, 0.5077361528286819), std::complex<double>(0.3847391427972284, 0.1136717951238837),
- std::complex<double>(-0.0756250564248881, 0.0056622911723172), std::complex<double>(-0.1267444401630109, -0.0349676272376008),
- std::complex<double>(-0.1793752883639813, 0.0720222655359702), std::complex<double>(-0.2678542619793421, 0.3152115802895427),
- std::complex<double>(-0.3718069213271066, 0.2275266747872172), std::complex<double>(-0.1372223722572021, 0.4314989948093362),
- std::complex<double>(-0.3316657641578328, -0.1655909947939444), std::complex<double>(-0.2158100484836540, 0.0614504774034524),
- std::complex<double>(-0.1901597954359592, -0.2294955549701665), std::complex<double>(-0.1864961465389693, -0.0486276177310768),
- std::complex<double>(-0.0000762326746410, 0.0000118155774181), std::complex<double>(0.0000118903581604, -0.0000251324432498),
- std::complex<double>(-0.0002204197663391, -0.0000213776348027), std::complex<double>(0.0001477083861977, 0.0000599750510518),
- std::complex<double>(-0.0003281057522772, -0.0000770207588466), std::complex<double>(0.0003478997686964, 0.0001481982639746),
- std::complex<double>(-0.0000625695757282, 0.0000249138990722), std::complex<double>(-0.0000960097542525, 0.0002521364065803),
- std::complex<double>(0.0001275344578325, 0.0000652362392482), std::complex<double>(-0.0003113309221942, 0.0001956734476566),
- std::complex<double>(0.0029807707669629, -0.0003262084082071), std::complex<double>(0.0001639620574332, -0.0000266272685197),
- std::complex<double>(0.0076282580587895, 0.0026614359017468), std::complex<double>(-0.0044850263974801, -0.0058337192660638),
- std::complex<double>(0.0124258438959177, 0.0067985224235178), std::complex<double>(-0.0126349778957970, -0.0100656881493938),
- std::complex<double>(0.0059031372522229, 0.0008660479915339), std::complex<double>(0.0039660364524413, -0.0100356333791398),
- std::complex<double>(-0.0020520685193773, -0.0028564379463666), std::complex<double>(0.0121039958869239, -0.0059701468961263),
- std::complex<double>(-0.0229975846564195, 0.0010565261888195), std::complex<double>(0.0019573207027441, 0.0050550600926414),
- std::complex<double>(-0.0682274156850413, -0.0758159820140411), std::complex<double>(0.0497303968865466, 0.1019681987654797),
- std::complex<double>(-0.1757936183439326, -0.1363710820472197), std::complex<double>(0.1765450269056824, 0.1555919358121995),
- std::complex<double>(-0.1541299429420569, -0.0281422177614844), std::complex<double>(0.0816399676454817, 0.0691599035109852),
- std::complex<double>(-0.0235110916473515, 0.0306385386726702), std::complex<double>(-0.0474273292450285, 0.0116831908947225),
- std::complex<double>(0.0333560984394624, -0.0009767086536162), std::complex<double>(0.0141704479374002, -0.0205386534626779),
- std::complex<double>(0.0562541280098909, 0.0743149092143081), std::complex<double>(-0.0226634801339250, -0.1439026188572270),
- std::complex<double>(0.1238595124159999, 0.1766108700786397), std::complex<double>(-0.1307647072780430, -0.2090615438301942),
- std::complex<double>(0.1557916917691289, 0.0646351862895731), std::complex<double>(0.0170294191358757, -0.1027926845803498),
- std::complex<double>(0.0543537332385954, -0.0366524906364179), std::complex<double>(0.1127180664279469, -0.0176607923511174),
- std::complex<double>(-0.0126732549319889, 0.0002042370658763), std::complex<double>(-0.0101360135082899, 0.0114084024114141),
- std::complex<double>(-0.0102147881225462, -0.0176848554302252), std::complex<double>(-0.0051268936720694, 0.0527621533959941),
- std::complex<double>(-0.0110701836450407, -0.0593085026046026), std::complex<double>(0.0140598301629874, 0.0738668439833535),
- std::complex<double>(-0.0389912915621699, -0.0301364165752433), std::complex<double>(-0.0462331759359031, 0.0405864871628086),
- std::complex<double>(-0.0251598701859194, 0.0115712688652445), std::complex<double>(-0.0563476280247398, 0.0079787883434624)
-};
+ std::complex<double>(0.9982445079290715, 0.0000650863154389),
+ std::complex<double>(1.0006230902158257, -0.0022053287681416),
+ std::complex<double>(0.0001002692200362, 0.0006838211278268),
+ std::complex<double>(-0.0003660049052840, -0.0008418920419220),
+ std::complex<double>(-0.0010581424498791, 0.0015237878543047),
+ std::complex<double>(0.0007398729642721, 0.0028468649470433),
+ std::complex<double>(-0.0039458389254656, -0.0007048354913730),
+ std::complex<double>(0.0007040177887611, 0.0007856369612188),
+ std::complex<double>(-0.0031701591723043, -0.0010521154166512),
+ std::complex<double>(-0.0007213036752903, -0.0007227764008022),
+ std::complex<double>(0.0550606068782634, 0.0011958385659938),
+ std::complex<double>(-0.0160912944232080, 0.0703645376267940),
+ std::complex<double>(0.0033849565901213, -0.0244636379385135),
+ std::complex<double>(0.0234264238829944, 0.0084068836453700),
+ std::complex<double>(0.0557107413978542, -0.0634701730653090),
+ std::complex<double>(-0.0139549526991330, -0.1175401658864208),
+ std::complex<double>(0.1336911750356096, 0.0202651327657687),
+ std::complex<double>(-0.0113385668361727, -0.0339262369086247),
+ std::complex<double>(0.0962263571740972, 0.0440074333288440),
+ std::complex<double>(0.0313595045238824, 0.0230763038515351),
+ std::complex<double>(-0.8624889445327827, -0.1522883072804402),
+ std::complex<double>(-0.0386800869486029, -0.7569350701887934),
+ std::complex<double>(0.0891332399420108, 0.1876527151756476),
+ std::complex<double>(-0.1012363483900640, -0.1975118891151966),
+ std::complex<double>(-0.6404795825927633, 0.7568775384981410),
+ std::complex<double>(0.0767245154665722, 1.3441875993523555),
+ std::complex<double>(-0.8758406699506004, 0.3350237639226141),
+ std::complex<double>(0.2824832769101577, 0.6821307442669313),
+ std::complex<double>(-0.3144282315609649, -0.2763869580286276),
+ std::complex<double>(-0.1705959031354030, -0.0712085950559831),
+ std::complex<double>(0.4039567648146965, 0.0810473144253429),
+ std::complex<double>(-0.0350803390479135, 0.5214591717801087),
+ std::complex<double>(0.2232030356124932, -0.2248154851829713),
+ std::complex<double>(0.4704343293662089, -0.3552101485419532),
+ std::complex<double>(0.9646419509627557, -0.8095088593139815),
+ std::complex<double>(0.1635280638865702, -1.4854352979459096),
+ std::complex<double>(1.0331569921006993, 0.0509705885336283),
+ std::complex<double>(0.1501121326521990, -0.5193414816770609),
+ std::complex<double>(0.4715775965513117, 0.5077361528286819),
+ std::complex<double>(0.3847391427972284, 0.1136717951238837),
+ std::complex<double>(-0.0756250564248881, 0.0056622911723172),
+ std::complex<double>(-0.1267444401630109, -0.0349676272376008),
+ std::complex<double>(-0.1793752883639813, 0.0720222655359702),
+ std::complex<double>(-0.2678542619793421, 0.3152115802895427),
+ std::complex<double>(-0.3718069213271066, 0.2275266747872172),
+ std::complex<double>(-0.1372223722572021, 0.4314989948093362),
+ std::complex<double>(-0.3316657641578328, -0.1655909947939444),
+ std::complex<double>(-0.2158100484836540, 0.0614504774034524),
+ std::complex<double>(-0.1901597954359592, -0.2294955549701665),
+ std::complex<double>(-0.1864961465389693, -0.0486276177310768),
+ std::complex<double>(-0.0000762326746410, 0.0000118155774181),
+ std::complex<double>(0.0000118903581604, -0.0000251324432498),
+ std::complex<double>(-0.0002204197663391, -0.0000213776348027),
+ std::complex<double>(0.0001477083861977, 0.0000599750510518),
+ std::complex<double>(-0.0003281057522772, -0.0000770207588466),
+ std::complex<double>(0.0003478997686964, 0.0001481982639746),
+ std::complex<double>(-0.0000625695757282, 0.0000249138990722),
+ std::complex<double>(-0.0000960097542525, 0.0002521364065803),
+ std::complex<double>(0.0001275344578325, 0.0000652362392482),
+ std::complex<double>(-0.0003113309221942, 0.0001956734476566),
+ std::complex<double>(0.0029807707669629, -0.0003262084082071),
+ std::complex<double>(0.0001639620574332, -0.0000266272685197),
+ std::complex<double>(0.0076282580587895, 0.0026614359017468),
+ std::complex<double>(-0.0044850263974801, -0.0058337192660638),
+ std::complex<double>(0.0124258438959177, 0.0067985224235178),
+ std::complex<double>(-0.0126349778957970, -0.0100656881493938),
+ std::complex<double>(0.0059031372522229, 0.0008660479915339),
+ std::complex<double>(0.0039660364524413, -0.0100356333791398),
+ std::complex<double>(-0.0020520685193773, -0.0028564379463666),
+ std::complex<double>(0.0121039958869239, -0.0059701468961263),
+ std::complex<double>(-0.0229975846564195, 0.0010565261888195),
+ std::complex<double>(0.0019573207027441, 0.0050550600926414),
+ std::complex<double>(-0.0682274156850413, -0.0758159820140411),
+ std::complex<double>(0.0497303968865466, 0.1019681987654797),
+ std::complex<double>(-0.1757936183439326, -0.1363710820472197),
+ std::complex<double>(0.1765450269056824, 0.1555919358121995),
+ std::complex<double>(-0.1541299429420569, -0.0281422177614844),
+ std::complex<double>(0.0816399676454817, 0.0691599035109852),
+ std::complex<double>(-0.0235110916473515, 0.0306385386726702),
+ std::complex<double>(-0.0474273292450285, 0.0116831908947225),
+ std::complex<double>(0.0333560984394624, -0.0009767086536162),
+ std::complex<double>(0.0141704479374002, -0.0205386534626779),
+ std::complex<double>(0.0562541280098909, 0.0743149092143081),
+ std::complex<double>(-0.0226634801339250, -0.1439026188572270),
+ std::complex<double>(0.1238595124159999, 0.1766108700786397),
+ std::complex<double>(-0.1307647072780430, -0.2090615438301942),
+ std::complex<double>(0.1557916917691289, 0.0646351862895731),
+ std::complex<double>(0.0170294191358757, -0.1027926845803498),
+ std::complex<double>(0.0543537332385954, -0.0366524906364179),
+ std::complex<double>(0.1127180664279469, -0.0176607923511174),
+ std::complex<double>(-0.0126732549319889, 0.0002042370658763),
+ std::complex<double>(-0.0101360135082899, 0.0114084024114141),
+ std::complex<double>(-0.0102147881225462, -0.0176848554302252),
+ std::complex<double>(-0.0051268936720694, 0.0527621533959941),
+ std::complex<double>(-0.0110701836450407, -0.0593085026046026),
+ std::complex<double>(0.0140598301629874, 0.0738668439833535),
+ std::complex<double>(-0.0389912915621699, -0.0301364165752433),
+ std::complex<double>(-0.0462331759359031, 0.0405864871628086),
+ std::complex<double>(-0.0251598701859194, 0.0115712688652445),
+ std::complex<double>(-0.0563476280247398, 0.0079787883434624)};
diff --git a/lobes/ElementResponse.cc b/lobes/ElementResponse.cc
index b2ff653464059aaf4e19a8d5b3b9e8003fa6af18..9be03164213b9986b6eccb44d83c00bedb41c164 100644
--- a/lobes/ElementResponse.cc
+++ b/lobes/ElementResponse.cc
@@ -27,128 +27,121 @@
// The coefficients are kept in an unnamed namespace which effectively makes
// them invisible outside this translation unit.
-namespace
-{
+namespace {
#include "DefaultCoeffLBA.cc"
#include "DefaultCoeffHBA.cc"
-}
+} // namespace
-namespace LOFAR
-{
+namespace LOFAR {
void element_response_lba(double freq, double theta, double phi,
- std::complex<double> (&response)[2][2])
-{
- element_response(freq, theta, phi, response, default_lba_freq_center,
- default_lba_freq_range, default_lba_coeff_shape, default_lba_coeff);
+ std::complex<double> (&response)[2][2]) {
+ element_response(freq, theta, phi, response, default_lba_freq_center,
+ default_lba_freq_range, default_lba_coeff_shape,
+ default_lba_coeff);
}
void element_response_hba(double freq, double theta, double phi,
- std::complex<double> (&response)[2][2])
-{
- element_response(freq, theta, phi, response, default_hba_freq_center,
- default_hba_freq_range, default_hba_coeff_shape, default_hba_coeff);
+ std::complex<double> (&response)[2][2]) {
+ element_response(freq, theta, phi, response, default_hba_freq_center,
+ default_hba_freq_range, default_hba_coeff_shape,
+ default_hba_coeff);
}
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[])
-{
- // Initialize the response to zero.
- response[0][0] = 0.0;
- response[0][1] = 0.0;
- response[1][0] = 0.0;
- response[1][1] = 0.0;
-
- // Clip directions below the horizon.
- if(theta >= M_PI_2)
- {
- return;
+ std::complex<double> (&response)[2][2],
+ double freq_center, double freq_range,
+ const unsigned int (&coeff_shape)[3],
+ const std::complex<double> coeff[]) {
+ // Initialize the response to zero.
+ response[0][0] = 0.0;
+ response[0][1] = 0.0;
+ response[1][0] = 0.0;
+ response[1][1] = 0.0;
+
+ // Clip directions below the horizon.
+ if (theta >= M_PI_2) {
+ return;
+ }
+
+ const unsigned int nHarmonics = coeff_shape[0];
+ const unsigned int nPowerTheta = coeff_shape[1];
+ const unsigned int nPowerFreq = coeff_shape[2];
+
+ // The model is parameterized in terms of a normalized frequency in the
+ // range [-1, 1]. The appropriate conversion is taken care of below.
+ freq = (freq - freq_center) / freq_range;
+
+ // The variables sign and kappa are used to compute the value of kappa
+ // mentioned in the description of the beam model [kappa = (-1)^k * (2 * k
+ //+ 1)] incrementally.
+ int sign = 1, kappa = 1;
+
+ std::complex<double> P[2], Pj[2];
+ for (unsigned int k = 0; k < nHarmonics; ++k) {
+ // Compute the (diagonal) projection matrix P for the current harmonic.
+ // This requires the evaluation of two polynomials in theta and freq (of
+ // degree nPowerTheta in theta and nPowerFreq in freq), one for each
+ // element of P. The polynomials are evaluated using Horner's rule.
+
+ // Horner's rule requires backward iteration of the coefficients, so
+ // move the iterator to the first location past the end of the block of
+ // coefficients for the current harmonic (k).
+ coeff += nPowerTheta * nPowerFreq * 2;
+
+ // Evaluate the polynomial. Note that the iterator is always decremented
+ // before it is dereferenced, using the prefix operator--. After
+ // evaluation of the polynomial, the iterator points exactly to the
+ // beginning of the block of coefficients for the current harmonic (k),
+ // that is, all the decrements together exactly cancel the increment
+ // aplied above.
+
+ // Evaluate the highest order term. Note that the order of the
+ // assigments is important because of the iterator decrement, i.e. P[1]
+ // should be assigned first.
+ P[1] = *--coeff;
+ P[0] = *--coeff;
+
+ for (unsigned int i = 0; i < nPowerFreq - 1; ++i) {
+ P[1] = P[1] * freq + *--coeff;
+ P[0] = P[0] * freq + *--coeff;
}
- const unsigned int nHarmonics = coeff_shape[0];
- const unsigned int nPowerTheta = coeff_shape[1];
- const unsigned int nPowerFreq = coeff_shape[2];
-
- // The model is parameterized in terms of a normalized frequency in the
- // range [-1, 1]. The appropriate conversion is taken care of below.
- freq = (freq - freq_center) / freq_range;
-
- // The variables sign and kappa are used to compute the value of kappa
- // mentioned in the description of the beam model [kappa = (-1)^k * (2 * k
- //+ 1)] incrementally.
- int sign = 1, kappa = 1;
-
- std::complex<double> P[2], Pj[2];
- for(unsigned int k = 0; k < nHarmonics; ++k)
- {
- // Compute the (diagonal) projection matrix P for the current harmonic.
- // This requires the evaluation of two polynomials in theta and freq (of
- // degree nPowerTheta in theta and nPowerFreq in freq), one for each
- // element of P. The polynomials are evaluated using Horner's rule.
-
- // Horner's rule requires backward iteration of the coefficients, so
- // move the iterator to the first location past the end of the block of
- // coefficients for the current harmonic (k).
- coeff += nPowerTheta * nPowerFreq * 2;
-
- // Evaluate the polynomial. Note that the iterator is always decremented
- // before it is dereferenced, using the prefix operator--. After
- // evaluation of the polynomial, the iterator points exactly to the
- // beginning of the block of coefficients for the current harmonic (k),
- // that is, all the decrements together exactly cancel the increment
- // aplied above.
-
- // Evaluate the highest order term. Note that the order of the
- // assigments is important because of the iterator decrement, i.e. P[1]
- // should be assigned first.
- P[1] = *--coeff;
- P[0] = *--coeff;
-
- for(unsigned int i = 0; i < nPowerFreq - 1; ++i)
- {
- P[1] = P[1] * freq + *--coeff;
- P[0] = P[0] * freq + *--coeff;
- }
-
- // Evaluate the remaining terms.
- for(unsigned int j = 0; j < nPowerTheta - 1; ++j)
- {
- Pj[1] = *--coeff;
- Pj[0] = *--coeff;
-
- for(unsigned int i = 0; i < nPowerFreq - 1; ++i)
- {
- Pj[1] = Pj[1] * freq + *--coeff;
- Pj[0] = Pj[0] * freq + *--coeff;
- }
-
- P[1] = P[1] * theta + Pj[1];
- P[0] = P[0] * theta + Pj[0];
- }
-
- // Because the increment and decrements cancel, the iterator points to
- // the same location as at the beginning of this iteration of the outer
- // loop. The next iteration should use the coefficients for the next
- // harmonic (k), so we move the iterator to the start of that block.
- coeff += nPowerTheta * nPowerFreq * 2;
-
- // Compute the Jones matrix for the current harmonic, by rotating P over
- // kappa * az, and add it to the result.
- const double angle = sign * kappa * phi;
- const double caz = std::cos(angle);
- const double saz = std::sin(angle);
-
- response[0][0] += caz * P[0];
- response[0][1] += -saz * P[1];
- response[1][0] += saz * P[0];
- response[1][1] += caz * P[1];
-
- // Update sign and kappa.
- sign = -sign;
- kappa += 2;
+ // Evaluate the remaining terms.
+ for (unsigned int j = 0; j < nPowerTheta - 1; ++j) {
+ Pj[1] = *--coeff;
+ Pj[0] = *--coeff;
+
+ for (unsigned int i = 0; i < nPowerFreq - 1; ++i) {
+ Pj[1] = Pj[1] * freq + *--coeff;
+ Pj[0] = Pj[0] * freq + *--coeff;
+ }
+
+ P[1] = P[1] * theta + Pj[1];
+ P[0] = P[0] * theta + Pj[0];
}
+
+ // Because the increment and decrements cancel, the iterator points to
+ // the same location as at the beginning of this iteration of the outer
+ // loop. The next iteration should use the coefficients for the next
+ // harmonic (k), so we move the iterator to the start of that block.
+ coeff += nPowerTheta * nPowerFreq * 2;
+
+ // Compute the Jones matrix for the current harmonic, by rotating P over
+ // kappa * az, and add it to the result.
+ const double angle = sign * kappa * phi;
+ const double caz = std::cos(angle);
+ const double saz = std::sin(angle);
+
+ response[0][0] += caz * P[0];
+ response[0][1] += -saz * P[1];
+ response[1][0] += saz * P[0];
+ response[1][1] += caz * P[1];
+
+ // Update sign and kappa.
+ sign = -sign;
+ kappa += 2;
+ }
}
-} //# namespace LOFAR
+} // namespace LOFAR
diff --git a/lobes/LOBESElementResponse.cc b/lobes/LOBESElementResponse.cc
index d54b18ed2e83f8f166592619c5d399fce8442199..a475f56929b0f8e28886fe4fab1cda7643e79423 100644
--- a/lobes/LOBESElementResponse.cc
+++ b/lobes/LOBESElementResponse.cc
@@ -5,30 +5,24 @@
namespace LOFAR {
namespace StationResponse {
-LOBESElementResponse::LOBESElementResponse(std::string name)
-{
-}
+LOBESElementResponse::LOBESElementResponse(std::string name) {}
void LOBESElementResponse::response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const
-{
-}
-
-std::shared_ptr<LOBESElementResponse> LOBESElementResponse::getInstance(std::string name)
-{
- static std::map<std::string, std::shared_ptr<LOBESElementResponse>> name_response_map;
-
- auto entry = name_response_map.find(name);
- if (entry == name_response_map.end()) {
- entry = name_response_map.insert(entry, {name, std::make_shared<LOBESElementResponse>(name)});
- }
- return entry->second;
-
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const {}
+
+std::shared_ptr<LOBESElementResponse> LOBESElementResponse::getInstance(
+ std::string name) {
+ static std::map<std::string, std::shared_ptr<LOBESElementResponse>>
+ name_response_map;
+
+ auto entry = name_response_map.find(name);
+ if (entry == name_response_map.end()) {
+ entry = name_response_map.insert(
+ entry, {name, std::make_shared<LOBESElementResponse>(name)});
+ }
+ return entry->second;
};
-} // namespace StationResponse
-} // namespace LOFAR
-
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/lobes/lobes.cc b/lobes/lobes.cc
index 01664be4dd6c64393ec282a82682ebfb77b4dd21..ff5db10e9bfa9145cde79b2b1ad05f8b5ba07385 100644
--- a/lobes/lobes.cc
+++ b/lobes/lobes.cc
@@ -12,282 +12,281 @@
using namespace std::complex_literals;
-
#include <tuple>
// python style enumerate function
// To make it possible to write:
// for (auto [i, v] : enumerate(iterable)) {...}
-template <typename T,
- typename TIter = decltype(std::begin(std::declval<T>())),
+template <typename T, typename TIter = decltype(std::begin(std::declval<T>())),
typename = decltype(std::end(std::declval<T>()))>
-constexpr auto enumerate(T && iterable)
-{
- struct iterator
- {
- size_t i;
- TIter iter;
- bool operator != (const iterator & other) const { return iter != other.iter; }
- void operator ++ () { ++i; ++iter; }
- auto operator * () const { return std::tie(i, *iter); }
- };
- struct iterable_wrapper
- {
- T iterable;
- auto begin() { return iterator{ 0, std::begin(iterable) }; }
- auto end() { return iterator{ 0, std::end(iterable) }; }
- };
- return iterable_wrapper{ std::forward<T>(iterable) };
+constexpr auto enumerate(T &&iterable) {
+ struct iterator {
+ size_t i;
+ TIter iter;
+ bool operator!=(const iterator &other) const { return iter != other.iter; }
+ void operator++() {
+ ++i;
+ ++iter;
+ }
+ auto operator*() const { return std::tie(i, *iter); }
+ };
+ struct iterable_wrapper {
+ T iterable;
+ auto begin() { return iterator{0, std::begin(iterable)}; }
+ auto end() { return iterator{0, std::end(iterable)}; }
+ };
+ return iterable_wrapper{std::forward<T>(iterable)};
}
+double P20(double x) { return 0.5 * (3 * x * x - 1); }
+double P21(double x) { return 3.0 * x * std::sqrt(1 - x * x); }
+double P22(double x) { return 3 * (1 - x * x); }
-double P20(double x) { return 0.5*(3*x*x-1); }
-double P21(double x) { return 3.0*x*std::sqrt(1-x*x); }
-double P22(double x) { return 3*(1-x*x); }
-
-int plustwo(int a)
-{
- std::cout << std::assoc_legendre(2, 0, 0.5) << '=' << P20(0.5) << '\n'
- << std::assoc_legendre(2, 1, 0.5) << '=' << P21(0.5) << '\n'
- << std::assoc_legendre(2, 2, 0.5) << '=' << P22(0.5) << '\n';
+int plustwo(int a) {
+ std::cout << std::assoc_legendre(2, 0, 0.5) << '=' << P20(0.5) << '\n'
+ << std::assoc_legendre(2, 1, 0.5) << '=' << P21(0.5) << '\n'
+ << std::assoc_legendre(2, 2, 0.5) << '=' << P22(0.5) << '\n';
- return a+2;
+ return a + 2;
}
+Eigen::MatrixXd big_mat() { return Eigen::MatrixXd::Zero(10000, 10000); }
-Eigen::MatrixXd big_mat()
-{
- return Eigen::MatrixXd::Zero(10000, 10000);
-}
-
+Eigen::ArrayXd P(int m, int n, py::EigenDRef<Eigen::ArrayXd> x) {
+ auto N = x.rows();
-Eigen::ArrayXd P(int m, int n, py::EigenDRef<Eigen::ArrayXd> x)
-{
- auto N = x.rows();
+ Eigen::ArrayXd result(N); // = Eigen::VectorXd::Zero(10, 1);
+ for (int i = 0; i < N; i++) {
+ result[i] = std::assoc_legendre(n, std::abs(m), x[i]);
+ }
- Eigen::ArrayXd result(N); // = Eigen::VectorXd::Zero(10, 1);
- for(int i = 0; i<N; i++)
- {
- result[i] = std::assoc_legendre(n, std::abs(m), x[i]);
- }
+ if (m < 0) {
+ result *=
+ std::pow(-1, -m) * std::tgamma(n + m + 1) / std::tgamma(n - m + 1);
+ }
- if (m<0)
- {
- result *= std::pow(-1, -m) * std::tgamma(n+m+1) / std::tgamma(n-m+1);
- }
-
- return std::move(result);
+ return std::move(result);
}
-Eigen::ArrayXd Pacc(int m, int n, py::EigenDRef<Eigen::ArrayXd> x)
-{
- auto N = x.rows();
+Eigen::ArrayXd Pacc(int m, int n, py::EigenDRef<Eigen::ArrayXd> x) {
+ auto N = x.rows();
- return (-(n+m) * (n-m+1.0) * Eigen::sqrt(1.0 - x*x) * P(m-1, n, x) - m * x * P(m, n, x))/(x*x-1.0);
+ return (-(n + m) * (n - m + 1.0) * Eigen::sqrt(1.0 - x * x) * P(m - 1, n, x) -
+ m * x * P(m, n, x)) /
+ (x * x - 1.0);
}
-
/** Compute Spherical Wave Function
- *
- * returns std::pair(Etheta, Ephi)
- */
-std::pair<Eigen::VectorXcd,Eigen::VectorXcd> F4far_new(int s, int m, int n, py::EigenDRef<const Eigen::ArrayXd> theta, py::EigenDRef<const Eigen::ArrayXd> phi, double beta)
-{
- int N = theta.rows();
- Eigen::VectorXd result(N);
-
- double C;
- if (m) {
- C = beta * std::sqrt(60.0) * 1.0/std::sqrt(n*(n+1.0)) * std::pow(-m/std::abs(m), m);
- }
- else {
- C = beta * std::sqrt(60.0) * 1.0/std::sqrt(n * (n+1.0));
- }
-
- Eigen::ArrayXcd q2;
- Eigen::ArrayXcd q3;
-
-
- Eigen::ArrayXd cos_theta = Eigen::cos(theta);
- Eigen::ArrayXd sin_theta = Eigen::sin(theta);
- Eigen::ArrayXd P_cos_theta = P(std::abs(m),n,cos_theta);
- Eigen::ArrayXd Pacc_cos_theta = Pacc(std::abs(m),n,cos_theta);
- Eigen::ArrayXcd exp_i_m_phi = Eigen::exp(1.0i*double(m)*phi);
-
- if (s == 1) {
-
- q2 = C * std::pow(-1.0i, -n-1)/beta * 1.0i * double(m) / (sin_theta) *
- std::sqrt((2.*n+1)/2.0 * std::tgamma(n-std::abs(m)+1) / std::tgamma(n+std::abs(m)+1)) *
- P_cos_theta * exp_i_m_phi;
-
- q3 = C * std::pow(-1.0i, -n-1)/beta * std::sqrt((2.*n+1) / 2.0 * std::tgamma(n-abs(m)+1) / std::tgamma(n+abs(m)+1)) *
- Pacc_cos_theta * sin_theta * exp_i_m_phi;
- }
- else if (s==2) {
- q2 = -C * std::pow(-1.0i, -n) / beta * std::sqrt((2.*n+1)/2.0 * std::tgamma(n-abs(m)+1) / std::tgamma(n+abs(m)+1)) *
- Pacc_cos_theta * sin_theta * exp_i_m_phi;
-
- q3 = C * std::pow(-1.0i, -n) / beta * 1.0i * double(m)/sin_theta * std::sqrt((2.*n+1) / 2.0 * std::tgamma(n-abs(m)+1) / std::tgamma(n+abs(m)+1)) *
- P_cos_theta * exp_i_m_phi;
- }
-
- return std::make_pair(q2,q3);
+ *
+ * returns std::pair(Etheta, Ephi)
+ */
+std::pair<Eigen::VectorXcd, Eigen::VectorXcd> F4far_new(
+ int s, int m, int n, py::EigenDRef<const Eigen::ArrayXd> theta,
+ py::EigenDRef<const Eigen::ArrayXd> phi, double beta) {
+ int N = theta.rows();
+ Eigen::VectorXd result(N);
+
+ double C;
+ if (m) {
+ C = beta * std::sqrt(60.0) * 1.0 / std::sqrt(n * (n + 1.0)) *
+ std::pow(-m / std::abs(m), m);
+ } else {
+ C = beta * std::sqrt(60.0) * 1.0 / std::sqrt(n * (n + 1.0));
+ }
+
+ Eigen::ArrayXcd q2;
+ Eigen::ArrayXcd q3;
+
+ Eigen::ArrayXd cos_theta = Eigen::cos(theta);
+ Eigen::ArrayXd sin_theta = Eigen::sin(theta);
+ Eigen::ArrayXd P_cos_theta = P(std::abs(m), n, cos_theta);
+ Eigen::ArrayXd Pacc_cos_theta = Pacc(std::abs(m), n, cos_theta);
+ Eigen::ArrayXcd exp_i_m_phi = Eigen::exp(1.0i * double(m) * phi);
+
+ if (s == 1) {
+
+ q2 = C * std::pow(-1.0i, -n - 1) / beta * 1.0i * double(m) /
+ (sin_theta)*std::sqrt((2. * n + 1) / 2.0 *
+ std::tgamma(n - std::abs(m) + 1) /
+ std::tgamma(n + std::abs(m) + 1)) *
+ P_cos_theta * exp_i_m_phi;
+
+ q3 = C * std::pow(-1.0i, -n - 1) / beta *
+ std::sqrt((2. * n + 1) / 2.0 * std::tgamma(n - abs(m) + 1) /
+ std::tgamma(n + abs(m) + 1)) *
+ Pacc_cos_theta * sin_theta * exp_i_m_phi;
+ } else if (s == 2) {
+ q2 = -C * std::pow(-1.0i, -n) / beta *
+ std::sqrt((2. * n + 1) / 2.0 * std::tgamma(n - abs(m) + 1) /
+ std::tgamma(n + abs(m) + 1)) *
+ Pacc_cos_theta * sin_theta * exp_i_m_phi;
+
+ q3 = C * std::pow(-1.0i, -n) / beta * 1.0i * double(m) / sin_theta *
+ std::sqrt((2. * n + 1) / 2.0 * std::tgamma(n - abs(m) + 1) /
+ std::tgamma(n + abs(m) + 1)) *
+ P_cos_theta * exp_i_m_phi;
+ }
+
+ return std::make_pair(q2, q3);
}
-Eigen::Array<std::complex<double>, 2, 2> element_response_lba(double freq, double theta, double phi)
-{
- Eigen::Array<std::complex<double>, 2, 2> response;
+Eigen::Array<std::complex<double>, 2, 2> element_response_lba(double freq,
+ double theta,
+ double phi) {
+ Eigen::Array<std::complex<double>, 2, 2> response;
- typedef std::complex<double> Response[2][2];
- Response &response_ = *((Response *) response.data());
+ typedef std::complex<double> Response[2][2];
+ Response &response_ = *((Response *)response.data());
- LOFAR::element_response_lba(freq, theta, phi, response_);
- return std::move(response);
+ LOFAR::element_response_lba(freq, theta, phi, response_);
+ return std::move(response);
}
-
-LobesBeamModel::LobesBeamModel(const std::string &data_file_name)
-{
- H5::H5File h5file( data_file_name.c_str(), H5F_ACC_RDONLY );
-
- H5::DataSet dataset = h5file.openDataSet( "coefficients" );
- H5::DataSpace dataspace = dataset.getSpace();
- int nr_elements = dataspace.getSimpleExtentNpoints();
-
- const std::string REAL( "r" );
- const std::string IMAG( "i" );
-
- H5::CompType h5_dcomplex( sizeof(std::complex<double>) );
- h5_dcomplex.insertMember( REAL, 0, H5::PredType::NATIVE_DOUBLE);
- h5_dcomplex.insertMember( IMAG, sizeof(double), H5::PredType::NATIVE_DOUBLE);
-
- /*
- * Get the number of dimensions in the dataspace.
- */
- int ndims_coefficients = dataspace.getSimpleExtentNdims();
- /*
- * Get the dimension size of each dimension in the dataspace and
- * display them.
- */
- hsize_t dims_coefficients[ndims_coefficients];
- dataspace.getSimpleExtentDims( dims_coefficients, NULL);
-
- std::cout << "Dimensions of coefficients:" << std::endl << " ";
- for(int i = 0; i<ndims_coefficients; i++)
- {
- std::cout << (size_t)(dims_coefficients[i]);
- if (i < ndims_coefficients-1)
- std::cout << " x ";
- else
- std::cout << std::endl;
- }
-
- // Reshape coefficients
- size_t nr_elements_first_dims = 1;
- for(int i = 0; i<ndims_coefficients-1; i++) nr_elements_first_dims *= dims_coefficients[i];
- std::cout << nr_elements_first_dims << " x " << (size_t)(dims_coefficients[ndims_coefficients-1]) << std::endl;
-
- m_coefficients = Eigen::ArrayXXcd(nr_elements_first_dims, (size_t)dims_coefficients[ndims_coefficients-1]);
- m_coefficients_shape = std::vector<size_t>(dims_coefficients, dims_coefficients + ndims_coefficients);
- dataset.read(m_coefficients.data(), h5_dcomplex);
-
-//=====================================
- dataset = h5file.openDataSet( "frequencies" );
- dataspace = dataset.getSpace();
- nr_elements = dataspace.getSimpleExtentNpoints();
-
- m_frequencies = std::vector<double> (nr_elements);
-
- dataset.read(m_frequencies.data(), H5::PredType::NATIVE_DOUBLE);
- std::cout << "Frequencies:" << std::endl;
- for(auto freq : m_frequencies) {
- std::cout << freq << " ";
- }
- std::cout << std::endl << std::endl;
-
-
- dataset = h5file.openDataSet( "nms" );
- dataspace = dataset.getSpace();
- nr_elements = dataspace.getSimpleExtentNpoints();
-
-
- /*
- * Get the number of dimensions in the dataspace.
- */
- int ndims_nms = dataspace.getSimpleExtentNdims();
- /*
- * Get the dimension size of each dimension in the dataspace and
- * display them.
- */
- hsize_t dims_nms[ndims_nms];
- dataspace.getSimpleExtentDims( dims_nms, NULL);
-
- std::cout << "Dimensions of nms:" << ndims_nms << std::endl << " ";
-
- for(int i = 0; i<ndims_nms; i++)
- {
- std::cout << (size_t)(dims_nms[i]);
- if (i < ndims_nms-1)
- std::cout << " x ";
- }
- std::cout << std::endl;
-
- m_nms = py::array_t<int>({dims_nms[0], dims_nms[1]});
- dataset.read(m_nms.mutable_data(), H5::PredType::NATIVE_INT);
+LobesBeamModel::LobesBeamModel(const std::string &data_file_name) {
+ H5::H5File h5file(data_file_name.c_str(), H5F_ACC_RDONLY);
+
+ H5::DataSet dataset = h5file.openDataSet("coefficients");
+ H5::DataSpace dataspace = dataset.getSpace();
+ int nr_elements = dataspace.getSimpleExtentNpoints();
+
+ const std::string REAL("r");
+ const std::string IMAG("i");
+
+ H5::CompType h5_dcomplex(sizeof(std::complex<double>));
+ h5_dcomplex.insertMember(REAL, 0, H5::PredType::NATIVE_DOUBLE);
+ h5_dcomplex.insertMember(IMAG, sizeof(double), H5::PredType::NATIVE_DOUBLE);
+
+ /*
+ * Get the number of dimensions in the dataspace.
+ */
+ int ndims_coefficients = dataspace.getSimpleExtentNdims();
+ /*
+ * Get the dimension size of each dimension in the dataspace and
+ * display them.
+ */
+ hsize_t dims_coefficients[ndims_coefficients];
+ dataspace.getSimpleExtentDims(dims_coefficients, NULL);
+
+ std::cout << "Dimensions of coefficients:" << std::endl << " ";
+ for (int i = 0; i < ndims_coefficients; i++) {
+ std::cout << (size_t)(dims_coefficients[i]);
+ if (i < ndims_coefficients - 1)
+ std::cout << " x ";
+ else
+ std::cout << std::endl;
+ }
+
+ // Reshape coefficients
+ size_t nr_elements_first_dims = 1;
+ for (int i = 0; i < ndims_coefficients - 1; i++)
+ nr_elements_first_dims *= dims_coefficients[i];
+ std::cout << nr_elements_first_dims << " x "
+ << (size_t)(dims_coefficients[ndims_coefficients - 1]) << std::endl;
+
+ m_coefficients =
+ Eigen::ArrayXXcd(nr_elements_first_dims,
+ (size_t)dims_coefficients[ndims_coefficients - 1]);
+ m_coefficients_shape = std::vector<size_t>(
+ dims_coefficients, dims_coefficients + ndims_coefficients);
+ dataset.read(m_coefficients.data(), h5_dcomplex);
+
+ //=====================================
+ dataset = h5file.openDataSet("frequencies");
+ dataspace = dataset.getSpace();
+ nr_elements = dataspace.getSimpleExtentNpoints();
+
+ m_frequencies = std::vector<double>(nr_elements);
+
+ dataset.read(m_frequencies.data(), H5::PredType::NATIVE_DOUBLE);
+ std::cout << "Frequencies:" << std::endl;
+ for (auto freq : m_frequencies) {
+ std::cout << freq << " ";
+ }
+ std::cout << std::endl << std::endl;
+
+ dataset = h5file.openDataSet("nms");
+ dataspace = dataset.getSpace();
+ nr_elements = dataspace.getSimpleExtentNpoints();
+
+ /*
+ * Get the number of dimensions in the dataspace.
+ */
+ int ndims_nms = dataspace.getSimpleExtentNdims();
+ /*
+ * Get the dimension size of each dimension in the dataspace and
+ * display them.
+ */
+ hsize_t dims_nms[ndims_nms];
+ dataspace.getSimpleExtentDims(dims_nms, NULL);
+
+ std::cout << "Dimensions of nms:" << ndims_nms << std::endl << " ";
+
+ for (int i = 0; i < ndims_nms; i++) {
+ std::cout << (size_t)(dims_nms[i]);
+ if (i < ndims_nms - 1) std::cout << " x ";
+ }
+ std::cout << std::endl;
+
+ m_nms = py::array_t<int>({dims_nms[0], dims_nms[1]});
+ dataset.read(m_nms.mutable_data(), H5::PredType::NATIVE_INT);
}
-py::array_t<std::complex<double>> LobesBeamModel::eval(py::EigenDRef<const Eigen::ArrayXd> theta, py::EigenDRef<const Eigen::ArrayXd> phi)
-{
+py::array_t<std::complex<double>> LobesBeamModel::eval(
+ py::EigenDRef<const Eigen::ArrayXd> theta,
+ py::EigenDRef<const Eigen::ArrayXd> phi) {
- double beta;
- for (auto [i, freq] : enumerate(m_frequencies)) {
- beta = 2.0 * M_PI * freq / 2.99792458e8;
- }
+ double beta;
+ for (auto [i, freq] : enumerate(m_frequencies)) {
+ beta = 2.0 * M_PI * freq / 2.99792458e8;
+ }
+ beta = 2.0 * M_PI * m_frequencies[0] / 2.99792458e8;
- beta = 2.0 * M_PI * m_frequencies[0] / 2.99792458e8;
+ int s, m, n;
- int s,m,n;
+ auto nms = m_nms.unchecked<2>();
- auto nms = m_nms.unchecked<2>();
+ // auto result = py::array_t<std::complex<double>>({2*theta.size(),
+ // m_nms.shape(0)}); Eigen::Map<Eigen::MatrixXcd>
+ // Base(result.mutable_data(), 2*theta.size(), m_nms.shape(0));
-// auto result = py::array_t<std::complex<double>>({2*theta.size(), m_nms.shape(0)});
-// Eigen::Map<Eigen::MatrixXcd> Base(result.mutable_data(), 2*theta.size(), m_nms.shape(0));
+ Eigen::MatrixXcd Base(2 * theta.size(), m_nms.shape(0));
- Eigen::MatrixXcd Base(2*theta.size(), m_nms.shape(0));
+ for (int i = 0; i < m_nms.shape(0); i++) {
+ int n = nms(i, 0);
+ int m = nms(i, 1);
+ int s = nms(i, 2);
+ // std::cout << i << ": (" << n << ", " << m << ", " << s << ")" <<
+ // std::endl;
+ Eigen::VectorXcd e_phi, e_theta;
+ // Compute SWF
+ std::tie(e_theta, e_phi) = F4far_new(s, m, n, theta, phi, beta);
+ // Interleave
+ Base.col(i)(Eigen::seq(0, Eigen::last, 2)) = e_theta;
+ Base.col(i)(Eigen::seq(1, Eigen::last, 2)) = e_phi;
+ }
- for(int i = 0; i < m_nms.shape(0); i++) {
- int n = nms(i,0);
- int m = nms(i,1);
- int s = nms(i,2);
-// std::cout << i << ": (" << n << ", " << m << ", " << s << ")" << std::endl;
- Eigen::VectorXcd e_phi, e_theta;
- // Compute SWF
- std::tie(e_theta, e_phi) = F4far_new(s, m, n, theta, phi, beta);
- // Interleave
- Base.col(i)(Eigen::seq(0,Eigen::last,2)) = e_theta;
- Base.col(i)(Eigen::seq(1,Eigen::last,2)) = e_phi;
- }
+ std::cout << Base.rows() << ", " << m_coefficients.rows() << std::endl;
+ std::cout << Base.cols() << ", " << m_coefficients.cols() << std::endl;
- std::cout << Base.rows() << ", " << m_coefficients.rows() << std::endl;
- std::cout << Base.cols() << ", " << m_coefficients.cols() << std::endl;
+ // Create array to return
+ auto result =
+ py::array_t<std::complex<double>>({Base.rows(), m_coefficients.rows()});
- // Create array to return
- auto result = py::array_t<std::complex<double>>({Base.rows(), m_coefficients.rows()});
+ // Map the result array to an Eigen array, to be able to assign it the result
+ // of an Eigen matrix multiplication
+ Eigen::Map<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic,
+ Eigen::RowMajor>>
+ result_matrix(result.mutable_data(), Base.rows(), m_coefficients.rows());
- // Map the result array to an Eigen array, to be able to assign it the result of an Eigen matrix multiplication
- Eigen::Map<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>
- result_matrix(result.mutable_data(), Base.rows(), m_coefficients.rows());
-
- // Multiply base functions by coefficients
- result_matrix = Base * m_coefficients.matrix().transpose();
-
- return result;
+ // Multiply base functions by coefficients
+ result_matrix = Base * m_coefficients.matrix().transpose();
+ return result;
}
-
-// py::array_t<std::complex<double>> LobesBeamModel::eval(py::array_t<double> theta, py::array_t<double> phi)
+// py::array_t<std::complex<double>> LobesBeamModel::eval(py::array_t<double>
+// theta, py::array_t<double> phi)
// {
// py::buffer_info theta_buffer = theta.request();
// py::buffer_info phi_buffer = phi.request();
@@ -307,7 +306,8 @@ py::array_t<std::complex<double>> LobesBeamModel::eval(py::EigenDRef<const Eigen
// Eigen::Map<Eigen::ArrayXd> theta_((double *) theta_buffer.ptr, size);
// Eigen::Map<Eigen::ArrayXd> phi_(size, (double *) phi_buffer.ptr);
//
-// // std::pair<Eigen::VectorXcd,Eigen::VectorXcd> F4far_new(s, m, n, theta, phi, beta)
+// // std::pair<Eigen::VectorXcd,Eigen::VectorXcd> F4far_new(s, m, n, theta,
+// phi, beta)
//
// int s,m,n;
//
@@ -326,26 +326,25 @@ py::array_t<std::complex<double>> LobesBeamModel::eval(py::EigenDRef<const Eigen
//
// }
-
PYBIND11_MODULE(lobes, m) {
- m.doc() = "LOBES module"; // optional module docstring
-
- m.def("plustwo", &plustwo, "A function which adds two to a number");
+ m.doc() = "LOBES module"; // optional module docstring
- m.def("assoc_legendre", &std::assoc_legendre<double>, "Associated legendre function");
+ m.def("plustwo", &plustwo, "A function which adds two to a number");
- m.def("scale", [](py::EigenDRef<Eigen::MatrixXd> m, double c) { m *= c; });
+ m.def("assoc_legendre", &std::assoc_legendre<double>,
+ "Associated legendre function");
- m.def("big_mat", &big_mat, "Big matrix");
+ m.def("scale", [](py::EigenDRef<Eigen::MatrixXd> m, double c) { m *= c; });
- m.def("P", &P, "Legendre");
- m.def("Pacc", &Pacc, "Legendre");
- m.def("F4far_new", &F4far_new, "F4far_new");
- m.def("element_response_lba", &element_response_lba, "element_response_lba");
+ m.def("big_mat", &big_mat, "Big matrix");
- py::class_<LobesBeamModel>(m, "LobesBeamModel")
- .def(py::init<const std::string &>())
- .def("eval", &LobesBeamModel::eval);
+ m.def("P", &P, "Legendre");
+ m.def("Pacc", &Pacc, "Legendre");
+ m.def("F4far_new", &F4far_new, "F4far_new");
+ m.def("element_response_lba", &element_response_lba, "element_response_lba");
+ py::class_<LobesBeamModel>(m, "LobesBeamModel")
+ .def(py::init<const std::string &>())
+ .def("eval", &LobesBeamModel::eval);
}
diff --git a/makeresponseimage.cc b/makeresponseimage.cc
index aa545ed0ca194e826fd40a80e320758ca6855a62..850bc7929edd683bfaf3aaa41e4878dddeb76364 100644
--- a/makeresponseimage.cc
+++ b/makeresponseimage.cc
@@ -61,576 +61,547 @@ using namespace LOFAR;
using namespace LOFAR::StationResponse;
using LOFAR::operator<<;
-namespace
-{
- /*!
- * \brief Map of ITRF directions required to compute an image of the
- * station beam.
- *
- * The station beam library uses the ITRF coordinate system to express
- * station positions and source directions. Since the Earth moves with
- * respect to the sky, the ITRF coordinates of a source vary with time.
- * This structure stores the ITRF coordinates for the station and tile beam
- * former reference directions, as well as for a grid of points on the sky,
- * along with the time for which these ITRF coordinates are valid.
- */
- struct ITRFDirectionMap
- {
- /*!
- * \brief The time for which this ITRF direction map is valid (MJD(UTC)
- * in seconds).
- */
- double_t time0;
-
- /*!
- * \brief Station beam former reference direction expressed in ITRF
- * coordinates.
- */
- vector3r_t station0;
-
- /*!
- * \brief Tile beam former reference direction expressed in ITRF
- * coordinates.
- */
- vector3r_t tile0;
-
- /*!
- * \brief ITRF coordinates for a grid of points on the sky.
- */
- Matrix<vector3r_t> directions;
- };
-
- /*!
- * \brief Create an ITRFDirectionMap.
- *
- * \param coordinates Sky coordinate system definition.
- * \param shape Number of points along the RA and DEC axis.
- * \param epoch Time for which to compute the ITRF coordinates.
- * \param position0 Station beam former reference position (phase
- * reference).
- * \param station0 Station beam former reference direction (pointing).
- * \param tile0 Tile beam former reference direction (pointing).
- */
- ITRFDirectionMap makeDirectionMap(const DirectionCoordinate &coordinates,
- const IPosition &shape, const MEpoch &epoch, const MPosition &position0,
- const MDirection &station0, const MDirection &tile0);
-
- /*!
- * \brief Create a DirectionCoordinate instance that defines an image
- * coordinate system on the sky (J2000).
- *
- * \param reference Direction that corresponds to the origin of the
- * coordinate system.
- * \param size Number of points along each axis (RA, DEC). The index of the
- * origin of the coordinate system is (size / 2, size / 2).
- * \param delta Angular step size in radians (assumed to be the same for
- * both axes).
- */
- DirectionCoordinate makeCoordinates(const MDirection &reference,
- unsigned int size, double delta);
-
- /*!
- * \brief Convert an ITRF position given as a StationResponse::vector3r_t
- * instance to a casacore::MPosition.
- */
- MPosition toMPositionITRF(const vector3r_t &position);
-
- /*!
- * \brief Convert a casacore::MPosition instance to a
- * StationResponse::vector3r_t instance.
- */
- vector3r_t fromMPosition(const MPosition &position);
-
- /*!
- * \brief Convert a casacore::MDirection instance to a
- * StationResponse::vector3r_t instance.
- */
- vector3r_t fromMDirection(const MDirection &direction);
-
- /*!
- * \brief Check if the specified column exists as a column of the specified
- * table.
- *
- * \param table The Table instance to check.
- * \param column The name of the column.
- */
- bool hasColumn(const Table &table, const string &column);
-
- /*!
- * \brief Check if the specified sub-table exists as a sub-table of the
- * specified table.
- *
- * \param table The Table instance to check.
- * \param name The name of the sub-table.
- */
- bool hasSubTable(const Table &table, const string &name);
-
- /*!
- * \brief Provide access to a sub-table by name.
- *
- * \param table The Table instance to which the sub-table is associated.
- * \param name The name of the sub-table.
- */
- Table getSubTable(const Table &table, const string &name);
-
- /*!
- * \brief Attempt to read the position of the observatory. If the
- * observatory position is unknown, the specified default position is
- * returned.
- *
- * \param ms MeasurementSet to read the observatory position from.
- * \param idObservation Identifier that determines of which observation the
- * observatory position should be read.
- * \param defaultPosition The position that will be returned if the
- * observatory position is unknown.
- */
- MPosition readObservatoryPosition(const MeasurementSet &ms,
- unsigned int idObservation, const MPosition &defaultPosition);
-
- /*!
- * \brief Read the list of unique timestamps.
- *
- * \param ms MeasurementSet to read the list of unique timestamps from.
- */
- Vector<Double> readUniqueTimes(const MeasurementSet &ms);
-
- /*!
- * \brief Read the reference frequency of the subband associated to the
- * specified data description identifier.
- *
- * \param ms MeasurementSet to read the reference frequency from.
- * \param idDataDescription Identifier that determines of which subband the
- * reference frequency should be read.
- */
- double readFreqReference(const MeasurementSet &ms,
- unsigned int idDataDescription);
-
- /*!
- * \brief Read the phase reference direction.
- *
- * \param ms MeasurementSet to read the phase reference direction from.
- * \param idField Identifier of the field of which the phase reference
- * direction should be read.
- */
- MDirection readPhaseReference(const MeasurementSet &ms,
- unsigned int idField);
-
- /*!
- * \brief Read the station beam former reference direction.
- *
- * \param ms MeasurementSet to read the station beam former reference
- * direction from.
- * \param idField Identifier of the field of which the station beam former
- * reference direction should be read.
- */
- MDirection readDelayReference(const MeasurementSet &ms,
- unsigned int idField);
-
- /*!
- * \brief Read the station beam former reference direction.
- *
- * \param ms MeasurementSet to read the tile beam former reference
- * direction from.
- * \param idField Identifier of the field of which the tile beam former
- * reference direction should be read.
- */
- MDirection readTileReference(const MeasurementSet &ms,
- unsigned int idField);
-
- /*!
- * \brief Store the specified cube of pixel data as a CASA image.
- *
- * \param data Pixel data. The third dimension is assumed to be of length
- * 4, referring to the correlation products XX, XY, YX, YY (in this order).
- * \param coordinates Sky coordinate system definition.
- * \param frequency Frequency for which the pixel data is valid (Hz).
- * \param name File name of the output image.
- */
- template <class T>
- void store(const Cube<T> &data, const DirectionCoordinate &coordinates,
- double frequency, const string &name);
-
- /*!
- * \brief Convert a string to a CASA Quantity (value with unit).
- */
- Quantity readQuantity(const String &in);
-
- /*!
- * \brief Remove all elements from the range [first, last) that fall
- * outside the interval [min, max].
- *
- * This function returns an iterator new_last such that the range [first,
- * new_last) contains no elements that fall outside the interval [min,
- * max]. The iterators in the range [new_last, last) are all still
- * dereferenceable, but the elements that they point to are unspecified.
- * The order of the elements that are not removed is unchanged.
- */
- template <typename T>
- T filter(T first, T last, int min, int max);
-} //# unnamed namespace
-
-
-int main(int argc, char *argv[])
-{
- TEST_SHOW_VERSION(argc, argv, StationResponse);
- INIT_LOGGER(basename(string(argv[0])));
- Version::show<StationResponseVersion>(cout);
-
- // Parse inputs.
- LOFAR::InputParSet inputs;
- inputs.create("ms", "", "Name of input MeasurementSet", "string");
- inputs.create("stations", "0", "IDs of stations to process", "int vector");
- inputs.create("cellsize", "60arcsec", "Angular pixel size",
- "quantity string");
- inputs.create("size", "256", "Number of pixels along each axis", "int");
- inputs.create("offset", "0s", "Time offset from the start of the MS",
- "quantity string");
- inputs.create("frames", "0", "Number of images that will be generated for"
- " each station (equally spaced over the duration of the MS)", "int");
- inputs.create("abs", "false", "If set to true, store the absolute value of"
- " the beam response instead of the complex value (intended for use with"
- " older versions of casaviewer)", "bool");
- inputs.readArguments(argc, argv);
-
- vector<int> stationIDs(inputs.getIntVector("stations"));
- Quantity cellsize = readQuantity(inputs.getString("cellsize"));
- unsigned int size = max(inputs.getInt("size"), 1);
- Quantity offset = readQuantity(inputs.getString("offset"));
- unsigned int nFrames = max(inputs.getInt("frames"), 1);
- Bool abs = inputs.getBool("abs");
-
- // Open MS.
- MeasurementSet ms(inputs.getString("ms"));
- unsigned int idObservation = 0, idField = 0, idDataDescription = 0;
-
- // Read station meta-data.
- vector<Station::Ptr> stations;
- readStations(ms, std::back_inserter(stations));
-
- // Remove illegal station indices.
- stationIDs.erase(filter(stationIDs.begin(), stationIDs.end(), 0,
- static_cast<int>(stations.size()) - 1), stationIDs.end());
-
- // Read unique timestamps
- Table selection =
- ms(ms.col("OBSERVATION_ID") == static_cast<Int>(idObservation)
- && ms.col("FIELD_ID") == static_cast<Int>(idField)
- && ms.col("DATA_DESC_ID") == static_cast<Int>(idDataDescription));
- Vector<Double> time = readUniqueTimes(selection);
-
- // Read reference frequency.
- double refFrequency = readFreqReference(ms, idDataDescription);
-
- // Use the position of the first selected station as the array reference
- // position if the observatory position cannot be found.
- MPosition refPosition = readObservatoryPosition(ms, idField,
- toMPositionITRF(stations.front()->position()));
-
- // Read phase reference direction.
- MDirection refPhase = readPhaseReference(ms, idField);
-
- // Read delay reference direction.
- MDirection refDelay = readDelayReference(ms, idField);
-
- // Read tile reference direction.
- MDirection refTile = readTileReference(ms, idField);
-
- // Create image coordinate system.
- IPosition shape(2, size, size);
- DirectionCoordinate coordinates = makeCoordinates(refPhase, size,
- cellsize.getValue("rad"));
-
- // Compute station response images.
- Cube<Complex> response(size, size, 4);
-
- MEpoch refEpoch;
- Quantity refTime(time(0), "s");
- refTime = refTime + offset;
- Quantity deltaTime((time(time.size() - 1) - time(0) - offset.getValue("s"))
- / (nFrames - 1), "s");
-
- for(size_t j = 0; j < nFrames; ++j)
- {
- cout << "[ Frame: " << (j + 1) << "/" << nFrames << " Offset: +"
- << refTime.getValue() - time(0) << " s ]" << endl;
-
- // Update reference epoch.
- refEpoch.set(refTime);
-
- cout << "Creating ITRF direction map... " << flush;
- ITRFDirectionMap directionMap = makeDirectionMap(coordinates, shape,
- refEpoch, refPosition, refDelay, refTile);
- cout << "done." << endl;
-
- cout << "Computing response images... " << flush;
- for(vector<int>::const_iterator it = stationIDs.begin(),
- end = stationIDs.end(); it != end; ++it)
- {
- Station::ConstPtr station = stations[*it];
- cout << *it << ":" << station->name() << " " << flush;
-
- for(size_t y = 0; y < size; ++y)
- {
- for(size_t x = 0; x < size; ++x)
- {
- matrix22c_t E = station->response(directionMap.time0,
- refFrequency, directionMap.directions(x, y),
- refFrequency, directionMap.station0,
- directionMap.tile0);
-
- response(x, y, 0) = E[0][0];
- response(x, y, 1) = E[0][1];
- response(x, y, 2) = E[1][0];
- response(x, y, 3) = E[1][1];
- }
- }
-
- std::ostringstream oss;
- oss << "response-" << station->name() << "-frame-" << (j + 1)
- << ".img";
-
- if(abs)
- {
- store(Cube<Float>(amplitude(response)), coordinates,
- refFrequency, oss.str());
- }
- else
- {
- store(response, coordinates, refFrequency, oss.str());
- }
+namespace {
+/*!
+ * \brief Map of ITRF directions required to compute an image of the
+ * station beam.
+ *
+ * The station beam library uses the ITRF coordinate system to express
+ * station positions and source directions. Since the Earth moves with
+ * respect to the sky, the ITRF coordinates of a source vary with time.
+ * This structure stores the ITRF coordinates for the station and tile beam
+ * former reference directions, as well as for a grid of points on the sky,
+ * along with the time for which these ITRF coordinates are valid.
+ */
+struct ITRFDirectionMap {
+ /*!
+ * \brief The time for which this ITRF direction map is valid (MJD(UTC)
+ * in seconds).
+ */
+ double_t time0;
+
+ /*!
+ * \brief Station beam former reference direction expressed in ITRF
+ * coordinates.
+ */
+ vector3r_t station0;
+
+ /*!
+ * \brief Tile beam former reference direction expressed in ITRF
+ * coordinates.
+ */
+ vector3r_t tile0;
+
+ /*!
+ * \brief ITRF coordinates for a grid of points on the sky.
+ */
+ Matrix<vector3r_t> directions;
+};
+
+/*!
+ * \brief Create an ITRFDirectionMap.
+ *
+ * \param coordinates Sky coordinate system definition.
+ * \param shape Number of points along the RA and DEC axis.
+ * \param epoch Time for which to compute the ITRF coordinates.
+ * \param position0 Station beam former reference position (phase
+ * reference).
+ * \param station0 Station beam former reference direction (pointing).
+ * \param tile0 Tile beam former reference direction (pointing).
+ */
+ITRFDirectionMap makeDirectionMap(const DirectionCoordinate &coordinates,
+ const IPosition &shape, const MEpoch &epoch,
+ const MPosition &position0,
+ const MDirection &station0,
+ const MDirection &tile0);
+
+/*!
+ * \brief Create a DirectionCoordinate instance that defines an image
+ * coordinate system on the sky (J2000).
+ *
+ * \param reference Direction that corresponds to the origin of the
+ * coordinate system.
+ * \param size Number of points along each axis (RA, DEC). The index of the
+ * origin of the coordinate system is (size / 2, size / 2).
+ * \param delta Angular step size in radians (assumed to be the same for
+ * both axes).
+ */
+DirectionCoordinate makeCoordinates(const MDirection &reference,
+ unsigned int size, double delta);
+
+/*!
+ * \brief Convert an ITRF position given as a StationResponse::vector3r_t
+ * instance to a casacore::MPosition.
+ */
+MPosition toMPositionITRF(const vector3r_t &position);
+
+/*!
+ * \brief Convert a casacore::MPosition instance to a
+ * StationResponse::vector3r_t instance.
+ */
+vector3r_t fromMPosition(const MPosition &position);
+
+/*!
+ * \brief Convert a casacore::MDirection instance to a
+ * StationResponse::vector3r_t instance.
+ */
+vector3r_t fromMDirection(const MDirection &direction);
+
+/*!
+ * \brief Check if the specified column exists as a column of the specified
+ * table.
+ *
+ * \param table The Table instance to check.
+ * \param column The name of the column.
+ */
+bool hasColumn(const Table &table, const string &column);
+
+/*!
+ * \brief Check if the specified sub-table exists as a sub-table of the
+ * specified table.
+ *
+ * \param table The Table instance to check.
+ * \param name The name of the sub-table.
+ */
+bool hasSubTable(const Table &table, const string &name);
+
+/*!
+ * \brief Provide access to a sub-table by name.
+ *
+ * \param table The Table instance to which the sub-table is associated.
+ * \param name The name of the sub-table.
+ */
+Table getSubTable(const Table &table, const string &name);
+
+/*!
+ * \brief Attempt to read the position of the observatory. If the
+ * observatory position is unknown, the specified default position is
+ * returned.
+ *
+ * \param ms MeasurementSet to read the observatory position from.
+ * \param idObservation Identifier that determines of which observation the
+ * observatory position should be read.
+ * \param defaultPosition The position that will be returned if the
+ * observatory position is unknown.
+ */
+MPosition readObservatoryPosition(const MeasurementSet &ms,
+ unsigned int idObservation,
+ const MPosition &defaultPosition);
+
+/*!
+ * \brief Read the list of unique timestamps.
+ *
+ * \param ms MeasurementSet to read the list of unique timestamps from.
+ */
+Vector<Double> readUniqueTimes(const MeasurementSet &ms);
+
+/*!
+ * \brief Read the reference frequency of the subband associated to the
+ * specified data description identifier.
+ *
+ * \param ms MeasurementSet to read the reference frequency from.
+ * \param idDataDescription Identifier that determines of which subband the
+ * reference frequency should be read.
+ */
+double readFreqReference(const MeasurementSet &ms,
+ unsigned int idDataDescription);
+
+/*!
+ * \brief Read the phase reference direction.
+ *
+ * \param ms MeasurementSet to read the phase reference direction from.
+ * \param idField Identifier of the field of which the phase reference
+ * direction should be read.
+ */
+MDirection readPhaseReference(const MeasurementSet &ms, unsigned int idField);
+
+/*!
+ * \brief Read the station beam former reference direction.
+ *
+ * \param ms MeasurementSet to read the station beam former reference
+ * direction from.
+ * \param idField Identifier of the field of which the station beam former
+ * reference direction should be read.
+ */
+MDirection readDelayReference(const MeasurementSet &ms, unsigned int idField);
+
+/*!
+ * \brief Read the station beam former reference direction.
+ *
+ * \param ms MeasurementSet to read the tile beam former reference
+ * direction from.
+ * \param idField Identifier of the field of which the tile beam former
+ * reference direction should be read.
+ */
+MDirection readTileReference(const MeasurementSet &ms, unsigned int idField);
+
+/*!
+ * \brief Store the specified cube of pixel data as a CASA image.
+ *
+ * \param data Pixel data. The third dimension is assumed to be of length
+ * 4, referring to the correlation products XX, XY, YX, YY (in this order).
+ * \param coordinates Sky coordinate system definition.
+ * \param frequency Frequency for which the pixel data is valid (Hz).
+ * \param name File name of the output image.
+ */
+template <class T>
+void store(const Cube<T> &data, const DirectionCoordinate &coordinates,
+ double frequency, const string &name);
+
+/*!
+ * \brief Convert a string to a CASA Quantity (value with unit).
+ */
+Quantity readQuantity(const String &in);
+
+/*!
+ * \brief Remove all elements from the range [first, last) that fall
+ * outside the interval [min, max].
+ *
+ * This function returns an iterator new_last such that the range [first,
+ * new_last) contains no elements that fall outside the interval [min,
+ * max]. The iterators in the range [new_last, last) are all still
+ * dereferenceable, but the elements that they point to are unspecified.
+ * The order of the elements that are not removed is unchanged.
+ */
+template <typename T>
+T filter(T first, T last, int min, int max);
+} // namespace
+
+int main(int argc, char *argv[]) {
+ TEST_SHOW_VERSION(argc, argv, StationResponse);
+ INIT_LOGGER(basename(string(argv[0])));
+ Version::show<StationResponseVersion>(cout);
+
+ // Parse inputs.
+ LOFAR::InputParSet inputs;
+ inputs.create("ms", "", "Name of input MeasurementSet", "string");
+ inputs.create("stations", "0", "IDs of stations to process", "int vector");
+ inputs.create("cellsize", "60arcsec", "Angular pixel size",
+ "quantity string");
+ inputs.create("size", "256", "Number of pixels along each axis", "int");
+ inputs.create("offset", "0s", "Time offset from the start of the MS",
+ "quantity string");
+ inputs.create("frames", "0",
+ "Number of images that will be generated for"
+ " each station (equally spaced over the duration of the MS)",
+ "int");
+ inputs.create(
+ "abs", "false",
+ "If set to true, store the absolute value of"
+ " the beam response instead of the complex value (intended for use with"
+ " older versions of casaviewer)",
+ "bool");
+ inputs.readArguments(argc, argv);
+
+ vector<int> stationIDs(inputs.getIntVector("stations"));
+ Quantity cellsize = readQuantity(inputs.getString("cellsize"));
+ unsigned int size = max(inputs.getInt("size"), 1);
+ Quantity offset = readQuantity(inputs.getString("offset"));
+ unsigned int nFrames = max(inputs.getInt("frames"), 1);
+ Bool abs = inputs.getBool("abs");
+
+ // Open MS.
+ MeasurementSet ms(inputs.getString("ms"));
+ unsigned int idObservation = 0, idField = 0, idDataDescription = 0;
+
+ // Read station meta-data.
+ vector<Station::Ptr> stations;
+ readStations(ms, std::back_inserter(stations));
+
+ // Remove illegal station indices.
+ stationIDs.erase(filter(stationIDs.begin(), stationIDs.end(), 0,
+ static_cast<int>(stations.size()) - 1),
+ stationIDs.end());
+
+ // Read unique timestamps
+ Table selection =
+ ms(ms.col("OBSERVATION_ID") == static_cast<Int>(idObservation) &&
+ ms.col("FIELD_ID") == static_cast<Int>(idField) &&
+ ms.col("DATA_DESC_ID") == static_cast<Int>(idDataDescription));
+ Vector<Double> time = readUniqueTimes(selection);
+
+ // Read reference frequency.
+ double refFrequency = readFreqReference(ms, idDataDescription);
+
+ // Use the position of the first selected station as the array reference
+ // position if the observatory position cannot be found.
+ MPosition refPosition = readObservatoryPosition(
+ ms, idField, toMPositionITRF(stations.front()->position()));
+
+ // Read phase reference direction.
+ MDirection refPhase = readPhaseReference(ms, idField);
+
+ // Read delay reference direction.
+ MDirection refDelay = readDelayReference(ms, idField);
+
+ // Read tile reference direction.
+ MDirection refTile = readTileReference(ms, idField);
+
+ // Create image coordinate system.
+ IPosition shape(2, size, size);
+ DirectionCoordinate coordinates =
+ makeCoordinates(refPhase, size, cellsize.getValue("rad"));
+
+ // Compute station response images.
+ Cube<Complex> response(size, size, 4);
+
+ MEpoch refEpoch;
+ Quantity refTime(time(0), "s");
+ refTime = refTime + offset;
+ Quantity deltaTime(
+ (time(time.size() - 1) - time(0) - offset.getValue("s")) / (nFrames - 1),
+ "s");
+
+ for (size_t j = 0; j < nFrames; ++j) {
+ cout << "[ Frame: " << (j + 1) << "/" << nFrames << " Offset: +"
+ << refTime.getValue() - time(0) << " s ]" << endl;
+
+ // Update reference epoch.
+ refEpoch.set(refTime);
+
+ cout << "Creating ITRF direction map... " << flush;
+ ITRFDirectionMap directionMap = makeDirectionMap(
+ coordinates, shape, refEpoch, refPosition, refDelay, refTile);
+ cout << "done." << endl;
+
+ cout << "Computing response images... " << flush;
+ for (vector<int>::const_iterator it = stationIDs.begin(),
+ end = stationIDs.end();
+ it != end; ++it) {
+ Station::ConstPtr station = stations[*it];
+ cout << *it << ":" << station->name() << " " << flush;
+
+ for (size_t y = 0; y < size; ++y) {
+ for (size_t x = 0; x < size; ++x) {
+ matrix22c_t E = station->response(
+ directionMap.time0, refFrequency, directionMap.directions(x, y),
+ refFrequency, directionMap.station0, directionMap.tile0);
+
+ response(x, y, 0) = E[0][0];
+ response(x, y, 1) = E[0][1];
+ response(x, y, 2) = E[1][0];
+ response(x, y, 3) = E[1][1];
}
- cout << endl;
+ }
- refTime = refTime + deltaTime;
+ std::ostringstream oss;
+ oss << "response-" << station->name() << "-frame-" << (j + 1) << ".img";
+
+ if (abs) {
+ store(Cube<Float>(amplitude(response)), coordinates, refFrequency,
+ oss.str());
+ } else {
+ store(response, coordinates, refFrequency, oss.str());
+ }
}
+ cout << endl;
- return 0;
-}
+ refTime = refTime + deltaTime;
+ }
-namespace
-{
- ITRFDirectionMap makeDirectionMap(const DirectionCoordinate &coordinates,
- const IPosition &shape, const MEpoch &epoch, const MPosition &position0,
- const MDirection &station0, const MDirection &tile0)
- {
- ITRFDirectionMap map;
-
- // Convert from MEpoch to a time in MJD(UTC) in seconds.
- MEpoch mEpochUTC = MEpoch::Convert(epoch, MEpoch::Ref(MEpoch::UTC))();
- MVEpoch mvEpochUTC = mEpochUTC.getValue();
- Quantity qEpochUTC = mvEpochUTC.getTime();
- map.time0 = qEpochUTC.getValue("s");
-
- // Create conversion engine J2000 => ITRF at the specified epoch.
- MDirection::Convert convertor = MDirection::Convert(MDirection::J2000,
- MDirection::Ref(MDirection::ITRF, MeasFrame(epoch, position0)));
-
- // Compute station and tile beam former reference directions in ITRF at
- // the specified epoch.
- map.station0 = fromMDirection(convertor(station0));
- map.tile0 = fromMDirection(convertor(tile0));
-
- // Pre-allocate space for the grid of ITRF directions.
- map.directions.resize(shape);
-
- // Compute ITRF directions.
- MDirection world;
- Vector<Double> pixel = coordinates.referencePixel();
- for(pixel(1) = 0.0; pixel(1) < shape(1); ++pixel(1))
- {
- for(pixel(0) = 0.0; pixel(0) < shape(0); ++pixel(0))
- {
- // CoordinateSystem::toWorld(): RA range [-pi,pi], DEC range
- // [-pi/2,pi/2].
- if(coordinates.toWorld(world, pixel))
- {
- map.directions(pixel(0), pixel(1)) =
- fromMDirection(convertor(world));
- }
- }
- }
+ return 0;
+}
- return map;
+namespace {
+ITRFDirectionMap makeDirectionMap(const DirectionCoordinate &coordinates,
+ const IPosition &shape, const MEpoch &epoch,
+ const MPosition &position0,
+ const MDirection &station0,
+ const MDirection &tile0) {
+ ITRFDirectionMap map;
+
+ // Convert from MEpoch to a time in MJD(UTC) in seconds.
+ MEpoch mEpochUTC = MEpoch::Convert(epoch, MEpoch::Ref(MEpoch::UTC))();
+ MVEpoch mvEpochUTC = mEpochUTC.getValue();
+ Quantity qEpochUTC = mvEpochUTC.getTime();
+ map.time0 = qEpochUTC.getValue("s");
+
+ // Create conversion engine J2000 => ITRF at the specified epoch.
+ MDirection::Convert convertor = MDirection::Convert(
+ MDirection::J2000,
+ MDirection::Ref(MDirection::ITRF, MeasFrame(epoch, position0)));
+
+ // Compute station and tile beam former reference directions in ITRF at
+ // the specified epoch.
+ map.station0 = fromMDirection(convertor(station0));
+ map.tile0 = fromMDirection(convertor(tile0));
+
+ // Pre-allocate space for the grid of ITRF directions.
+ map.directions.resize(shape);
+
+ // Compute ITRF directions.
+ MDirection world;
+ Vector<Double> pixel = coordinates.referencePixel();
+ for (pixel(1) = 0.0; pixel(1) < shape(1); ++pixel(1)) {
+ for (pixel(0) = 0.0; pixel(0) < shape(0); ++pixel(0)) {
+ // CoordinateSystem::toWorld(): RA range [-pi,pi], DEC range
+ // [-pi/2,pi/2].
+ if (coordinates.toWorld(world, pixel)) {
+ map.directions(pixel(0), pixel(1)) = fromMDirection(convertor(world));
+ }
}
+ }
- DirectionCoordinate makeCoordinates(const MDirection &reference,
- unsigned int size, double delta)
- {
- MDirection referenceJ2000 = MDirection::Convert(reference,
- MDirection::J2000)();
- Quantum<Vector<Double> > referenceAngles = referenceJ2000.getAngle();
- double ra = referenceAngles.getBaseValue()(0);
- double dec = referenceAngles.getBaseValue()(1);
-
- Matrix<Double> xform(2,2);
- xform = 0.0; xform.diagonal() = 1.0;
- return DirectionCoordinate(MDirection::J2000,
- Projection(Projection::SIN), ra, dec, -delta, delta, xform,
- size / 2, size / 2);
- }
+ return map;
+}
- MPosition toMPositionITRF(const vector3r_t &position)
- {
- MVPosition mvITRF(position[0], position[1], position[2]);
- return MPosition(mvITRF, MPosition::ITRF);
- }
+DirectionCoordinate makeCoordinates(const MDirection &reference,
+ unsigned int size, double delta) {
+ MDirection referenceJ2000 =
+ MDirection::Convert(reference, MDirection::J2000)();
+ Quantum<Vector<Double> > referenceAngles = referenceJ2000.getAngle();
+ double ra = referenceAngles.getBaseValue()(0);
+ double dec = referenceAngles.getBaseValue()(1);
+
+ Matrix<Double> xform(2, 2);
+ xform = 0.0;
+ xform.diagonal() = 1.0;
+ return DirectionCoordinate(MDirection::J2000, Projection(Projection::SIN), ra,
+ dec, -delta, delta, xform, size / 2, size / 2);
+}
- vector3r_t fromMPosition(const MPosition &position)
- {
- MVPosition mvPosition = position.getValue();
- vector3r_t result = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
- return result;
- }
+MPosition toMPositionITRF(const vector3r_t &position) {
+ MVPosition mvITRF(position[0], position[1], position[2]);
+ return MPosition(mvITRF, MPosition::ITRF);
+}
- vector3r_t fromMDirection(const MDirection &direction)
- {
- MVDirection mvDirection = direction.getValue();
- vector3r_t result = {{mvDirection(0), mvDirection(1), mvDirection(2)}};
- return result;
- }
+vector3r_t fromMPosition(const MPosition &position) {
+ MVPosition mvPosition = position.getValue();
+ vector3r_t result = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
+ return result;
+}
- bool hasColumn(const Table &table, const string &column)
- {
- return table.tableDesc().isColumn(column);
- }
+vector3r_t fromMDirection(const MDirection &direction) {
+ MVDirection mvDirection = direction.getValue();
+ vector3r_t result = {{mvDirection(0), mvDirection(1), mvDirection(2)}};
+ return result;
+}
- bool hasSubTable(const Table &table, const string &name)
- {
- return table.keywordSet().isDefined(name);
- }
+bool hasColumn(const Table &table, const string &column) {
+ return table.tableDesc().isColumn(column);
+}
- Table getSubTable(const Table &table, const string &name)
- {
- return table.keywordSet().asTable(name);
- }
+bool hasSubTable(const Table &table, const string &name) {
+ return table.keywordSet().isDefined(name);
+}
- MPosition readObservatoryPosition(const MeasurementSet &ms,
- unsigned int idObservation, const MPosition &defaultPosition)
- {
- // Get the instrument position in ITRF coordinates, or use the centroid
- // of the station positions if the instrument position is unknown.
- ROMSObservationColumns observation(ms.observation());
- ASSERT(observation.nrow() > idObservation);
- ASSERT(!observation.flagRow()(idObservation));
-
- // Read observatory name and try to look-up its position.
- const string observatory = observation.telescopeName()(idObservation);
-
- // Look-up observatory position, default to specified default position.
- MPosition position(defaultPosition);
- MeasTable::Observatory(position, observatory);
- return position;
- }
+Table getSubTable(const Table &table, const string &name) {
+ return table.keywordSet().asTable(name);
+}
- Vector<Double> readUniqueTimes(const MeasurementSet &ms)
- {
- Table tab_sorted = ms.sort("TIME", Sort::Ascending, Sort::HeapSort
- | Sort::NoDuplicates);
+MPosition readObservatoryPosition(const MeasurementSet &ms,
+ unsigned int idObservation,
+ const MPosition &defaultPosition) {
+ // Get the instrument position in ITRF coordinates, or use the centroid
+ // of the station positions if the instrument position is unknown.
+ ROMSObservationColumns observation(ms.observation());
+ ASSERT(observation.nrow() > idObservation);
+ ASSERT(!observation.flagRow()(idObservation));
+
+ // Read observatory name and try to look-up its position.
+ const string observatory = observation.telescopeName()(idObservation);
+
+ // Look-up observatory position, default to specified default position.
+ MPosition position(defaultPosition);
+ MeasTable::Observatory(position, observatory);
+ return position;
+}
- ROScalarColumn<Double> c_time(tab_sorted, "TIME");
- return c_time.getColumn();
- }
+Vector<Double> readUniqueTimes(const MeasurementSet &ms) {
+ Table tab_sorted =
+ ms.sort("TIME", Sort::Ascending, Sort::HeapSort | Sort::NoDuplicates);
- double readFreqReference(const MeasurementSet &ms,
- unsigned int idDataDescription)
- {
- ROMSDataDescColumns desc(ms.dataDescription());
- ASSERT(desc.nrow() > idDataDescription);
- ASSERT(!desc.flagRow()(idDataDescription));
- uInt idWindow = desc.spectralWindowId()(idDataDescription);
+ ROScalarColumn<Double> c_time(tab_sorted, "TIME");
+ return c_time.getColumn();
+}
- ROMSSpWindowColumns window(ms.spectralWindow());
- ASSERT(window.nrow() > idWindow);
- ASSERT(!window.flagRow()(idWindow));
+double readFreqReference(const MeasurementSet &ms,
+ unsigned int idDataDescription) {
+ ROMSDataDescColumns desc(ms.dataDescription());
+ ASSERT(desc.nrow() > idDataDescription);
+ ASSERT(!desc.flagRow()(idDataDescription));
+ uInt idWindow = desc.spectralWindowId()(idDataDescription);
- return window.refFrequency()(idWindow);
- }
+ ROMSSpWindowColumns window(ms.spectralWindow());
+ ASSERT(window.nrow() > idWindow);
+ ASSERT(!window.flagRow()(idWindow));
- MDirection readPhaseReference(const MeasurementSet &ms,
- unsigned int idField)
- {
- ROMSFieldColumns field(ms.field());
- ASSERT(field.nrow() > idField);
- ASSERT(!field.flagRow()(idField));
+ return window.refFrequency()(idWindow);
+}
- return field.phaseDirMeas(idField);
- }
+MDirection readPhaseReference(const MeasurementSet &ms, unsigned int idField) {
+ ROMSFieldColumns field(ms.field());
+ ASSERT(field.nrow() > idField);
+ ASSERT(!field.flagRow()(idField));
- MDirection readDelayReference(const MeasurementSet &ms,
- unsigned int idField)
- {
- ROMSFieldColumns field(ms.field());
- ASSERT(field.nrow() > idField);
- ASSERT(!field.flagRow()(idField));
+ return field.phaseDirMeas(idField);
+}
- return field.delayDirMeas(idField);
- }
+MDirection readDelayReference(const MeasurementSet &ms, unsigned int idField) {
+ ROMSFieldColumns field(ms.field());
+ ASSERT(field.nrow() > idField);
+ ASSERT(!field.flagRow()(idField));
- MDirection readTileReference(const MeasurementSet &ms, unsigned int idField)
- {
- // The MeasurementSet class does not support LOFAR specific columns, so
- // we use ROArrayMeasColumn to read the tile beam reference direction.
- Table tab_field = getSubTable(ms, "FIELD");
-
- static const String columnName = "LOFAR_TILE_BEAM_DIR";
- if(hasColumn(tab_field, columnName))
- {
- ROArrayMeasColumn<MDirection> c_direction(tab_field, columnName);
- if(c_direction.isDefined(idField))
- {
- return c_direction(idField)(IPosition(1, 0));
- }
- }
+ return field.delayDirMeas(idField);
+}
- // By default, the tile beam reference direction is assumed to be equal
- // to the station beam reference direction (for backward compatibility,
- // and for non-HBA measurements).
- return readDelayReference(ms, idField);
- }
+MDirection readTileReference(const MeasurementSet &ms, unsigned int idField) {
+ // The MeasurementSet class does not support LOFAR specific columns, so
+ // we use ROArrayMeasColumn to read the tile beam reference direction.
+ Table tab_field = getSubTable(ms, "FIELD");
- template <class T>
- void store(const Cube<T> &data, const DirectionCoordinate &coordinates,
- double frequency, const string &name)
- {
- ASSERT(data.shape()(2) == 4);
-
- Vector<Int> stokes(4);
- stokes(0) = Stokes::XX;
- stokes(1) = Stokes::XY;
- stokes(2) = Stokes::YX;
- stokes(3) = Stokes::YY;
-
- CoordinateSystem csys;
- csys.addCoordinate(coordinates);
- csys.addCoordinate(StokesCoordinate(stokes));
- csys.addCoordinate(SpectralCoordinate(MFrequency::TOPO, frequency, 0.0,
- 0.0));
-
- PagedImage<T> im(TiledShape(IPosition(4, data.shape()(0),
- data.shape()(1), 4, 1)), csys, name);
- im.putSlice(data, IPosition(4, 0, 0, 0, 0));
+ static const String columnName = "LOFAR_TILE_BEAM_DIR";
+ if (hasColumn(tab_field, columnName)) {
+ ROArrayMeasColumn<MDirection> c_direction(tab_field, columnName);
+ if (c_direction.isDefined(idField)) {
+ return c_direction(idField)(IPosition(1, 0));
}
+ }
- Quantity readQuantity(const String &in)
- {
- Quantity result;
- bool status = Quantity::read(result, in);
- ASSERT(status);
- return result;
- }
+ // By default, the tile beam reference direction is assumed to be equal
+ // to the station beam reference direction (for backward compatibility,
+ // and for non-HBA measurements).
+ return readDelayReference(ms, idField);
+}
- template <typename T>
- T filter(T first, T last, int min, int max)
- {
- T new_last = first;
- for(; first != last; ++first)
- {
- if(*first >= min && *first <= max)
- {
- *new_last++ = *first;
- }
- }
+template <class T>
+void store(const Cube<T> &data, const DirectionCoordinate &coordinates,
+ double frequency, const string &name) {
+ ASSERT(data.shape()(2) == 4);
+
+ Vector<Int> stokes(4);
+ stokes(0) = Stokes::XX;
+ stokes(1) = Stokes::XY;
+ stokes(2) = Stokes::YX;
+ stokes(3) = Stokes::YY;
+
+ CoordinateSystem csys;
+ csys.addCoordinate(coordinates);
+ csys.addCoordinate(StokesCoordinate(stokes));
+ csys.addCoordinate(SpectralCoordinate(MFrequency::TOPO, frequency, 0.0, 0.0));
+
+ PagedImage<T> im(
+ TiledShape(IPosition(4, data.shape()(0), data.shape()(1), 4, 1)), csys,
+ name);
+ im.putSlice(data, IPosition(4, 0, 0, 0, 0));
+}
- return new_last;
+Quantity readQuantity(const String &in) {
+ Quantity result;
+ bool status = Quantity::read(result, in);
+ ASSERT(status);
+ return result;
+}
+
+template <typename T>
+T filter(T first, T last, int min, int max) {
+ T new_last = first;
+ for (; first != last; ++first) {
+ if (*first >= min && *first <= max) {
+ *new_last++ = *first;
}
-} // unnamed namespace.
+ }
+
+ return new_last;
+}
+} // unnamed namespace.
diff --git a/oskar/OSKARDatafile.cc b/oskar/OSKARDatafile.cc
index e2aef0d4ff25f99775efeb47e72436009a061ba5..419103bff658b8933d39b68085e64bd58f8d2810 100644
--- a/oskar/OSKARDatafile.cc
+++ b/oskar/OSKARDatafile.cc
@@ -2,33 +2,29 @@
#include <iostream>
-Datafile::Datafile(
- const std::string& filename)
-{
- // Open file
- std::cout << "read oskar datafile: " << filename << std::endl;
- m_h5_file.reset(new H5::H5File(filename, H5F_ACC_RDONLY));
+Datafile::Datafile(const std::string& filename) {
+ // Open file
+ std::cout << "read oskar datafile: " << filename << std::endl;
+ m_h5_file.reset(new H5::H5File(filename, H5F_ACC_RDONLY));
- // Disable HDF5 error prints
- H5::Exception::dontPrint();
+ // Disable HDF5 error prints
+ H5::Exception::dontPrint();
};
-std::shared_ptr<Dataset> Datafile::get(
- const unsigned int freq)
-{
- std::lock_guard<std::mutex> lock(m_mutex);
+std::shared_ptr<Dataset> Datafile::get(const unsigned int freq) {
+ std::lock_guard<std::mutex> lock(m_mutex);
- // Find dataset for frequency
- auto entry = m_map.find(freq);
+ // Find dataset for frequency
+ auto entry = m_map.find(freq);
- // If found, retrieve pointer to dataset
- if (entry != m_map.end()) {
- return entry->second;
- }
+ // If found, retrieve pointer to dataset
+ if (entry != m_map.end()) {
+ return entry->second;
+ }
- // Read and return dataset
- std::shared_ptr<Dataset> dataset_ptr;
- dataset_ptr.reset(new Dataset(*m_h5_file, freq));
- m_map.insert({freq, dataset_ptr});
- return dataset_ptr;
+ // Read and return dataset
+ std::shared_ptr<Dataset> dataset_ptr;
+ dataset_ptr.reset(new Dataset(*m_h5_file, freq));
+ m_map.insert({freq, dataset_ptr});
+ return dataset_ptr;
}
diff --git a/oskar/OSKARDataset.cc b/oskar/OSKARDataset.cc
index abc73fd48a300da5d6df797339c5b82d6cff7aa2..fdc3e8d7ba56bb806a23d73b8edf26aa22dca41c 100644
--- a/oskar/OSKARDataset.cc
+++ b/oskar/OSKARDataset.cc
@@ -3,77 +3,71 @@
#include "OSKARDataset.h"
-Dataset::Dataset(
- H5::H5File& h5_file,
- const unsigned int freq)
-{
- // Try to read coefficients for this frequency
- std::string dataset_name = std::to_string((int) (freq / 1e6));
+Dataset::Dataset(H5::H5File& h5_file, const unsigned int freq) {
+ // Try to read coefficients for this frequency
+ std::string dataset_name = std::to_string((int)(freq / 1e6));
- try {
- // Open dataset
- H5::DataSet dataset = h5_file.openDataSet(dataset_name);
+ try {
+ // Open dataset
+ H5::DataSet dataset = h5_file.openDataSet(dataset_name);
- // Read dataset dimensions
- H5::DataSpace dataspace = dataset.getSpace();
- unsigned int rank = dataspace.getSimpleExtentNdims();
- assert(rank == m_dataset_rank);
+ // Read dataset dimensions
+ H5::DataSpace dataspace = dataset.getSpace();
+ unsigned int rank = dataspace.getSimpleExtentNdims();
+ assert(rank == m_dataset_rank);
- // Get dimensions
- hsize_t dims[rank];
- dataspace.getSimpleExtentDims(dims, NULL);
- auto nr_elements = dims[0];
- auto nr_coeffs = dims[1];
- assert(dims[2] == 4); // tetm*pol
+ // Get dimensions
+ hsize_t dims[rank];
+ dataspace.getSimpleExtentDims(dims, NULL);
+ auto nr_elements = dims[0];
+ auto nr_coeffs = dims[1];
+ assert(dims[2] == 4); // tetm*pol
- // Coefficient data stored as:
- // [nr_elements][nr_coefficients][4],
- // with inner dimension:
- // (x_te_re, x_te_im), (x_tm_re, x_tm_im),
- // (y_te_re, y_te_im), (y_tm_re, y_tm_im)
- #if defined(DEBUG)
- std::cout << "nr_elements: " << nr_elements << std::endl;
- std::cout << "nr_coeffs: " << nr_coeffs << std::endl;
- #endif
+// Coefficient data stored as:
+// [nr_elements][nr_coefficients][4],
+// with inner dimension:
+// (x_te_re, x_te_im), (x_tm_re, x_tm_im),
+// (y_te_re, y_te_im), (y_tm_re, y_tm_im)
+#if defined(DEBUG)
+ std::cout << "nr_elements: " << nr_elements << std::endl;
+ std::cout << "nr_coeffs: " << nr_coeffs << std::endl;
+#endif
- // Check total number of coefficients to find l_max
- auto l_max_d = sqrt(nr_coeffs + 1) - 1;
- auto l_max = (int)round(l_max_d);
- #if defined(DEBUG)
- std::cout << "l_max: " << l_max << std::endl;
- #endif
+ // Check total number of coefficients to find l_max
+ auto l_max_d = sqrt(nr_coeffs + 1) - 1;
+ auto l_max = (int)round(l_max_d);
+#if defined(DEBUG)
+ std::cout << "l_max: " << l_max << std::endl;
+#endif
- // Sanity check
- assert(l_max * (l_max+2) == nr_coeffs);
+ // Sanity check
+ assert(l_max * (l_max + 2) == nr_coeffs);
- // Set members
- m_nr_elements = nr_elements;
- m_nr_coeffs = nr_coeffs;
- m_l_max = l_max;
+ // Set members
+ m_nr_elements = nr_elements;
+ m_nr_coeffs = nr_coeffs;
+ m_l_max = l_max;
- // Read coefficients into data vector
- m_data.resize(nr_elements * nr_coeffs * 4);
- assert(dims[0]*dims[1]*dims[2]==m_data.size());
- H5::DataType data_type = dataset.getDataType();
- assert(data_type.getSize() == sizeof(std::complex<double>));
- dataset.read(m_data.data(), data_type, dataspace);
- } catch (H5::FileIException& e) {
- std::stringstream message;
- message << "Could not load dataset for frequency " << dataset_name << " Mhz";
- throw std::runtime_error(message.str());
- }
+ // Read coefficients into data vector
+ m_data.resize(nr_elements * nr_coeffs * 4);
+ assert(dims[0] * dims[1] * dims[2] == m_data.size());
+ H5::DataType data_type = dataset.getDataType();
+ assert(data_type.getSize() == sizeof(std::complex<double>));
+ dataset.read(m_data.data(), data_type, dataspace);
+ } catch (H5::FileIException& e) {
+ std::stringstream message;
+ message << "Could not load dataset for frequency " << dataset_name
+ << " Mhz";
+ throw std::runtime_error(message.str());
+ }
}
-size_t Dataset::get_index(
- const unsigned int element) const
-{
- return element * m_nr_coeffs * 4;
+size_t Dataset::get_index(const unsigned int element) const {
+ return element * m_nr_coeffs * 4;
}
-std::complex<double>* Dataset::get_alpha_ptr(
- const unsigned int element)
-{
- assert(element < get_nr_elements());
- size_t index = get_index(element);
- return m_data.data() + index;
+std::complex<double>* Dataset::get_alpha_ptr(const unsigned int element) {
+ assert(element < get_nr_elements());
+ size_t index = get_index(element);
+ return m_data.data() + index;
}
diff --git a/oskar/OSKARElementResponse.cc b/oskar/OSKARElementResponse.cc
index 655a8e6324f9bead791e9b4852e4722d749c261a..7cd0c84cd6a0e7c512d467b9138435ac05e8d6e2 100644
--- a/oskar/OSKARElementResponse.cc
+++ b/oskar/OSKARElementResponse.cc
@@ -4,69 +4,61 @@
#include <iostream>
namespace LOFAR {
-namespace StationResponse{
+namespace StationResponse {
void OSKARElementResponseDipole::response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const
-{
- double dipole_length_m = 1; // TODO
- std::complex<double>* response_ptr = (std::complex<double> *) response;
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const {
+ double dipole_length_m = 1; // TODO
+ std::complex<double>* response_ptr = (std::complex<double>*)response;
- double phi_x = phi;
- double phi_y = phi + M_PI_2;
- oskar_evaluate_dipole_pattern_double(1, &theta, &phi_x, freq, dipole_length_m, response_ptr);
- oskar_evaluate_dipole_pattern_double(1, &theta, &phi_y, freq, dipole_length_m, response_ptr + 2);
+ double phi_x = phi;
+ double phi_y = phi + M_PI_2;
+ oskar_evaluate_dipole_pattern_double(1, &theta, &phi_x, freq, dipole_length_m,
+ response_ptr);
+ oskar_evaluate_dipole_pattern_double(1, &theta, &phi_y, freq, dipole_length_m,
+ response_ptr + 2);
}
-OSKARElementResponseSphericalWave::OSKARElementResponseSphericalWave()
-{
- std::string path = get_path("oskar.h5");
- m_datafile.reset(new Datafile(path));
+OSKARElementResponseSphericalWave::OSKARElementResponseSphericalWave() {
+ std::string path = get_path("oskar.h5");
+ m_datafile.reset(new Datafile(path));
}
void OSKARElementResponseSphericalWave::response(
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const
-{
- // This ElementResponse model is element specific, so an element_id is required
- // to know for what element the response needs to be evaluated
- // A std::invalid_argument exception is thrown although strictly speaking
- // it are not the given arguments that are invalid, but the response(...) method with
- // a different signature should have been called.
- throw std::invalid_argument("OSKARElementResponseSphericalWave: missing argument element_id");
+ double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const {
+ // This ElementResponse model is element specific, so an element_id is
+ // required to know for what element the response needs to be evaluated A
+ // std::invalid_argument exception is thrown although strictly speaking it are
+ // not the given arguments that are invalid, but the response(...) method with
+ // a different signature should have been called.
+ throw std::invalid_argument(
+ "OSKARElementResponseSphericalWave: missing argument element_id");
}
void OSKARElementResponseSphericalWave::response(
- int element_id,
- double freq,
- double theta,
- double phi,
- std::complex<double> (&response)[2][2]) const
-{
- auto dataset = m_datafile->get(freq);
- auto l_max = dataset->get_l_max();
+ int element_id, double freq, double theta, double phi,
+ std::complex<double> (&response)[2][2]) const {
+ auto dataset = m_datafile->get(freq);
+ auto l_max = dataset->get_l_max();
- std::complex<double>* response_ptr = (std::complex<double> *) response;
- std::complex<double>* alpha_ptr = dataset->get_alpha_ptr(element_id);
+ std::complex<double>* response_ptr = (std::complex<double>*)response;
+ std::complex<double>* alpha_ptr = dataset->get_alpha_ptr(element_id);
- double phi_x = phi;
- double phi_y = phi + M_PI_2;
- oskar_evaluate_spherical_wave_sum_double(1, &theta, &phi_x, &phi_y, l_max, alpha_ptr, response_ptr);
+ double phi_x = phi;
+ double phi_y = phi + M_PI_2;
+ oskar_evaluate_spherical_wave_sum_double(1, &theta, &phi_x, &phi_y, l_max,
+ alpha_ptr, response_ptr);
}
std::string OSKARElementResponseSphericalWave::get_path(
- const char* filename) const
-{
- std::stringstream ss;
- ss << LOFARBEAM_DATA_DIR << "/";
- ss << filename;
- return ss.str();
+ const char* filename) const {
+ std::stringstream ss;
+ ss << LOFARBEAM_DATA_DIR << "/";
+ ss << filename;
+ return ss.str();
}
-} // namespace StationResponse
-} // namespace LOFAR
+} // namespace StationResponse
+} // namespace LOFAR
diff --git a/oskar/oskar_evaluate_dipole_pattern.cc b/oskar/oskar_evaluate_dipole_pattern.cc
index 69ecf8776b11b19c7a37d927b0b371e32c4f825b..e09ebcc9d7b0042d970b58b6a0741b0d022d808c 100644
--- a/oskar/oskar_evaluate_dipole_pattern.cc
+++ b/oskar/oskar_evaluate_dipole_pattern.cc
@@ -31,86 +31,72 @@
#include "oskar.h"
-template<typename FP, typename FP2>
-inline void oskar_dipole(
- FP theta,
- FP phi,
- FP kL,
- FP2& e_theta,
- FP2& e_phi)
-{
- FP sin_phi, cos_phi;
- FP sin_theta, cos_theta;
- oskar_sincos(phi, &sin_phi, &cos_phi);
- oskar_sincos(theta, &sin_theta, &cos_theta);
- FP cos_kL = (FP) cos(kL);
- const FP denom = (FP)1 + cos_phi*cos_phi * (cos_theta*cos_theta - (FP)1);
- if (denom == (FP)0)
- e_theta.x = e_theta.y = e_phi.x = e_phi.y = (FP)0;
- else {
- const FP q = kL * cos_phi * sin_theta;
- const FP cos_q = cos(q);
- const FP t = (cos_q - cos_kL) / denom;
- e_theta.x = -cos_phi * cos_theta * t;
- e_phi.x = sin_phi * t;
- e_phi.y = e_theta.y = (FP)0;
- }
+template <typename FP, typename FP2>
+inline void oskar_dipole(FP theta, FP phi, FP kL, FP2& e_theta, FP2& e_phi) {
+ FP sin_phi, cos_phi;
+ FP sin_theta, cos_theta;
+ oskar_sincos(phi, &sin_phi, &cos_phi);
+ oskar_sincos(theta, &sin_theta, &cos_theta);
+ FP cos_kL = (FP)cos(kL);
+ const FP denom = (FP)1 + cos_phi * cos_phi * (cos_theta * cos_theta - (FP)1);
+ if (denom == (FP)0)
+ e_theta.x = e_theta.y = e_phi.x = e_phi.y = (FP)0;
+ else {
+ const FP q = kL * cos_phi * sin_theta;
+ const FP cos_q = cos(q);
+ const FP t = (cos_q - cos_kL) / denom;
+ e_theta.x = -cos_phi * cos_theta * t;
+ e_phi.x = sin_phi * t;
+ e_phi.y = e_theta.y = (FP)0;
+ }
}
-template<typename FP, typename FP2>
-void oskar_evaluate_dipole_pattern(
- const int num_points,
- const FP* theta,
- const FP* phi,
- const FP kL,
- const int stride,
- const int E_theta_offset,
- const int E_phi_offset,
- FP2* e_theta,
- FP2* e_phi)
-{
- for (int i = 0; i < num_points; i++) {
- const int i_out = i * stride;\
- const int theta_out = i_out + E_theta_offset;\
- const int phi_out = i_out + E_phi_offset;\
- oskar_dipole<FP, FP2>(theta[i], phi[i], kL, e_theta[theta_out], e_phi[phi_out]);
- }
+template <typename FP, typename FP2>
+void oskar_evaluate_dipole_pattern(const int num_points, const FP* theta,
+ const FP* phi, const FP kL, const int stride,
+ const int E_theta_offset,
+ const int E_phi_offset, FP2* e_theta,
+ FP2* e_phi) {
+ for (int i = 0; i < num_points; i++) {
+ const int i_out = i * stride;
+ const int theta_out = i_out + E_theta_offset;
+ const int phi_out = i_out + E_phi_offset;
+ oskar_dipole<FP, FP2>(theta[i], phi[i], kL, e_theta[theta_out],
+ e_phi[phi_out]);
+ }
}
-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)
-{
- const int stride = 4;
- const int offset = 0;
- const int E_theta_offset = offset;
- const int E_phi_offset = offset + 1;
- const double kL = dipole_length_m * (M_PI * freq_hz / 299792458);
- double2* pattern_ptr = (double2 *) 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) {
+ const int stride = 4;
+ const int offset = 0;
+ const int E_theta_offset = offset;
+ const int E_phi_offset = offset + 1;
+ const double kL = dipole_length_m * (M_PI * freq_hz / 299792458);
+ double2* pattern_ptr = (double2*)pattern;
- oskar_evaluate_dipole_pattern<double, double2>(
- num_points, theta, phi, kL, stride, E_theta_offset, E_phi_offset, pattern_ptr, pattern_ptr);
+ oskar_evaluate_dipole_pattern<double, double2>(
+ num_points, theta, phi, kL, stride, E_theta_offset, E_phi_offset,
+ pattern_ptr, pattern_ptr);
}
-void oskar_evaluate_dipole_pattern_float(
- const int num_points,
- const float* theta,
- const float* phi,
- const float freq_hz,
- const float dipole_length_m,
- std::complex<double>* pattern)
-{
- const int stride = 4;
- const int offset = 0;
- const int E_theta_offset = offset;
- const int E_phi_offset = offset + 1;
- const float kL = dipole_length_m * (M_PI * freq_hz / 299792458);
- float2* pattern_ptr = (float2 *) pattern;
+void oskar_evaluate_dipole_pattern_float(const int num_points,
+ const float* theta, const float* phi,
+ const float freq_hz,
+ const float dipole_length_m,
+ std::complex<double>* pattern) {
+ const int stride = 4;
+ const int offset = 0;
+ const int E_theta_offset = offset;
+ const int E_phi_offset = offset + 1;
+ const float kL = dipole_length_m * (M_PI * freq_hz / 299792458);
+ float2* pattern_ptr = (float2*)pattern;
- oskar_evaluate_dipole_pattern<float, float2>(
- num_points, theta, phi, kL, stride, E_theta_offset, E_phi_offset, pattern_ptr, pattern_ptr);
+ oskar_evaluate_dipole_pattern<float, float2>(
+ num_points, theta, phi, kL, stride, E_theta_offset, E_phi_offset,
+ pattern_ptr, pattern_ptr);
}
diff --git a/oskar/oskar_evaluate_spherical_wave_sum.cc b/oskar/oskar_evaluate_spherical_wave_sum.cc
index 8d543a2888e056ac8189d2b774e2d55821118d0c..26f40781c604a3ef74dfbebb7dbebf3e43ffa73f 100644
--- a/oskar/oskar_evaluate_spherical_wave_sum.cc
+++ b/oskar/oskar_evaluate_spherical_wave_sum.cc
@@ -32,100 +32,98 @@
#include "oskar_vector_types.h"
#include "oskar_helper.h"
-template<typename FP, typename FP2, typename FP4c>
-void oskar_evaluate_spherical_wave_sum(
- int num_points, const FP* theta,
- const FP* phi_x, const FP* phi_y, int l_max,
- const FP4c* alpha, FP4c* pattern)
-{
- FP2 Xp, Xt, Yp, Yt;
- make_zero2(Xp); make_zero2(Xt);
- make_zero2(Yp); make_zero2(Yt);
+template <typename FP, typename FP2, typename FP4c>
+void oskar_evaluate_spherical_wave_sum(int num_points, const FP* theta,
+ const FP* phi_x, const FP* phi_y,
+ int l_max, const FP4c* alpha,
+ FP4c* pattern) {
+ FP2 Xp, Xt, Yp, Yt;
+ make_zero2(Xp);
+ make_zero2(Xt);
+ make_zero2(Yp);
+ make_zero2(Yt);
- #pragma omp parallel for
- for (int i = 0; i < num_points; i++) {
- const FP phi_x_ = phi_x[i], theta_ = theta[i];
- const FP phi_y_ = phi_y[i];
- /* Propagate NAN. */
- if (phi_x_ != phi_x_) {
- Xp.x = Xp.y = Xt.x = Xt.y = phi_x_;
- Yp.x = Yp.y = Yt.x = Yt.y = phi_x_;
+#pragma omp parallel for
+ for (int i = 0; i < num_points; i++) {
+ const FP phi_x_ = phi_x[i], theta_ = theta[i];
+ const FP phi_y_ = phi_y[i];
+ /* Propagate NAN. */
+ if (phi_x_ != phi_x_) {
+ Xp.x = Xp.y = Xt.x = Xt.y = phi_x_;
+ Yp.x = Yp.y = Yt.x = Yt.y = phi_x_;
+ } else {
+ FP sin_t, cos_t;
+ oskar_sincos(theta_, &sin_t, &cos_t);
+ for (int l = 1; l <= l_max; ++l) {
+ const int ind0 = l * l - 1 + l;
+ const FP f_ = (2 * l + 1) / (4 * ((FP)M_PI) * l * (l + 1));
+ for (int abs_m = l; abs_m >= 0; --abs_m) {
+ FP p, pds, dpms, sin_p, cos_p;
+ oskar_legendre2(l, abs_m, cos_t, sin_t, p, pds, dpms);
+ if (abs_m == 0) {
+ sin_p = (FP)0;
+ cos_p = sqrt(f_);
+ const FP4c alpha_ = alpha[ind0];
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.a, alpha_.b, Xt,
+ Xp);
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.c, alpha_.d, Yt,
+ Yp);
+ } else {
+ FP d_fact = (FP)1, s_fact = (FP)1;
+ const int d_ = l - abs_m, s_ = l + abs_m;
+ d_fact = std::tgamma(d_ + 1);
+ s_fact = std::tgamma(s_ + 1);
+ const FP ff = f_ * d_fact / s_fact;
+ const FP nf = sqrt(ff);
+ const FP4c alpha_m = alpha[ind0 - abs_m];
+ const FP4c alpha_p = alpha[ind0 + abs_m];
+ p = -abs_m * phi_x_;
+ oskar_sincos(p, &sin_p, &cos_p);
+ sin_p *= nf;
+ cos_p *= nf;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, -abs_m, alpha_m.a,
+ alpha_m.b, Xt, Xp);
+ sin_p = -sin_p;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.a, alpha_p.b,
+ Xt, Xp);
+ p = -abs_m * phi_y_;
+ oskar_sincos(p, &sin_p, &cos_p);
+ sin_p *= nf;
+ cos_p *= nf;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, -abs_m, alpha_m.c,
+ alpha_m.d, Yt, Yp);
+ sin_p = -sin_p;
+ oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.c, alpha_p.d,
+ Yt, Yp);
+ }
}
- else {
- FP sin_t, cos_t;
- oskar_sincos(theta_, &sin_t, &cos_t);
- for (int l = 1; l <= l_max; ++l) {
- const int ind0 = l * l - 1 + l;
- const FP f_ = (2 * l + 1) / (4 * ((FP)M_PI) * l * (l + 1));
- for (int abs_m = l; abs_m >=0; --abs_m) {
- FP p, pds, dpms, sin_p, cos_p;
- oskar_legendre2(l, abs_m, cos_t, sin_t, p, pds, dpms);
- if (abs_m == 0) {
- sin_p = (FP)0; cos_p = sqrt(f_);
- const FP4c alpha_ = alpha[ind0];
- oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.a, alpha_.b, Xt, Xp);
- oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.c, alpha_.d, Yt, Yp);
- }
- else {
- FP d_fact = (FP)1, s_fact = (FP)1;
- const int d_ = l - abs_m, s_ = l + abs_m;
- d_fact = std::tgamma(d_ + 1);
- s_fact = std::tgamma(s_ + 1);
- const FP ff = f_ * d_fact / s_fact;
- const FP nf = sqrt(ff);
- const FP4c alpha_m = alpha[ind0 - abs_m];
- const FP4c alpha_p = alpha[ind0 + abs_m];
- p = -abs_m * phi_x_;
- oskar_sincos(p, &sin_p, &cos_p);
- sin_p *= nf; cos_p *= nf;
- oskar_sph_wave(pds, dpms, sin_p, cos_p, -abs_m, alpha_m.a, alpha_m.b, Xt, Xp);
- sin_p = -sin_p;
- oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.a, alpha_p.b, Xt, Xp);
- p = -abs_m * phi_y_;
- oskar_sincos(p, &sin_p, &cos_p);
- sin_p *= nf; cos_p *= nf;
- oskar_sph_wave(pds, dpms, sin_p, cos_p, -abs_m, alpha_m.c, alpha_m.d, Yt, Yp);
- sin_p = -sin_p;
- oskar_sph_wave(pds, dpms, sin_p, cos_p, abs_m, alpha_p.c, alpha_p.d, Yt, Yp);
- }
- }
- }
- }
- pattern[i].a = Xt;
- pattern[i].b = Xp;
- pattern[i].c = Yt;
- pattern[i].d = Yp;
+ }
}
+ pattern[i].a = Xt;
+ pattern[i].b = Xp;
+ pattern[i].c = Yt;
+ pattern[i].d = Yp;
+ }
}
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,
- std::complex<double>* pattern)
-{
- double4c* alpha_ptr = (double4c *) alpha;
- double4c* pattern_ptr = (double4c *) pattern;
+ 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) {
+ double4c* alpha_ptr = (double4c*)alpha;
+ double4c* pattern_ptr = (double4c*)pattern;
- oskar_evaluate_spherical_wave_sum<double, double2, double4c>(
- num_points, theta, phi_x, phi_y, l_max, alpha_ptr, pattern_ptr);
+ oskar_evaluate_spherical_wave_sum<double, double2, double4c>(
+ num_points, theta, phi_x, phi_y, l_max, alpha_ptr, pattern_ptr);
}
void oskar_evaluate_spherical_wave_sum_float(
- const int num_points,
- const float* theta,
- const float* phi_x,
- const float* phi_y,
- const int l_max,
- const std::complex<float>* alpha,
- std::complex<float>* pattern)
-{
- float4c* alpha_ptr = (float4c *) alpha;
- float4c* pattern_ptr = (float4c *) pattern;
+ const int num_points, const float* theta, const float* phi_x,
+ const float* phi_y, const int l_max, const std::complex<float>* alpha,
+ std::complex<float>* pattern) {
+ float4c* alpha_ptr = (float4c*)alpha;
+ float4c* pattern_ptr = (float4c*)pattern;
- oskar_evaluate_spherical_wave_sum<float, float2, float4c>(
- num_points, theta, phi_x, phi_y, l_max, alpha_ptr, pattern_ptr);
+ oskar_evaluate_spherical_wave_sum<float, float2, float4c>(
+ num_points, theta, phi_x, phi_y, l_max, alpha_ptr, pattern_ptr);
}
diff --git a/pystationresponse.cc b/pystationresponse.cc
old mode 100755
new mode 100644
index 4b8b2bb2199983fa48a3fc0b806edd40516ad493..960ab4d924855c90a8133351e5761ccfe978622b
--- a/pystationresponse.cc
+++ b/pystationresponse.cc
@@ -19,7 +19,6 @@
//#
//# $Id: pystationresponse.cc 33141 2015-12-16 15:10:19Z dijkema $
-
#include "ITRFDirection.h"
#include "LofarMetaDataUtil.h"
#include "Station.h"
@@ -48,439 +47,387 @@ using namespace casacore;
using namespace boost::python;
using namespace LOFAR::StationResponse;
-namespace LOFAR
-{
-namespace BBS
-{
- namespace
- {
- /*!
- * \brief Convert an ITRF position given as a StationResponse::vector3r_t
- * instance to a casacore::MPosition.
- */
- MPosition toMPositionITRF(const vector3r_t &position);
-
- /*!
- * \brief Convert a casacore::MPosition instance to a
- # StationResponse::vector3r_t instance.
- */
- vector3r_t fromMPosition(const MPosition &position);
-
- /*!
- * \brief Convert a casacore::MDirection instance to a
- * StationResponse::vector3r_t instance.
- */
- vector3r_t fromMDirection(const MDirection &direction);
-
- /*!
- * \brief Check if the specified column exists as a column of the specified
- * table.
- *
- * \param table The Table instance to check.
- * \param column The name of the column.
- */
- bool hasColumn(const Table &table, const string &column);
-
- /*!
- * \brief Check if the specified sub-table exists as a sub-table of the
- * specified table.
- *
- * \param table The Table instance to check.
- * \param name The name of the sub-table.
- */
- bool hasSubTable(const Table &table, const string &name);
-
- /*!
- * \brief Provide access to a sub-table by name.
- *
- * \param table The Table instance to which the sub-table is associated.
- * \param name The name of the sub-table.
- */
- Table getSubTable(const Table &table, const string &name);
-
- /*!
- * \brief Attempt to read the position of the observatory. If the
- * observatory position is unknown, the specified default position is
- * returned.
- *
- * \param ms MeasurementSet to read the observatory position from.
- * \param idObservation Identifier that determines of which observation the
- * observatory position should be read.
- * \param defaultPosition The position that will be returned if the
- * observatory position is unknown.
- */
- MPosition readObservatoryPosition(const MeasurementSet &ms,
- unsigned int idObservation, const MPosition &defaultPosition);
-
- /*!
- * \brief Read the phase reference direction.
- *
- * \param ms MeasurementSet to read the phase reference direction from.
- * \param idField Identifier of the field of which the phase reference
- * direction should be read.
- */
- MDirection readPhaseReference(const MeasurementSet &ms,
- unsigned int idField);
-
- /*!
- * \brief Read the station beam former reference direction.
- *
- * \param ms MeasurementSet to read the station beam former reference
- * direction from.
- * \param idField Identifier of the field of which the station beam former
- * reference direction should be read.
- */
- MDirection readDelayReference(const MeasurementSet &ms,
- unsigned int idField);
-
- /*!
- * \brief Read the station beam former reference direction.
- *
- * \param ms MeasurementSet to read the tile beam former reference direction
- * from.
- * \param idField Identifier of the field of which the tile beam former
- * reference direction should be read.
- */
- MDirection readTileReference(const MeasurementSet &ms,
- unsigned int idField);
- } //# unnamed namespace
-
- class PyStationResponse
- {
- public:
- PyStationResponse(const string& msName, bool inverse = false,
- bool useElementResponse = true, bool useArrayFactor = true,
- bool useChanFreq = false);
-
- // Get the software version.
- string version(const string& type) const;
-
- // Set the delay reference direction in radians, J2000. The delay reference
- // direction is the direction used by the station beamformer.
- void setRefDelay(double ra, double dec);
-
- // Get the delay reference direction in meters, ITRF. The delay reference
- // direction is the direction used by the station beamformer.
- ValueHolder getRefDelay(real_t time);
-
- // Set the tile reference direction in radians, J2000. The tile reference
- // direction is the direction used by the analog tile beamformer and is
- // relevant only for HBA observations.
- void setRefTile(double ra, double dec);
-
- // Get the tile reference direction in meters, ITRF. The delay reference
- // direction is the direction used by the analog tile beamformer and is
- // relevant only for HBA observations.
- ValueHolder getRefTile(real_t time);
-
- // Set the direction of interest in radians, J2000. Can and often will be
- // different than the delay and/or tile reference direction.
- void setDirection(double ra, double dec);
-
- // Get the direction of intereset in meters, ITRF.
- ValueHolder getDirection(real_t time);
-
- // Compute the LOFAR beam Jones matrices for the given time, station, and/or
- // channel.
- ValueHolder evaluate0(double time);
- ValueHolder evaluate1(double time, int station);
- ValueHolder evaluate2(double time, int station, int channel);
- ValueHolder evaluate3(double time, int station, double freq);
- ValueHolder evaluate4(double time, int station, double freq, const ValueHolder& direction, const ValueHolder& station0, const ValueHolder& tile0);
-
- private:
- Matrix<DComplex> evaluate_itrf(
- const Station::ConstPtr &station, double time, double freq, double freq0,
- const vector3r_t &direction, const vector3r_t &station0,
- const vector3r_t &tile0) const;
-
- Matrix<DComplex> evaluate(const Station::ConstPtr &station, double time,
- double freq, double freq0) const;
-
- Cube<DComplex> evaluate(const Station::ConstPtr &station, double time,
- const Vector<Double> &freq, const Vector<Double> &freq0) const;
-
- void invert(matrix22c_t &in) const;
- void invert(diag22c_t &in) const;
-
- //# Data members.
- bool itsInverse;
- bool itsUseElementResponse;
- bool itsUseArrayFactor;
- bool itsUseChanFreq;
- Vector<Station::Ptr> itsStations;
- Vector<Double> itsChanFreq;
- Vector<Double> itsRefFreq;
-
- vector3r_t itsRefPosition;
- ITRFDirection::Ptr itsRefDelay;
- ITRFDirection::Ptr itsRefTile;
-
- ITRFDirection::Ptr itsDirection;
- };
-
- PyStationResponse::PyStationResponse(const string &name, bool inverse,
- bool useElementResponse, bool useArrayFactor, bool useChanFreq)
+namespace LOFAR {
+namespace BBS {
+namespace {
+/*!
+ * \brief Convert an ITRF position given as a StationResponse::vector3r_t
+ * instance to a casacore::MPosition.
+ */
+MPosition toMPositionITRF(const vector3r_t &position);
+
+/*!
+ * \brief Convert a casacore::MPosition instance to a
+ # StationResponse::vector3r_t instance.
+ */
+vector3r_t fromMPosition(const MPosition &position);
+
+/*!
+ * \brief Convert a casacore::MDirection instance to a
+ * StationResponse::vector3r_t instance.
+ */
+vector3r_t fromMDirection(const MDirection &direction);
+
+/*!
+ * \brief Check if the specified column exists as a column of the specified
+ * table.
+ *
+ * \param table The Table instance to check.
+ * \param column The name of the column.
+ */
+bool hasColumn(const Table &table, const string &column);
+
+/*!
+ * \brief Check if the specified sub-table exists as a sub-table of the
+ * specified table.
+ *
+ * \param table The Table instance to check.
+ * \param name The name of the sub-table.
+ */
+bool hasSubTable(const Table &table, const string &name);
+
+/*!
+ * \brief Provide access to a sub-table by name.
+ *
+ * \param table The Table instance to which the sub-table is associated.
+ * \param name The name of the sub-table.
+ */
+Table getSubTable(const Table &table, const string &name);
+
+/*!
+ * \brief Attempt to read the position of the observatory. If the
+ * observatory position is unknown, the specified default position is
+ * returned.
+ *
+ * \param ms MeasurementSet to read the observatory position from.
+ * \param idObservation Identifier that determines of which observation the
+ * observatory position should be read.
+ * \param defaultPosition The position that will be returned if the
+ * observatory position is unknown.
+ */
+MPosition readObservatoryPosition(const MeasurementSet &ms,
+ unsigned int idObservation,
+ const MPosition &defaultPosition);
+
+/*!
+ * \brief Read the phase reference direction.
+ *
+ * \param ms MeasurementSet to read the phase reference direction from.
+ * \param idField Identifier of the field of which the phase reference
+ * direction should be read.
+ */
+MDirection readPhaseReference(const MeasurementSet &ms, unsigned int idField);
+
+/*!
+ * \brief Read the station beam former reference direction.
+ *
+ * \param ms MeasurementSet to read the station beam former reference
+ * direction from.
+ * \param idField Identifier of the field of which the station beam former
+ * reference direction should be read.
+ */
+MDirection readDelayReference(const MeasurementSet &ms, unsigned int idField);
+
+/*!
+ * \brief Read the station beam former reference direction.
+ *
+ * \param ms MeasurementSet to read the tile beam former reference direction
+ * from.
+ * \param idField Identifier of the field of which the tile beam former
+ * reference direction should be read.
+ */
+MDirection readTileReference(const MeasurementSet &ms, unsigned int idField);
+} // namespace
+
+class PyStationResponse {
+ public:
+ PyStationResponse(const string &msName, bool inverse = false,
+ bool useElementResponse = true, bool useArrayFactor = true,
+ bool useChanFreq = false);
+
+ // Get the software version.
+ string version(const string &type) const;
+
+ // Set the delay reference direction in radians, J2000. The delay reference
+ // direction is the direction used by the station beamformer.
+ void setRefDelay(double ra, double dec);
+
+ // Get the delay reference direction in meters, ITRF. The delay reference
+ // direction is the direction used by the station beamformer.
+ ValueHolder getRefDelay(real_t time);
+
+ // Set the tile reference direction in radians, J2000. The tile reference
+ // direction is the direction used by the analog tile beamformer and is
+ // relevant only for HBA observations.
+ void setRefTile(double ra, double dec);
+
+ // Get the tile reference direction in meters, ITRF. The delay reference
+ // direction is the direction used by the analog tile beamformer and is
+ // relevant only for HBA observations.
+ ValueHolder getRefTile(real_t time);
+
+ // Set the direction of interest in radians, J2000. Can and often will be
+ // different than the delay and/or tile reference direction.
+ void setDirection(double ra, double dec);
+
+ // Get the direction of intereset in meters, ITRF.
+ ValueHolder getDirection(real_t time);
+
+ // Compute the LOFAR beam Jones matrices for the given time, station, and/or
+ // channel.
+ ValueHolder evaluate0(double time);
+ ValueHolder evaluate1(double time, int station);
+ ValueHolder evaluate2(double time, int station, int channel);
+ ValueHolder evaluate3(double time, int station, double freq);
+ ValueHolder evaluate4(double time, int station, double freq,
+ const ValueHolder &direction,
+ const ValueHolder &station0, const ValueHolder &tile0);
+
+ private:
+ Matrix<DComplex> evaluate_itrf(const Station::ConstPtr &station, double time,
+ double freq, double freq0,
+ const vector3r_t &direction,
+ const vector3r_t &station0,
+ const vector3r_t &tile0) const;
+
+ Matrix<DComplex> evaluate(const Station::ConstPtr &station, double time,
+ double freq, double freq0) const;
+
+ Cube<DComplex> evaluate(const Station::ConstPtr &station, double time,
+ const Vector<Double> &freq,
+ const Vector<Double> &freq0) const;
+
+ void invert(matrix22c_t &in) const;
+ void invert(diag22c_t &in) const;
+
+ //# Data members.
+ bool itsInverse;
+ bool itsUseElementResponse;
+ bool itsUseArrayFactor;
+ bool itsUseChanFreq;
+ Vector<Station::Ptr> itsStations;
+ Vector<Double> itsChanFreq;
+ Vector<Double> itsRefFreq;
+
+ vector3r_t itsRefPosition;
+ ITRFDirection::Ptr itsRefDelay;
+ ITRFDirection::Ptr itsRefTile;
+
+ ITRFDirection::Ptr itsDirection;
+};
+
+PyStationResponse::PyStationResponse(const string &name, bool inverse,
+ bool useElementResponse,
+ bool useArrayFactor, bool useChanFreq)
: itsInverse(inverse),
itsUseElementResponse(useElementResponse),
itsUseArrayFactor(useArrayFactor),
- itsUseChanFreq(useChanFreq)
- {
- MeasurementSet ms(name);
-
- // Read spectral window id.
- const unsigned int idDataDescription = 0;
- ROMSDataDescColumns desc(ms.dataDescription());
- assert(desc.nrow() > idDataDescription);
- assert(!desc.flagRow()(idDataDescription));
-
- // Read the spectral information.
- const unsigned int idWindow = desc.spectralWindowId()(idDataDescription);
- ROMSSpWindowColumns window(ms.spectralWindow());
- assert(window.nrow() > idWindow);
- assert(!window.flagRow()(idWindow));
-
- itsChanFreq = window.chanFreq()(idWindow);
- itsRefFreq = Vector<Double>(itsChanFreq.size(),
- window.refFrequency()(idWindow));
-
- // Read the station information.
- ROMSAntennaColumns antenna(ms.antenna());
- itsStations.resize(antenna.nrow());
- for(unsigned int i = 0; i < antenna.nrow(); ++i)
- {
- itsStations(i) = readStation(ms, i);
- }
+ itsUseChanFreq(useChanFreq) {
+ MeasurementSet ms(name);
+
+ // Read spectral window id.
+ const unsigned int idDataDescription = 0;
+ ROMSDataDescColumns desc(ms.dataDescription());
+ assert(desc.nrow() > idDataDescription);
+ assert(!desc.flagRow()(idDataDescription));
+
+ // Read the spectral information.
+ const unsigned int idWindow = desc.spectralWindowId()(idDataDescription);
+ ROMSSpWindowColumns window(ms.spectralWindow());
+ assert(window.nrow() > idWindow);
+ assert(!window.flagRow()(idWindow));
+
+ itsChanFreq = window.chanFreq()(idWindow);
+ itsRefFreq =
+ Vector<Double>(itsChanFreq.size(), window.refFrequency()(idWindow));
+
+ // Read the station information.
+ ROMSAntennaColumns antenna(ms.antenna());
+ itsStations.resize(antenna.nrow());
+ for (unsigned int i = 0; i < antenna.nrow(); ++i) {
+ itsStations(i) = readStation(ms, i);
+ }
- // Read observatory position. If unknown, default to the position of the
- // first station.
- unsigned int idObservation = 0;
- MPosition refPosition = readObservatoryPosition(ms, idObservation,
- toMPositionITRF(itsStations(0)->position()));
- itsRefPosition = fromMPosition(MPosition::Convert(refPosition,
- MPosition::ITRF)());
-
- // Read the reference directions.
- unsigned int idField = 0;
- itsRefDelay.reset(new ITRFDirection(itsRefPosition,
+ // Read observatory position. If unknown, default to the position of the
+ // first station.
+ unsigned int idObservation = 0;
+ MPosition refPosition = readObservatoryPosition(
+ ms, idObservation, toMPositionITRF(itsStations(0)->position()));
+ itsRefPosition =
+ fromMPosition(MPosition::Convert(refPosition, MPosition::ITRF)());
+
+ // Read the reference directions.
+ unsigned int idField = 0;
+ itsRefDelay.reset(new ITRFDirection(
+ itsRefPosition,
fromMDirection(MDirection::Convert(readDelayReference(ms, idField),
- MDirection::J2000)())));
+ MDirection::J2000)())));
- itsRefTile.reset(new ITRFDirection(itsRefPosition,
+ itsRefTile.reset(new ITRFDirection(
+ itsRefPosition,
fromMDirection(MDirection::Convert(readTileReference(ms, idField),
- MDirection::J2000)())));
+ MDirection::J2000)())));
- itsDirection.reset(new ITRFDirection(itsRefPosition,
+ itsDirection.reset(new ITRFDirection(
+ itsRefPosition,
fromMDirection(MDirection::Convert(readPhaseReference(ms, idField),
- MDirection::J2000)())));
- }
+ MDirection::J2000)())));
+}
- string PyStationResponse::version(const string& type) const
- {
-// return Version::getInfo<pystationresponseVersion>("stationresponse", type);
- return "0.1";
- }
+string PyStationResponse::version(const string &type) const {
+ // return Version::getInfo<pystationresponseVersion>("stationresponse",
+ // type);
+ return "0.1";
+}
- void PyStationResponse::setRefDelay(double ra, double dec)
- {
- vector2r_t direction = {{ra, dec}};
- itsRefDelay.reset(new ITRFDirection(itsRefPosition, direction));
- }
+void PyStationResponse::setRefDelay(double ra, double dec) {
+ vector2r_t direction = {{ra, dec}};
+ itsRefDelay.reset(new ITRFDirection(itsRefPosition, direction));
+}
- ValueHolder PyStationResponse::getRefDelay(real_t time)
- {
- vector3r_t refDelay=itsRefDelay->at(time);
- Vector<Double> result(3);
- result(0)=refDelay[0]; result(1)=refDelay[1]; result(2)=refDelay[2];
+ValueHolder PyStationResponse::getRefDelay(real_t time) {
+ vector3r_t refDelay = itsRefDelay->at(time);
+ Vector<Double> result(3);
+ result(0) = refDelay[0];
+ result(1) = refDelay[1];
+ result(2) = refDelay[2];
- return ValueHolder(result);
- }
+ return ValueHolder(result);
+}
- void PyStationResponse::setRefTile(double ra, double dec)
- {
- vector2r_t direction = {{ra, dec}};
- itsRefTile.reset(new ITRFDirection(itsRefPosition, direction));
- }
+void PyStationResponse::setRefTile(double ra, double dec) {
+ vector2r_t direction = {{ra, dec}};
+ itsRefTile.reset(new ITRFDirection(itsRefPosition, direction));
+}
- ValueHolder PyStationResponse::getRefTile(real_t time)
- {
- vector3r_t refTile=itsRefTile->at(time);
- Vector<Double> result(3);
- result(0)=refTile[0]; result(1)=refTile[1]; result(2)=refTile[2];
+ValueHolder PyStationResponse::getRefTile(real_t time) {
+ vector3r_t refTile = itsRefTile->at(time);
+ Vector<Double> result(3);
+ result(0) = refTile[0];
+ result(1) = refTile[1];
+ result(2) = refTile[2];
- return ValueHolder(result);
- }
+ return ValueHolder(result);
+}
- void PyStationResponse::setDirection(double ra, double dec)
- {
- vector2r_t direction = {{ra, dec}};
- itsDirection.reset(new ITRFDirection(itsRefPosition, direction));
- }
+void PyStationResponse::setDirection(double ra, double dec) {
+ vector2r_t direction = {{ra, dec}};
+ itsDirection.reset(new ITRFDirection(itsRefPosition, direction));
+}
- ValueHolder PyStationResponse::getDirection(real_t time)
- {
- vector3r_t direction=itsDirection->at(time);
- Vector<Double> result(3);
- result(0)=direction[0]; result(1)=direction[1]; result(2)=direction[2];
+ValueHolder PyStationResponse::getDirection(real_t time) {
+ vector3r_t direction = itsDirection->at(time);
+ Vector<Double> result(3);
+ result(0) = direction[0];
+ result(1) = direction[1];
+ result(2) = direction[2];
- return ValueHolder(result);
- }
+ return ValueHolder(result);
+}
- ValueHolder PyStationResponse::evaluate0(double time)
- {
- Array<DComplex> result(IPosition(4, 2, 2, itsChanFreq.size(),
- itsStations.size()));
-
- for(unsigned int i = 0; i < itsStations.size(); ++i)
- {
- IPosition start(4, 0, 0, 0, i);
- IPosition end(4, 1, 1, itsChanFreq.size() - 1, i);
- Cube<DComplex> slice = result(start, end).nonDegenerate();
- if(itsUseChanFreq)
- {
- slice = evaluate(itsStations(i), time, itsChanFreq, itsChanFreq);
- }
- else
- {
- slice = evaluate(itsStations(i), time, itsChanFreq, itsRefFreq);
- }
+ValueHolder PyStationResponse::evaluate0(double time) {
+ Array<DComplex> result(
+ IPosition(4, 2, 2, itsChanFreq.size(), itsStations.size()));
+
+ for (unsigned int i = 0; i < itsStations.size(); ++i) {
+ IPosition start(4, 0, 0, 0, i);
+ IPosition end(4, 1, 1, itsChanFreq.size() - 1, i);
+ Cube<DComplex> slice = result(start, end).nonDegenerate();
+ if (itsUseChanFreq) {
+ slice = evaluate(itsStations(i), time, itsChanFreq, itsChanFreq);
+ } else {
+ slice = evaluate(itsStations(i), time, itsChanFreq, itsRefFreq);
}
-
- return ValueHolder(result);
}
- ValueHolder PyStationResponse::evaluate1(double time, int station)
- {
-// assertSTR(station >= 0 && static_cast<size_t>(station)
-// < itsStations.size(), "invalid station number: " << station);
+ return ValueHolder(result);
+}
- if(itsUseChanFreq)
- {
- return ValueHolder(evaluate(itsStations(station), time, itsChanFreq,
- itsChanFreq));
- }
+ValueHolder PyStationResponse::evaluate1(double time, int station) {
+ // assertSTR(station >= 0 && static_cast<size_t>(station)
+ // < itsStations.size(), "invalid station number: " << station);
- return ValueHolder(evaluate(itsStations(station), time, itsChanFreq,
- itsRefFreq));
+ if (itsUseChanFreq) {
+ return ValueHolder(
+ evaluate(itsStations(station), time, itsChanFreq, itsChanFreq));
}
- ValueHolder PyStationResponse::evaluate2(double time, int station,
- int channel)
- {
-// assertSTR(station >= 0 && static_cast<size_t>(station)
-// < itsStations.size(), "invalid station number: " << station);
-// assertSTR(channel >= 0 && static_cast<size_t>(channel)
-// < itsChanFreq.size(), "invalid channel number: " << channel);
-
+ return ValueHolder(
+ evaluate(itsStations(station), time, itsChanFreq, itsRefFreq));
+}
- double freq = itsChanFreq(channel);
- if(itsUseChanFreq)
- {
- return ValueHolder(evaluate(itsStations(station), time, freq, freq));
- }
+ValueHolder PyStationResponse::evaluate2(double time, int station,
+ int channel) {
+ // assertSTR(station >= 0 && static_cast<size_t>(station)
+ // < itsStations.size(), "invalid station number: " << station);
+ // assertSTR(channel >= 0 && static_cast<size_t>(channel)
+ // < itsChanFreq.size(), "invalid channel number: " << channel);
- double freq0 = itsRefFreq(channel);
- return ValueHolder(evaluate(itsStations(station), time, freq, freq0));
+ double freq = itsChanFreq(channel);
+ if (itsUseChanFreq) {
+ return ValueHolder(evaluate(itsStations(station), time, freq, freq));
}
- ValueHolder PyStationResponse::evaluate3(double time, int station,
- double freq)
- {
-// assertSTR(station >= 0 && static_cast<size_t>(station)
-// < itsStations.size(), "invalid station number: " << station);
+ double freq0 = itsRefFreq(channel);
+ return ValueHolder(evaluate(itsStations(station), time, freq, freq0));
+}
- if(itsUseChanFreq)
- {
- return ValueHolder(evaluate(itsStations(station), time, freq, freq));
- }
+ValueHolder PyStationResponse::evaluate3(double time, int station,
+ double freq) {
+ // assertSTR(station >= 0 && static_cast<size_t>(station)
+ // < itsStations.size(), "invalid station number: " << station);
- double freq0 = itsRefFreq(0);
- return ValueHolder(evaluate(itsStations(station), time, freq, freq0));
+ if (itsUseChanFreq) {
+ return ValueHolder(evaluate(itsStations(station), time, freq, freq));
}
- ValueHolder PyStationResponse::evaluate4(double time, int station, double freq, const ValueHolder& vh_direction, const ValueHolder& vh_station0, const ValueHolder& vh_tile0)
- {
- assert (vh_direction.dataType() == TpArrayDouble);
- assert (vh_station0.dataType() == TpArrayDouble);
- assert (vh_tile0.dataType() == TpArrayDouble);
- Array<Double> arr_dir(vh_direction.asArrayDouble());
- Array<Double> st0_dir(vh_station0.asArrayDouble());
- Array<Double> tile_dir(vh_tile0.asArrayDouble());
- vector3r_t direction={{arr_dir.data()[0],arr_dir.data()[1],arr_dir.data()[2]}};
- vector3r_t station0 ={{st0_dir.data()[0],st0_dir.data()[1],st0_dir.data()[2]}};
- vector3r_t tile0 ={{tile_dir.data()[0],tile_dir.data()[1],tile_dir.data()[2]}};
-
- if(itsUseChanFreq)
- {
- return ValueHolder(evaluate_itrf(itsStations(station), time, freq, freq,
- direction, station0, tile0));
- }
+ double freq0 = itsRefFreq(0);
+ return ValueHolder(evaluate(itsStations(station), time, freq, freq0));
+}
- double freq0 = itsRefFreq(0);
- return ValueHolder(evaluate_itrf(itsStations(station), time, freq, freq0,
+ValueHolder PyStationResponse::evaluate4(double time, int station, double freq,
+ const ValueHolder &vh_direction,
+ const ValueHolder &vh_station0,
+ const ValueHolder &vh_tile0) {
+ assert(vh_direction.dataType() == TpArrayDouble);
+ assert(vh_station0.dataType() == TpArrayDouble);
+ assert(vh_tile0.dataType() == TpArrayDouble);
+ Array<Double> arr_dir(vh_direction.asArrayDouble());
+ Array<Double> st0_dir(vh_station0.asArrayDouble());
+ Array<Double> tile_dir(vh_tile0.asArrayDouble());
+ vector3r_t direction = {
+ {arr_dir.data()[0], arr_dir.data()[1], arr_dir.data()[2]}};
+ vector3r_t station0 = {
+ {st0_dir.data()[0], st0_dir.data()[1], st0_dir.data()[2]}};
+ vector3r_t tile0 = {
+ {tile_dir.data()[0], tile_dir.data()[1], tile_dir.data()[2]}};
+
+ if (itsUseChanFreq) {
+ return ValueHolder(evaluate_itrf(itsStations(station), time, freq, freq,
direction, station0, tile0));
}
- Cube<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
- double time, const Vector<Double> &freq, const Vector<Double> &freq0) const
- {
- Cube<DComplex> result(2, 2, freq.size(), 0.0);
- if(itsUseArrayFactor)
- {
- vector3r_t direction = itsDirection->at(time);
- vector3r_t station0 = itsRefDelay->at(time);
- vector3r_t tile0 = itsRefTile->at(time);
-
- if(itsUseElementResponse)
- {
- for(unsigned int i = 0; i < freq.size(); ++i)
- {
- matrix22c_t response = station->response(time, freq(i), direction,
- freq0(i), station0, tile0);
-
- if(itsInverse)
- {
- invert(response);
- }
-
- result(0, 0, i) = response[0][0];
- result(1, 0, i) = response[0][1];
- result(0, 1, i) = response[1][0];
- result(1, 1, i) = response[1][1];
- }
- }
- else
- {
- for(unsigned int i = 0; i < freq.size(); ++i)
- {
- diag22c_t af = station->arrayFactor(time, freq(i), direction,
- freq0(i), station0, tile0);
-
- if(itsInverse)
- {
- invert(af);
- }
-
- result(0, 0, i) = af[0];
- result(1, 1, i) = af[1];
- }
- }
- }
- else if(itsUseElementResponse)
- {
- // For a station with multiple antenna fields, need to select for which
- // field the element response will be evaluated. Here the first field of the
- // station is always selected.
- AntennaField::ConstPtr field = *station->beginFields();
-
- vector3r_t direction = itsDirection->at(time);
- for(unsigned int i = 0; i < freq.size(); ++i)
- {
- matrix22c_t response = field->elementResponse(time, freq(i),
- direction);
-
- if(itsInverse)
- {
+ double freq0 = itsRefFreq(0);
+ return ValueHolder(evaluate_itrf(itsStations(station), time, freq, freq0,
+ direction, station0, tile0));
+}
+
+Cube<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
+ double time,
+ const Vector<Double> &freq,
+ const Vector<Double> &freq0) const {
+ Cube<DComplex> result(2, 2, freq.size(), 0.0);
+ if (itsUseArrayFactor) {
+ vector3r_t direction = itsDirection->at(time);
+ vector3r_t station0 = itsRefDelay->at(time);
+ vector3r_t tile0 = itsRefTile->at(time);
+
+ if (itsUseElementResponse) {
+ for (unsigned int i = 0; i < freq.size(); ++i) {
+ matrix22c_t response = station->response(time, freq(i), direction,
+ freq0(i), station0, tile0);
+
+ if (itsInverse) {
invert(response);
}
@@ -489,68 +436,59 @@ namespace BBS
result(0, 1, i) = response[1][0];
result(1, 1, i) = response[1][1];
}
- }
- else
- {
- for(unsigned int i = 0; i < freq.size(); ++i)
- {
- result(0, 0, i) = 1.0;
- result(1, 1, i) = 1.0;
- }
- }
-
- return result;
- }
+ } else {
+ for (unsigned int i = 0; i < freq.size(); ++i) {
+ diag22c_t af = station->arrayFactor(time, freq(i), direction, freq0(i),
+ station0, tile0);
- Matrix<DComplex> PyStationResponse::evaluate_itrf(
- const Station::ConstPtr &station, double time, double freq, double freq0,
- const vector3r_t &direction, const vector3r_t &station0,
- const vector3r_t &tile0) const
- {
- Matrix<DComplex> result(2, 2, 0.0);
- if(itsUseArrayFactor)
- {
- if(itsUseElementResponse)
- {
- matrix22c_t response = station->response(time, freq, direction, freq0,
- station0, tile0);
-
- if(itsInverse)
- {
- invert(response);
+ if (itsInverse) {
+ invert(af);
}
- result(0, 0) = response[0][0];
- result(1, 0) = response[0][1];
- result(0, 1) = response[1][0];
- result(1, 1) = response[1][1];
+ result(0, 0, i) = af[0];
+ result(1, 1, i) = af[1];
}
- else
- {
- diag22c_t af = station->arrayFactor(time, freq, direction, freq0,
- station0, tile0);
+ }
+ } else if (itsUseElementResponse) {
+ // For a station with multiple antenna fields, need to select for which
+ // field the element response will be evaluated. Here the first field of the
+ // station is always selected.
+ AntennaField::ConstPtr field = *station->beginFields();
- if(itsInverse)
- {
- invert(af);
- }
+ vector3r_t direction = itsDirection->at(time);
+ for (unsigned int i = 0; i < freq.size(); ++i) {
+ matrix22c_t response = field->elementResponse(time, freq(i), direction);
- result(0, 0) = af[0];
- result(1, 1) = af[1];
+ if (itsInverse) {
+ invert(response);
}
+
+ result(0, 0, i) = response[0][0];
+ result(1, 0, i) = response[0][1];
+ result(0, 1, i) = response[1][0];
+ result(1, 1, i) = response[1][1];
}
- else if(itsUseElementResponse)
- {
- // For a station with multiple antenna fields, need to select for which
- // field the element response will be evaluated. Here the first field of
- // the station is always selected.
- AntennaField::ConstPtr field = *station->beginFields();
-
- matrix22c_t response = field->elementResponse(time, freq,
- direction);
-
- if(itsInverse)
- {
+ } else {
+ for (unsigned int i = 0; i < freq.size(); ++i) {
+ result(0, 0, i) = 1.0;
+ result(1, 1, i) = 1.0;
+ }
+ }
+
+ return result;
+}
+
+Matrix<DComplex> PyStationResponse::evaluate_itrf(
+ const Station::ConstPtr &station, double time, double freq, double freq0,
+ const vector3r_t &direction, const vector3r_t &station0,
+ const vector3r_t &tile0) const {
+ Matrix<DComplex> result(2, 2, 0.0);
+ if (itsUseArrayFactor) {
+ if (itsUseElementResponse) {
+ matrix22c_t response =
+ station->response(time, freq, direction, freq0, station0, tile0);
+
+ if (itsInverse) {
invert(response);
}
@@ -558,198 +496,203 @@ namespace BBS
result(1, 0) = response[0][1];
result(0, 1) = response[1][0];
result(1, 1) = response[1][1];
+ } else {
+ diag22c_t af =
+ station->arrayFactor(time, freq, direction, freq0, station0, tile0);
+
+ if (itsInverse) {
+ invert(af);
+ }
+
+ result(0, 0) = af[0];
+ result(1, 1) = af[1];
}
- else
- {
- result(0, 0) = 1.0;
- result(1, 1) = 1.0;
+ } else if (itsUseElementResponse) {
+ // For a station with multiple antenna fields, need to select for which
+ // field the element response will be evaluated. Here the first field of
+ // the station is always selected.
+ AntennaField::ConstPtr field = *station->beginFields();
+
+ matrix22c_t response = field->elementResponse(time, freq, direction);
+
+ if (itsInverse) {
+ invert(response);
}
- return result;
+ result(0, 0) = response[0][0];
+ result(1, 0) = response[0][1];
+ result(0, 1) = response[1][0];
+ result(1, 1) = response[1][1];
+ } else {
+ result(0, 0) = 1.0;
+ result(1, 1) = 1.0;
}
- Matrix<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
- double time, double freq, double freq0) const
- {
- vector3r_t direction;
- vector3r_t station0;
- vector3r_t tile0;
- if (itsUseArrayFactor) {
- direction = itsDirection->at(time);
- station0 = itsRefDelay->at(time);
- tile0 = itsRefTile->at(time);
- } else if (itsUseElementResponse) {
- direction = itsDirection->at(time);
- }
- return evaluate_itrf(station, time, freq, freq0, direction, station0, tile0);
+ return result;
+}
+
+Matrix<DComplex> PyStationResponse::evaluate(const Station::ConstPtr &station,
+ double time, double freq,
+ double freq0) const {
+ vector3r_t direction;
+ vector3r_t station0;
+ vector3r_t tile0;
+ if (itsUseArrayFactor) {
+ direction = itsDirection->at(time);
+ station0 = itsRefDelay->at(time);
+ tile0 = itsRefTile->at(time);
+ } else if (itsUseElementResponse) {
+ direction = itsDirection->at(time);
}
+ return evaluate_itrf(station, time, freq, freq0, direction, station0, tile0);
+}
- void PyStationResponse::invert(matrix22c_t &in) const
- {
- complex_t invDet = 1.0 / (in[0][0] * in[1][1] - in[0][1] * in[1][0]);
+void PyStationResponse::invert(matrix22c_t &in) const {
+ complex_t invDet = 1.0 / (in[0][0] * in[1][1] - in[0][1] * in[1][0]);
- complex_t tmp = in[1][1];
- in[1][1] = in[0][0];
- in[0][0] = tmp;
+ complex_t tmp = in[1][1];
+ in[1][1] = in[0][0];
+ in[0][0] = tmp;
- in[0][0] *= invDet;
- in[0][1] *= -invDet;
- in[1][0] *= -invDet;
- in[1][1] *= invDet;
- }
+ in[0][0] *= invDet;
+ in[0][1] *= -invDet;
+ in[1][0] *= -invDet;
+ in[1][1] *= invDet;
+}
- void PyStationResponse::invert(diag22c_t &in) const
- {
- DComplex invDet = 1.0 / (in[0] * in[1]);
- DComplex tmp = in[1];
- in[1] = in[0];
- in[0] = tmp;
+void PyStationResponse::invert(diag22c_t &in) const {
+ DComplex invDet = 1.0 / (in[0] * in[1]);
+ DComplex tmp = in[1];
+ in[1] = in[0];
+ in[0] = tmp;
- in[0] *= invDet;
- in[1] *= invDet;
- }
+ in[0] *= invDet;
+ in[1] *= invDet;
+}
- // Now define the interface in Boost-Python.
- void pystationresponse()
- {
- class_<PyStationResponse> ("StationResponse",
- init<std::string, bool, bool, bool, bool>())
- .def ("version", &PyStationResponse::version,
- (boost::python::arg("type")="other"))
- .def ("setRefDelay", &PyStationResponse::setRefDelay,
- (boost::python::arg("ra"), boost::python::arg("dec")))
- .def ("getRefDelay", &PyStationResponse::getRefDelay,
- (boost::python::arg("time")))
- .def ("setRefTile", &PyStationResponse::setRefTile,
- (boost::python::arg("ra"), boost::python::arg("dec")))
- .def ("getRefTile", &PyStationResponse::getRefTile,
- (boost::python::arg("time")))
- .def ("setDirection", &PyStationResponse::setDirection,
- (boost::python::arg("ra"), boost::python::arg("dec")))
- .def ("getDirection", &PyStationResponse::getDirection,
- (boost::python::arg("time")))
- .def ("evaluate0", &PyStationResponse::evaluate0,
- (boost::python::arg("time")))
- .def ("evaluate1", &PyStationResponse::evaluate1,
- (boost::python::arg("time"), boost::python::arg("station")))
- .def ("evaluate2", &PyStationResponse::evaluate2,
- (boost::python::arg("time"), boost::python::arg("station"),
- boost::python::arg("channel")))
- .def ("evaluate3", &PyStationResponse::evaluate3,
- (boost::python::arg("time"), boost::python::arg("station"),
- boost::python::arg("freq")))
- .def ("evaluate4", &PyStationResponse::evaluate4,
- (boost::python::arg("time"), boost::python::arg("station"),
- boost::python::arg("freq"), boost::python::arg("direction"),
- boost::python::arg("station0"), boost::python::arg("tile0")))
- ;
- }
+// Now define the interface in Boost-Python.
+void pystationresponse() {
+ class_<PyStationResponse>("StationResponse",
+ init<std::string, bool, bool, bool, bool>())
+ .def("version", &PyStationResponse::version,
+ (boost::python::arg("type") = "other"))
+ .def("setRefDelay", &PyStationResponse::setRefDelay,
+ (boost::python::arg("ra"), boost::python::arg("dec")))
+ .def("getRefDelay", &PyStationResponse::getRefDelay,
+ (boost::python::arg("time")))
+ .def("setRefTile", &PyStationResponse::setRefTile,
+ (boost::python::arg("ra"), boost::python::arg("dec")))
+ .def("getRefTile", &PyStationResponse::getRefTile,
+ (boost::python::arg("time")))
+ .def("setDirection", &PyStationResponse::setDirection,
+ (boost::python::arg("ra"), boost::python::arg("dec")))
+ .def("getDirection", &PyStationResponse::getDirection,
+ (boost::python::arg("time")))
+ .def("evaluate0", &PyStationResponse::evaluate0,
+ (boost::python::arg("time")))
+ .def("evaluate1", &PyStationResponse::evaluate1,
+ (boost::python::arg("time"), boost::python::arg("station")))
+ .def("evaluate2", &PyStationResponse::evaluate2,
+ (boost::python::arg("time"), boost::python::arg("station"),
+ boost::python::arg("channel")))
+ .def("evaluate3", &PyStationResponse::evaluate3,
+ (boost::python::arg("time"), boost::python::arg("station"),
+ boost::python::arg("freq")))
+ .def("evaluate4", &PyStationResponse::evaluate4,
+ (boost::python::arg("time"), boost::python::arg("station"),
+ boost::python::arg("freq"), boost::python::arg("direction"),
+ boost::python::arg("station0"), boost::python::arg("tile0")));
+}
- namespace
- {
- MPosition toMPositionITRF(const vector3r_t &position)
- {
- MVPosition mvITRF(position[0], position[1], position[2]);
- return MPosition(mvITRF, MPosition::ITRF);
- }
+namespace {
+MPosition toMPositionITRF(const vector3r_t &position) {
+ MVPosition mvITRF(position[0], position[1], position[2]);
+ return MPosition(mvITRF, MPosition::ITRF);
+}
- vector3r_t fromMPosition(const MPosition &position)
- {
- MVPosition mvPosition = position.getValue();
- vector3r_t result = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
- return result;
- }
+vector3r_t fromMPosition(const MPosition &position) {
+ MVPosition mvPosition = position.getValue();
+ vector3r_t result = {{mvPosition(0), mvPosition(1), mvPosition(2)}};
+ return result;
+}
- vector3r_t fromMDirection(const MDirection &direction)
- {
- MVDirection mvDirection = direction.getValue();
- vector3r_t result = {{mvDirection(0), mvDirection(1), mvDirection(2)}};
- return result;
- }
+vector3r_t fromMDirection(const MDirection &direction) {
+ MVDirection mvDirection = direction.getValue();
+ vector3r_t result = {{mvDirection(0), mvDirection(1), mvDirection(2)}};
+ return result;
+}
- bool hasColumn(const Table &table, const string &column)
- {
- return table.tableDesc().isColumn(column);
- }
+bool hasColumn(const Table &table, const string &column) {
+ return table.tableDesc().isColumn(column);
+}
- bool hasSubTable(const Table &table, const string &name)
- {
- return table.keywordSet().isDefined(name);
- }
+bool hasSubTable(const Table &table, const string &name) {
+ return table.keywordSet().isDefined(name);
+}
- Table getSubTable(const Table &table, const string &name)
- {
- return table.keywordSet().asTable(name);
- }
+Table getSubTable(const Table &table, const string &name) {
+ return table.keywordSet().asTable(name);
+}
- MPosition readObservatoryPosition(const MeasurementSet &ms,
- unsigned int idObservation, const MPosition &defaultPosition)
- {
- // Get the instrument position in ITRF coordinates, or use the centroid
- // of the station positions if the instrument position is unknown.
- ROMSObservationColumns observation(ms.observation());
- assert(observation.nrow() > idObservation);
- assert(!observation.flagRow()(idObservation));
-
- // Read observatory name and try to look-up its position.
- const string observatory = observation.telescopeName()(idObservation);
-
- // Look-up observatory position, default to specified default position.
- MPosition position(defaultPosition);
- MeasTable::Observatory(position, observatory);
- return position;
- }
+MPosition readObservatoryPosition(const MeasurementSet &ms,
+ unsigned int idObservation,
+ const MPosition &defaultPosition) {
+ // Get the instrument position in ITRF coordinates, or use the centroid
+ // of the station positions if the instrument position is unknown.
+ ROMSObservationColumns observation(ms.observation());
+ assert(observation.nrow() > idObservation);
+ assert(!observation.flagRow()(idObservation));
+
+ // Read observatory name and try to look-up its position.
+ const string observatory = observation.telescopeName()(idObservation);
+
+ // Look-up observatory position, default to specified default position.
+ MPosition position(defaultPosition);
+ MeasTable::Observatory(position, observatory);
+ return position;
+}
- MDirection readPhaseReference(const MeasurementSet &ms,
- unsigned int idField)
- {
- ROMSFieldColumns field(ms.field());
- assert(field.nrow() > idField);
- assert(!field.flagRow()(idField));
+MDirection readPhaseReference(const MeasurementSet &ms, unsigned int idField) {
+ ROMSFieldColumns field(ms.field());
+ assert(field.nrow() > idField);
+ assert(!field.flagRow()(idField));
- return field.phaseDirMeas(idField);
- }
+ return field.phaseDirMeas(idField);
+}
- MDirection readDelayReference(const MeasurementSet &ms,
- unsigned int idField)
- {
- ROMSFieldColumns field(ms.field());
- assert(field.nrow() > idField);
- assert(!field.flagRow()(idField));
+MDirection readDelayReference(const MeasurementSet &ms, unsigned int idField) {
+ ROMSFieldColumns field(ms.field());
+ assert(field.nrow() > idField);
+ assert(!field.flagRow()(idField));
- return field.delayDirMeas(idField);
- }
+ return field.delayDirMeas(idField);
+}
- MDirection readTileReference(const MeasurementSet &ms,
- unsigned int idField)
- {
- // The MeasurementSet class does not support LOFAR specific columns, so we
- // use ROArrayMeasColumn to read the tile beam reference direction.
- Table tab_field = getSubTable(ms, "FIELD");
-
- static const String columnName = "LOFAR_TILE_BEAM_DIR";
- if(hasColumn(tab_field, columnName))
- {
- ROArrayMeasColumn<MDirection> c_direction(tab_field, columnName);
- if(c_direction.isDefined(idField))
- {
- return c_direction(idField)(IPosition(1, 0));
- }
- }
+MDirection readTileReference(const MeasurementSet &ms, unsigned int idField) {
+ // The MeasurementSet class does not support LOFAR specific columns, so we
+ // use ROArrayMeasColumn to read the tile beam reference direction.
+ Table tab_field = getSubTable(ms, "FIELD");
- // By default, the tile beam reference direction is assumed to be equal
- // to the station beam reference direction (for backward compatibility,
- // and for non-HBA measurements).
- return readDelayReference(ms, idField);
+ static const String columnName = "LOFAR_TILE_BEAM_DIR";
+ if (hasColumn(tab_field, columnName)) {
+ ROArrayMeasColumn<MDirection> c_direction(tab_field, columnName);
+ if (c_direction.isDefined(idField)) {
+ return c_direction(idField)(IPosition(1, 0));
}
- } //# unnamed namespace
+ }
+
+ // By default, the tile beam reference direction is assumed to be equal
+ // to the station beam reference direction (for backward compatibility,
+ // and for non-HBA measurements).
+ return readDelayReference(ms, idField);
+}
+} // namespace
-} //# namespace BBS
-} //# namespace LOFAR
+} // namespace BBS
+} // namespace LOFAR
// Define the python module itself.
-BOOST_PYTHON_MODULE(_stationresponse)
-{
+BOOST_PYTHON_MODULE(_stationresponse) {
casacore::python::register_convert_excp();
casacore::python::register_convert_basicdata();
casacore::python::register_convert_casa_valueholder();
diff --git a/scripts/clang-format-check.sh b/scripts/clang-format-check.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a750894f39e445da27378ccdc9a0fdeab1cd03ba
--- /dev/null
+++ b/scripts/clang-format-check.sh
@@ -0,0 +1,36 @@
+#!/bin/sh
+# Do check on clang-format, script is inspired on:
+# - https://dev.to/10xlearner/formatting-cpp-c-javascript-and-other-stuff-2pof
+# - https://gitlab.freedesktop.org/monado/monado/-/blob/master/scripts/format-and-spellcheck.sh
+#
+# Author: Jakob Maljaars
+# Email: jakob.maljaars_@_stcorp.nl
+#
+# To hook this script to pre-commit include the line "./scripts/clang-format-check.sh" to .git/hooks/pre-commit
+# and make sure pre-commit is an executable shell script.
+
+set -e
+PATCH_NAME=clang-fixes.diff
+
+echo "Running clang-format..."
+echo
+
+# Run clang-format from run-clang-format.sh
+SCRIPT_PATH=$(dirname "$0")
+$SCRIPT_PATH/run-clang-format.sh
+
+# Can't use tee because it hides the exit code
+if git diff --patch --exit-code > $PATCH_NAME; then
+ echo
+ # print in bold-face green
+ echo "\e[1m\e[32mclang-format and codespell changed nothing."
+ rm -f $PATCH_NAME
+ exit 0;
+else
+ echo
+ # Print in bold-face red
+ echo "\e[1m\e[31mclang-format made at least one change. Run clang-format before pushing!"
+ rm -f $PATCH_NAME
+ exit 1;
+fi
+
diff --git a/scripts/run-clang-format.sh b/scripts/run-clang-format.sh
new file mode 100755
index 0000000000000000000000000000000000000000..cc985fcbfb53db8f3aacfa7281012fc2cdde1631
--- /dev/null
+++ b/scripts/run-clang-format.sh
@@ -0,0 +1,15 @@
+
+#!/bin/sh
+#
+# Author: Jakob Maljaars
+# Email: jakob.maljaars_@_stcorp.nl
+
+set -e
+
+SCRIPT_PATH=$(dirname "$0")
+cd $SCRIPT_PATH
+# Move up to parent folder which contains the source
+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 \{\} +
diff --git a/test/tElementBeamHamaker.cc b/test/tElementBeamHamaker.cc
index 00f1865cca07a663b6c6de301af31ae420fddbd4..bab531d5c305a79248c0942970837753a7762c92 100644
--- a/test/tElementBeamHamaker.cc
+++ b/test/tElementBeamHamaker.cc
@@ -4,11 +4,10 @@
using namespace LOFAR::StationResponse;
-int main(int argc, char** argv)
-{
- ElementResponseModel model(Hamaker);
- double frequency = 132e6; // Mhz
- std::string input_filename(TEST_MEASUREMENTSET);
- std::string output_filename("element-beams-hamaker.fits");
- run(model, frequency, input_filename, output_filename);
+int main(int argc, char** argv) {
+ ElementResponseModel model(Hamaker);
+ double frequency = 132e6; // Mhz
+ std::string input_filename(TEST_MEASUREMENTSET);
+ std::string output_filename("element-beams-hamaker.fits");
+ run(model, frequency, input_filename, output_filename);
}
\ No newline at end of file
diff --git a/test/tElementBeamOSKARDipole.cc b/test/tElementBeamOSKARDipole.cc
index 65a89f934fbb22c18bcb98cc5f0c955acfca05cf..848fbe3fd91d60aad6e3532b04d96cbc759ec0c0 100644
--- a/test/tElementBeamOSKARDipole.cc
+++ b/test/tElementBeamOSKARDipole.cc
@@ -4,11 +4,10 @@
using namespace LOFAR::StationResponse;
-int main(int argc, char** argv)
-{
- ElementResponseModel model(OSKARDipole);
- double frequency = 140e6; // Mhz
- std::string input_filename(TEST_MEASUREMENTSET);
- std::string output_filename("element-beams-oskar-dipole.fits");
- run(model, frequency, input_filename, output_filename);
+int main(int argc, char** argv) {
+ ElementResponseModel model(OSKARDipole);
+ double frequency = 140e6; // Mhz
+ std::string input_filename(TEST_MEASUREMENTSET);
+ std::string output_filename("element-beams-oskar-dipole.fits");
+ run(model, frequency, input_filename, output_filename);
}
\ No newline at end of file
diff --git a/test/tElementBeamOSKARSphericalWave.cc b/test/tElementBeamOSKARSphericalWave.cc
index b6297bb64b9a8a7bd1a08b69673bfdcce099dde6..e7e22f3793fcd362b912d118be01617b9003c3db 100644
--- a/test/tElementBeamOSKARSphericalWave.cc
+++ b/test/tElementBeamOSKARSphericalWave.cc
@@ -4,11 +4,10 @@
using namespace LOFAR::StationResponse;
-int main(int argc, char** argv)
-{
- ElementResponseModel model(OSKARSphericalWave);
- double frequency = 140e6; // Mhz
- std::string input_filename(TEST_MEASUREMENTSET);
- std::string output_filename("element-beams-oskar-sphericalwave.fits");
- run(model, frequency, input_filename, output_filename);
+int main(int argc, char** argv) {
+ ElementResponseModel model(OSKARSphericalWave);
+ double frequency = 140e6; // Mhz
+ std::string input_filename(TEST_MEASUREMENTSET);
+ std::string output_filename("element-beams-oskar-sphericalwave.fits");
+ run(model, frequency, input_filename, output_filename);
}
\ No newline at end of file
diff --git a/test/tStation.cc b/test/tStation.cc
index ac7e5b3973d43427a7478212a573699918057853..94d589eb1969f70b1fee7e59a7f0efb159352e77 100644
--- a/test/tStation.cc
+++ b/test/tStation.cc
@@ -10,62 +10,64 @@
using namespace LOFAR::StationResponse;
-int main()
-{
- std::string name = "CS001LBA";
- vector3r_t position;
-// ElementResponseModel model = ElementResponseModel::Hamaker;
+int main() {
+ std::string name = "CS001LBA";
+ vector3r_t position;
+ // ElementResponseModel model = ElementResponseModel::Hamaker;
-// Station station(name, position, model);
-//
-// double time;
-// double freq;
-// vector3r_t direction = {0.0, 0.0, 1.0};
-// matrix22c_t response;
+ // Station station(name, position, model);
+ //
+ // double time;
+ // double freq;
+ // vector3r_t direction = {0.0, 0.0, 1.0};
+ // matrix22c_t response;
-// auto antenna0 = Element::Ptr(new Element(station.get_element_response(),0));
-// auto antenna1 = Element::Ptr(new Element(station.get_element_response(),1));
-//
-// station.set_antenna(antenna0);
-//
-// std::cout << response[0][0] << std::endl;
-//
-// response = station.response(time, freq, direction);
-//
-// std::cout << response[0][0] << std::endl;
-//
-// auto beam_former = BeamFormer::Ptr(new BeamFormer());
-// beam_former->add_antenna(antenna0);
-// beam_former->add_antenna(antenna1);
-// station.set_antenna(beam_former);
-// response = station.response(time, freq, direction);
-// std::cout << response[0][0] << std::endl;
-//
+ // auto antenna0 = Element::Ptr(new
+ // Element(station.get_element_response(),0)); auto antenna1 =
+ // Element::Ptr(new Element(station.get_element_response(),1));
+ //
+ // station.set_antenna(antenna0);
+ //
+ // std::cout << response[0][0] << std::endl;
+ //
+ // response = station.response(time, freq, direction);
+ //
+ // std::cout << response[0][0] << std::endl;
+ //
+ // auto beam_former = BeamFormer::Ptr(new BeamFormer());
+ // beam_former->add_antenna(antenna0);
+ // beam_former->add_antenna(antenna1);
+ // station.set_antenna(beam_former);
+ // response = station.response(time, freq, direction);
+ // std::cout << response[0][0] << std::endl;
+ //
- casacore::MeasurementSet ms(TEST_MEASUREMENTSET);
+ casacore::MeasurementSet ms(TEST_MEASUREMENTSET);
- double firstTime = casacore::ScalarColumn<double>(ms, "TIME")(0);
- std::cout << firstTime << std::endl;
- double time = firstTime;
- double freq = 55e6;
- matrix22c_t response;
- vector3r_t direction = {0.0, 0.0, 1.0};
+ double firstTime = casacore::ScalarColumn<double>(ms, "TIME")(0);
+ std::cout << firstTime << std::endl;
+ double time = firstTime;
+ double freq = 55e6;
+ matrix22c_t response;
+ vector3r_t direction = {0.0, 0.0, 1.0};
+ // Station::Ptr station = readStation(ms, 0,
+ // ElementResponseModel::OSKARDipole);
+ Station::Ptr station = readStation(ms, 0);
-// Station::Ptr station = readStation(ms, 0, ElementResponseModel::OSKARDipole);
- Station::Ptr station = readStation(ms, 0);
+ auto freq_beamformer = freq;
+ auto station_pointing = direction;
+ auto tile_pointing = direction;
- auto freq_beamformer = freq;
- auto station_pointing = direction;
- auto tile_pointing = direction;
+ for (int i = 0; i < 10; i++) {
+ auto d = direction;
+ d[1] = -0.2 + 0.04 * i;
+ d = normalize(d);
+ response = station->response(time, freq, d, freq_beamformer,
+ station_pointing, tile_pointing);
+ std::cout << response[0][0] << " " << response[0][1] << " "
+ << response[1][0] << " " << response[1][1] << " " << std::endl;
+ }
- for(int i=0; i<10; i++) {
- auto d = direction;
- d[1] = -0.2 + 0.04*i;
- d = normalize(d);
- response = station->response(time, freq, d, freq_beamformer, station_pointing, tile_pointing);
- std::cout << response[0][0] << " " << response[0][1] << " " << response[1][0] << " " << response[1][1] << " " << std::endl;
- }
-
- return 0;
+ return 0;
}
diff --git a/test/tStationBeamHamaker.cc b/test/tStationBeamHamaker.cc
index 3cfc22deaa66bdf5153acda11b029fe205630eea..28f2a52dc3b41ccb07e328559dce74621664538a 100644
--- a/test/tStationBeamHamaker.cc
+++ b/test/tStationBeamHamaker.cc
@@ -4,11 +4,10 @@
using namespace LOFAR::StationResponse;
-int main(int argc, char** argv)
-{
- ElementResponseModel model(Hamaker);
- double frequency = 132e6; // Mhz
- std::string input_filename(TEST_MEASUREMENTSET);
- std::string output_filename("station-beams-hamaker.fits");
- run(model, frequency, input_filename, output_filename);
+int main(int argc, char** argv) {
+ ElementResponseModel model(Hamaker);
+ double frequency = 132e6; // Mhz
+ std::string input_filename(TEST_MEASUREMENTSET);
+ std::string output_filename("station-beams-hamaker.fits");
+ run(model, frequency, input_filename, output_filename);
}
\ No newline at end of file
diff --git a/test/tStationBeamOSKARDipole.cc b/test/tStationBeamOSKARDipole.cc
index bd6c97a773eaad5a3a36f52d6fbe26d39fab19aa..1ea3dca2d0ebb4d98e7943a4ab661262fdb957e6 100644
--- a/test/tStationBeamOSKARDipole.cc
+++ b/test/tStationBeamOSKARDipole.cc
@@ -4,11 +4,10 @@
using namespace LOFAR::StationResponse;
-int main(int argc, char** argv)
-{
- ElementResponseModel model(OSKARDipole);
- double frequency = 140e6; // Mhz
- std::string input_filename(TEST_MEASUREMENTSET);
- std::string output_filename("station-beams-oskar-dipole.fits");
- run(model, frequency, input_filename, output_filename);
+int main(int argc, char** argv) {
+ ElementResponseModel model(OSKARDipole);
+ double frequency = 140e6; // Mhz
+ std::string input_filename(TEST_MEASUREMENTSET);
+ std::string output_filename("station-beams-oskar-dipole.fits");
+ run(model, frequency, input_filename, output_filename);
}
\ No newline at end of file
diff --git a/test/tStationBeamOSKARSphericalWave.cc b/test/tStationBeamOSKARSphericalWave.cc
index d68c9b21f896bafc6570ec4b008cb985d5a391cb..4af0dfabbf67a27297f2bc599199d5f5f006faec 100644
--- a/test/tStationBeamOSKARSphericalWave.cc
+++ b/test/tStationBeamOSKARSphericalWave.cc
@@ -4,11 +4,10 @@
using namespace LOFAR::StationResponse;
-int main(int argc, char** argv)
-{
- ElementResponseModel model(OSKARSphericalWave);
- double frequency = 140e6; // Mhz
- std::string input_filename(TEST_MEASUREMENTSET);
- std::string output_filename("station-beams-oskar-sphericalwave.fits");
- run(model, frequency, input_filename, output_filename);
+int main(int argc, char** argv) {
+ ElementResponseModel model(OSKARSphericalWave);
+ double frequency = 140e6; // Mhz
+ std::string input_filename(TEST_MEASUREMENTSET);
+ std::string output_filename("station-beams-oskar-sphericalwave.fits");
+ run(model, frequency, input_filename, output_filename);
}
\ No newline at end of file