diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index c20ab45cc8d8c57c69f39403a8c84286e505510c..d491e0474dc9e9579634e1d0350ff0f4c19f341f 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -13,6 +13,7 @@ add_library(everybeam SHARED antenna.cc elementresponse.cc beamformer.cc + beamformeridenticalantennas.cc element.cc load.cc coords/itrfconverter.cc @@ -44,6 +45,7 @@ install ( install (FILES antenna.h beamformer.h + beamformeridenticalantennas.h element.h elementresponse.h lofarreadutils.h diff --git a/cpp/beamformer.cc b/cpp/beamformer.cc index c3cc34d87a9835d07ad8aa7ed46ecaab17e70a46..daf3048c00ccee7b6f88067b620a5da9cdc6c653 100644 --- a/cpp/beamformer.cc +++ b/cpp/beamformer.cc @@ -43,7 +43,8 @@ std::vector<std::complex<double>> BeamFormer::ComputeGeometricResponse( std::vector<std::pair<std::complex<double>, std::complex<double>>> BeamFormer::ComputeWeights(const vector3r_t &pointing, double freq) const { // Get geometric response for pointing direction - auto geometric_response = ComputeGeometricResponse(freq, pointing); + std::vector<std::complex<double>> geometric_response = + ComputeGeometricResponse(freq, pointing); // Initialize and fill result double weight_sum[2] = {0.0, 0.0}; @@ -51,7 +52,7 @@ BeamFormer::ComputeWeights(const vector3r_t &pointing, double freq) const { geometric_response.size()); for (std::size_t idx = 0; idx < antennas_.size(); ++idx) { // Compute conjugate of geometric response - auto phasor_conj = std::conj(geometric_response[idx]); + std::complex<double> phasor_conj = std::conj(geometric_response[idx]); // Compute the delays in x/y direction result[idx] = {phasor_conj * (1.0 * antennas_[idx]->enabled_[0]), phasor_conj * (1.0 * antennas_[idx]->enabled_[1])}; @@ -72,16 +73,20 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq, const vector3r_t &direction, const Options &options) const { // Weights based on pointing direction of beam - auto weights = ComputeWeights(options.station0, options.freq0); + std::vector<std::pair<std::complex<double>, std::complex<double>>> weights = + ComputeWeights(options.station0, options.freq0); // Weights based on direction of interest - auto geometric_response = ComputeGeometricResponse(freq, direction); + std::vector<std::complex<double>> geometric_response = + ComputeGeometricResponse(freq, direction); matrix22c_t result = {0}; for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size(); ++antenna_idx) { - auto antenna = antennas_[antenna_idx]; - auto antenna_weight = weights[antenna_idx]; - auto antenna_geometric_reponse = geometric_response[antenna_idx]; + Antenna::Ptr antenna = antennas_[antenna_idx]; + std::pair<std::complex<double>, std::complex<double>> antenna_weight = + weights[antenna_idx]; + std::complex<double> antenna_geometric_reponse = + geometric_response[antenna_idx]; matrix22c_t antenna_response = antenna->Response(time, freq, direction, options); @@ -96,4 +101,30 @@ matrix22c_t BeamFormer::LocalResponse(real_t time, real_t freq, } return result; } + +diag22c_t BeamFormer::LocalArrayFactor(real_t time, real_t freq, + const vector3r_t &direction, + const Options &options) const { + // Weights based on pointing direction of beam + std::vector<std::pair<std::complex<double>, std::complex<double>>> weights = + ComputeWeights(options.station0, options.freq0); + // Weights based on direction of interest + std::vector<std::complex<double>> geometric_response = + ComputeGeometricResponse(freq, direction); + + diag22c_t result = {0}; + for (std::size_t antenna_idx = 0; antenna_idx < antennas_.size(); + ++antenna_idx) { + Antenna::Ptr antenna = antennas_[antenna_idx]; + std::pair<std::complex<double>, std::complex<double>> antenna_weight = + weights[antenna_idx]; + std::complex<double> antenna_geometric_reponse = + geometric_response[antenna_idx]; + + result[0] += antenna_weight.first * antenna_geometric_reponse; + result[1] += antenna_weight.second * antenna_geometric_reponse; + } + return result; +} + } // namespace everybeam diff --git a/cpp/beamformer.h b/cpp/beamformer.h index 86131a9d0c1f2583fb4c0ff9ecdcc89aa172f6f7..6bd52d91bce5e077aab5e90109b0fffc0ba85d56 100644 --- a/cpp/beamformer.h +++ b/cpp/beamformer.h @@ -59,7 +59,7 @@ class BeamFormer : public Antenna { */ void AddAntenna(Antenna::Ptr antenna) { antennas_.push_back(antenna); } - private: + protected: vector3r_t local_phase_reference_position_; // in coordinate system of Antenna @@ -76,9 +76,7 @@ class BeamFormer : public Antenna { // "representation" of Jones matrix virtual diag22c_t LocalArrayFactor(real_t time, real_t freq, const vector3r_t &direction, - const Options &options) const override { - return {1.0, 1.0}; - } + const Options &options) const override; // Compute the geometric response for all the antennas in the BeamFormer based // on the probing frequency and a specified direction (either pointing dir or diff --git a/cpp/beamformeridenticalantennas.cc b/cpp/beamformeridenticalantennas.cc new file mode 100644 index 0000000000000000000000000000000000000000..3aa01a7a8722639d0cfa56ea2c13cfafedde2400 --- /dev/null +++ b/cpp/beamformeridenticalantennas.cc @@ -0,0 +1,30 @@ +#include "beamformeridenticalantennas.h" + +#include "common/constants.h" +#include "common/mathutils.h" + +#include <cmath> + +namespace everybeam { + +matrix22c_t BeamFormerIdenticalAntennas::LocalResponse( + real_t time, real_t freq, const vector3r_t &direction, + const Options &options) const { + matrix22c_t result; + + auto antenna = antennas_[0]; + + matrix22c_t antenna_response = + antenna->Response(time, freq, direction, options); + + diag22c_t array_factor = LocalArrayFactor(time, freq, direction, options); + + result[0][0] = array_factor[0] * antenna_response[0][0]; + result[0][1] = array_factor[0] * antenna_response[0][1]; + result[1][0] = array_factor[1] * antenna_response[1][0]; + result[1][1] = array_factor[1] * antenna_response[1][1]; + + return result; +} + +} // namespace everybeam diff --git a/cpp/beamformeridenticalantennas.h b/cpp/beamformeridenticalantennas.h new file mode 100644 index 0000000000000000000000000000000000000000..31185fd9abb3cb46bd3e46bed97ac63f0ccb9615 --- /dev/null +++ b/cpp/beamformeridenticalantennas.h @@ -0,0 +1,46 @@ +#ifndef EVERYBEAM_BEAMFORMERIDENTICALANTENNAS_H +#define EVERYBEAM_BEAMFORMERIDENTICALANTENNAS_H + +#include "beamformer.h" + +namespace everybeam { +class BeamFormerIdenticalAntennas : public BeamFormer { + public: + /** + * @brief Construct a new BeamFormerIdenticalAntennas object + * + */ + BeamFormerIdenticalAntennas() : BeamFormer() {} + + /** + * @brief Construct a new BeamFormerIdenticalAntennas object given a + * coordinate system. + * + * @param coordinate_system + */ + BeamFormerIdenticalAntennas(const CoordinateSystem &coordinate_system) + : BeamFormer(coordinate_system) {} + + /** + * @brief Construct a new BeamFormer object given a coordinate system and a + * phase reference position + * + * @param coordinate_system + * @param phase_reference_position + */ + BeamFormerIdenticalAntennas(CoordinateSystem coordinate_system, + vector3r_t phase_reference_position) + : BeamFormer(coordinate_system, phase_reference_position) {} + + BeamFormerIdenticalAntennas(vector3r_t phase_reference_position) + : BeamFormer(phase_reference_position) {} + + private: + // Compute the BeamFormer response in certain direction of arrival (ITRF, m) + // and return (Jones) matrix of response + virtual matrix22c_t LocalResponse(real_t time, real_t freq, + const vector3r_t &direction, + const Options &options) const override; +}; +} // namespace everybeam +#endif diff --git a/cpp/lofarreadutils.cc b/cpp/lofarreadutils.cc index c2ee5e74e246078743ca8139ce9bcbcd8eff9a55..2fa2b397d181d28af8657b47becec9b1dd52abd7 100644 --- a/cpp/lofarreadutils.cc +++ b/cpp/lofarreadutils.cc @@ -22,6 +22,7 @@ // $Id$ #include "lofarreadutils.h" +#include "beamformeridenticalantennas.h" #include "common/mathutils.h" #include "common/casautils.h" @@ -146,7 +147,8 @@ vector3r_t TransformToFieldCoordinates( BeamFormer::Ptr MakeTile(unsigned int id, const vector3r_t &position, const TileConfig &tile_config, ElementResponse::Ptr element_response) { - BeamFormer::Ptr tile = BeamFormer::Ptr(new BeamFormer(position)); + BeamFormer::Ptr tile = + BeamFormer::Ptr(new BeamFormerIdenticalAntennas(position)); for (unsigned int id = 0; id < tile_config.size(); id++) { vector3r_t antenna_position = tile_config[id]; @@ -175,7 +177,8 @@ BeamFormer::Ptr ReadAntennaField(const Table &table, unsigned int id, // 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)); + BeamFormer::Ptr beam_former( + new BeamFormerIdenticalAntennas(coordinate_system)); ROScalarColumn<String> c_name(table, "NAME"); ROArrayQuantColumn<Double> c_offset(table, "ELEMENT_OFFSET", "m");