diff --git a/CMakeLists.txt b/CMakeLists.txt
index aefd9597eff049bef8f6e3793edbef91c2302403..adf9cbce91caa7c2dab9f8428338d4e9dfe0e682 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -7,7 +7,7 @@ cmake_minimum_required(VERSION 3.7)
#------------------------------------------------------------------------------
# Set version name and project number
-set(EVERYBEAM_VERSION 0.2.3)
+set(EVERYBEAM_VERSION 0.3.0)
if(EVERYBEAM_VERSION MATCHES "^([0-9]+)\\.([0-9]+)\\.([0-9]+)")
set(EVERYBEAM_VERSION_MAJOR "${CMAKE_MATCH_1}")
set(EVERYBEAM_VERSION_MINOR "${CMAKE_MATCH_2}")
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index ae1bfe3f41f527c147aeeb5157ba0942d7612030..52b87e80e138f87aca7df9afdb3d724360825a77 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -46,6 +46,7 @@ add_library(
coords/itrfdirection.cc
msreadutils.cc
station.cc
+ phasedarrayresponse.cc
telescope/lofar.cc
telescope/dish.cc
telescope/mwa.cc
@@ -102,6 +103,7 @@ install(
station.h
load.h
options.h
+ phasedarrayresponse.h
DESTINATION "include/${CMAKE_PROJECT_NAME}")
# TODO: not sure whether this is needed at all
diff --git a/cpp/aterms/atermconfig.cc b/cpp/aterms/atermconfig.cc
index 0021456d1aa011d334293f20d0950d625cd2c86d..c2bd0642c5cfebb43b04911c6bf08841958256bd 100644
--- a/cpp/aterms/atermconfig.cc
+++ b/cpp/aterms/atermconfig.cc
@@ -289,19 +289,6 @@ std::unique_ptr<ATermBeam> ATermConfig::GetATermBeam(
new EveryBeamATerm(ms, coordinate_system, options));
}
-std::unique_ptr<ATermBeam> ATermConfig::GetATermBeam(
- const casacore::MeasurementSet& ms,
- const CoordinateSystem& coordinate_system, const ATermSettings& settings,
- bool frequency_interpolation, bool use_differential_beam,
- bool use_channel_frequency, const std::string& element_response_model,
- const std::string& beam_mode) {
- everybeam::Options options = ConvertToEBOptions(
- ms, settings, frequency_interpolation, use_differential_beam,
- use_channel_frequency, element_response_model, beam_mode);
- return std::unique_ptr<ATermBeam>(
- new EveryBeamATerm(ms, coordinate_system, options));
-}
-
everybeam::Options ATermConfig::ConvertToEBOptions(
const casacore::MeasurementSet& ms, const ATermSettings& settings,
bool frequency_interpolation, const std::string& beam_normalisation_mode,
@@ -330,21 +317,5 @@ everybeam::Options ATermConfig::ConvertToEBOptions(
options.beam_normalisation_mode = beam_normalisation_mode_enum;
return options;
}
-
-everybeam::Options ATermConfig::ConvertToEBOptions(
- const casacore::MeasurementSet& ms, const ATermSettings& settings,
- bool frequency_interpolation, bool use_differential_beam,
- bool use_channel_frequency, const std::string& element_response_model,
- const std::string& beam_mode) {
- const std::string beam_normalisation_mode_none_str = "none";
- everybeam::Options options = ConvertToEBOptions(
- ms, settings, frequency_interpolation, beam_normalisation_mode_none_str,
- use_channel_frequency, element_response_model, beam_mode);
- options.beam_normalisation_mode = use_differential_beam
- ? BeamNormalisationMode::kFull
- : BeamNormalisationMode::kPreApplied;
- return options;
-}
-
} // namespace aterms
} // namespace everybeam
diff --git a/cpp/aterms/atermconfig.h b/cpp/aterms/atermconfig.h
index 21391e483bff29ffff61738bd87bfcad162d0764..7bd62737524383d41075eb87fcf42077b1515e2f 100644
--- a/cpp/aterms/atermconfig.h
+++ b/cpp/aterms/atermconfig.h
@@ -55,30 +55,6 @@ class ATermConfig final : public ATermBase {
return avgTime;
}
- /**
- * @brief Static method for constructing an (EveryBeam)ATerm
- *
- * @param ms Measurement set
- * @param coordinate_system struct with image settings
- * @param settings aterm specific settings
- * @param frequency_interpolation Interpolate between frequencies? Relevant
- * for MWA only.
- * @param use_differential_beam Apply differential beam
- * @param use_channel_frequency Use channel frequency
- * @param element_response_model Element response model
- * @return std::unique_ptr<ATermBeam>
- */
- [[deprecated(
- "Use other overload of GetATermBeam() instead. Function will be removed "
- "in version "
- "0.3.0.")]] static std::unique_ptr<ATermBeam>
- GetATermBeam(const casacore::MeasurementSet& ms,
- const coords::CoordinateSystem& coordinate_system,
- const ATermSettings& settings, bool frequency_interpolation,
- bool use_differential_beam, bool use_channel_frequency,
- const std::string& element_response_model,
- const std::string& beam_mode = "full");
-
/**
* @brief Static method for constructing an (EveryBeam)ATerm
*
@@ -119,17 +95,6 @@ class ATermConfig final : public ATermBase {
bool use_channel_frequency, const std::string& element_response_model,
const std::string& beam_mode = "full");
- [[deprecated(
- "Use other overload of GetATermBeam() instead. Function will be removed "
- "in version "
- "0.3.0.")]] static everybeam::Options
- ConvertToEBOptions(const casacore::MeasurementSet& ms,
- const ATermSettings& settings,
- bool frequency_interpolation, bool use_differential_beam,
- bool use_channel_frequency,
- const std::string& element_response_model,
- const std::string& beam_mode = "full");
-
private:
const size_t n_antennas_;
const coords::CoordinateSystem coordinate_system_;
diff --git a/cpp/griddedresponse/aartfaacgrid.cc b/cpp/griddedresponse/aartfaacgrid.cc
index b16989582a4a9a18c5709ec1943fa27e54b7e18a..82f2985a367ef49b2cd3cddc513363320a9aa05b 100644
--- a/cpp/griddedresponse/aartfaacgrid.cc
+++ b/cpp/griddedresponse/aartfaacgrid.cc
@@ -3,6 +3,8 @@
#include "aartfaacgrid.h"
+#include <aocommon/uvector.h>
+
namespace everybeam {
namespace griddedresponse {
@@ -27,7 +29,7 @@ void AartfaacGrid::MakeIntegratedSnapshot(
ResponseAllStations(beam_mode, buffer_undersampled.data(), time, frequency,
field_id);
- // For Aartfaac, we can simply weight a (time) snapshot with the accumulated
+ // For Aartfaac, simply weight a (time) snapshot with the accumulated
// baseline weights
const size_t n_baselines = n_stations * (n_stations + 1) / 2;
double snapshot_weight = 0.;
diff --git a/cpp/griddedresponse/dishgrid.cc b/cpp/griddedresponse/dishgrid.cc
index af502e513780ac863450874897e74929a19db244..c78e5dbcae8295b67a3632ac9da2104febf379c9 100644
--- a/cpp/griddedresponse/dishgrid.cc
+++ b/cpp/griddedresponse/dishgrid.cc
@@ -24,19 +24,18 @@ void DishGrid::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
const telescope::Dish& dish_telescope =
static_cast<const telescope::Dish&>(*telescope_);
- const double pdir_ra =
- dish_telescope.ms_properties_.field_pointing[field_id].first;
- const double pdir_dec =
- dish_telescope.ms_properties_.field_pointing[field_id].second;
+ double pdir_ra;
+ double pdir_dec;
+ std::tie(pdir_ra, pdir_dec) = dish_telescope.GetFieldPointing()[field_id];
const double max_radius_arc_min =
- dish_telescope.coefficients_->MaxRadiusInArcMin();
+ dish_telescope.GetDishCoefficients()->MaxRadiusInArcMin();
const double reference_frequency =
- dish_telescope.coefficients_->ReferenceFrequency();
+ dish_telescope.GetDishCoefficients()->ReferenceFrequency();
circularsymmetric::VoltagePattern vp(
- dish_telescope.coefficients_->GetFrequencies(frequency),
+ dish_telescope.GetDishCoefficients()->GetFrequencies(frequency),
max_radius_arc_min, reference_frequency);
const aocommon::UVector<double> coefs_vec =
- dish_telescope.coefficients_->GetCoefficients(frequency);
+ dish_telescope.GetDishCoefficients()->GetCoefficients(frequency);
vp.EvaluatePolynomial(coefs_vec, false);
vp.Render(buffer, width_, height_, dl_, dm_, ra_, dec_, pdir_ra, pdir_dec,
phase_centre_dl_, phase_centre_dm_, frequency);
diff --git a/cpp/griddedresponse/griddedresponse.cc b/cpp/griddedresponse/griddedresponse.cc
index e564a12e45ab236002cc6e33d3046e012a606c47..1b0d3b81b9f5135dc00660573e1daa228867b10f 100644
--- a/cpp/griddedresponse/griddedresponse.cc
+++ b/cpp/griddedresponse/griddedresponse.cc
@@ -6,8 +6,10 @@
#include "../telescope/telescope.h"
+#include <aocommon/uvector.h>
#include <aocommon/matrix2x2.h>
#include <aocommon/matrix4x4.h>
+#include <aocommon/hmatrix4x4.h>
#include <vector>
#include <complex>
#include <stdexcept>
@@ -16,9 +18,9 @@ using aocommon::HMC4x4;
using aocommon::MC2x2;
using aocommon::MC4x4;
using aocommon::UVector;
-using everybeam::common::FFTResampler;
-using everybeam::griddedresponse::GriddedResponse;
+namespace everybeam {
+namespace griddedresponse {
void GriddedResponse::IntegratedResponse(
BeamMode beam_mode, float* buffer, double time, double frequency,
size_t field_id, size_t undersampling_factor,
@@ -163,7 +165,7 @@ void GriddedResponse::DoFFTResampling(
float* buffer, int width_in, int height_in, int width_out, int height_out,
const std::vector<aocommon::HMC4x4>& matrices) {
// (FFT) resampling, run multi-threaded?
- FFTResampler resampler(width_in, height_in, width_out, height_out);
+ common::FFTResampler resampler(width_in, height_in, width_out, height_out);
resampler.SetWindowFunction(aocommon::WindowFunction::RaisedHann, true);
UVector<float> lowres_input(width_in * height_in);
// Loop over the "double" representation of the HMC4x4 Hermitian matrix
@@ -176,3 +178,5 @@ void GriddedResponse::DoFFTResampling(
buffer + p * width_out * height_out);
}
}
+} // namespace griddedresponse
+} // namespace everybeam
\ No newline at end of file
diff --git a/cpp/griddedresponse/griddedresponse.h b/cpp/griddedresponse/griddedresponse.h
index 7f99c1cc6cc584433d6af1720c8f98b07951c737..9a422a28ac17e4413777c85987ffc036c70b5876 100644
--- a/cpp/griddedresponse/griddedresponse.h
+++ b/cpp/griddedresponse/griddedresponse.h
@@ -11,14 +11,12 @@
#include "./../coords/itrfconverter.h"
#include "../beammode.h"
+#include <aocommon/hmatrix4x4.h>
+
#include <memory>
#include <vector>
#include <thread>
-#include <aocommon/lane.h>
-#include <aocommon/uvector.h>
-#include <aocommon/hmatrix4x4.h>
-
#include <casacore/measures/Measures/MDirection.h>
namespace everybeam {
@@ -37,37 +35,6 @@ class GriddedResponse {
public:
virtual ~GriddedResponse() {}
- /**
- * @brief See Response()
- */
- [[deprecated(
- "Use Response() instead. Function will be removed in version "
- "0.3.0.")]] void
- CalculateStation(std::complex<float>* buffer, double time, double freq,
- size_t station_idx, size_t field_id) {
- Response(BeamMode::kFull, buffer, time, freq, station_idx, field_id);
- };
-
- /**
- * @brief Compute the full beam response for a single station. Result is
- * stored in the output buffer, which should accommodate a Jones matrix (4
- * complex valued floats) per pixel.
- *
- * @param buffer Output buffer, compute and set size with
- * GriddedResponse::GetStationBufferSize(1)
- * @param station_idx Station index, must be smaller than number of stations
- * in the Telescope
- * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
- * @param frequency Frequency (Hz)
- */
- [[deprecated(
- "Use Response() instead. Function will be removed in version "
- "0.3.0.")]] void
- FullResponse(std::complex<float>* buffer, double time, double freq,
- size_t station_idx, size_t field_id) {
- Response(BeamMode::kFull, buffer, time, freq, station_idx, field_id);
- }
-
/**
* @brief Compute the beam for a single station, given a prescribed beam mode.
* Result is stored in the output buffer, which should accommodate a Jones
@@ -86,37 +53,6 @@ class GriddedResponse {
double time, double freq, size_t station_idx,
size_t field_id) = 0;
- /**
- * @brief See ResponseAllStations()
- */
- [[deprecated(
- "Use ResponseAllStations() instead. Function will be removed in "
- "version "
- "0.3.0.")]] void
- CalculateAllStations(std::complex<float>* buffer, double time,
- double frequency, size_t field_id) {
- ResponseAllStations(BeamMode::kFull, buffer, time, frequency, field_id);
- };
-
- /**
- * @brief Compute the full beam response for all stations in a Telescope.
- * Result is stored in the output buffer, which should accommodate a Jones
- * matrix (4 complex valued floats) per pixel for each station.
- *
- * @param buffer Output buffer, compute and set size with
- * GriddedResponse::GetStationBufferSize()
- * @param time Time, modified Julian date, UTC, in seconds (MJD(UTC), s).
- * @param frequency Frequency (Hz)
- */
- [[deprecated(
- "Use ResponseAllStations() instead. Function will be removed in "
- "version "
- "0.3.0.")]] void
- FullResponseAllStations(std::complex<float>* buffer, double time,
- double frequency, size_t field_id) {
- ResponseAllStations(BeamMode::kFull, buffer, time, frequency, field_id);
- }
-
/**
* @brief Compute the array factor for all stations in a Telescope.
* Result is stored in the output buffer, which should accommodate a Jones
@@ -169,19 +105,6 @@ class GriddedResponse {
size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights);
- [[deprecated(
- "Overload will be kept alive to allow WSClean migrating to float "
- "implementation")]] void
- IntegratedResponse(BeamMode beam_mode, double* buffer, double time,
- double frequency, size_t field_id,
- size_t undersampling_factor,
- const std::vector<double>& baseline_weights) {
- std::vector<float> buffer_float(GetIntegratedBufferSize());
- IntegratedResponse(beam_mode, buffer_float.data(), time, frequency,
- field_id, undersampling_factor, baseline_weights);
- std::copy_n(buffer_float.begin(), GetIntegratedBufferSize(), buffer);
- }
-
/**
* @brief Calculate integrated beam over multiple time steps.
* This function makes use of @ref MakeIntegratedSnapshot(). Subclasses
@@ -206,19 +129,6 @@ class GriddedResponse {
size_t undersampling_factor,
const std::vector<double>& baseline_weights);
- [[deprecated(
- "Overload will be kept alive to allow WSClean migrating to float "
- "implementation")]] void
- IntegratedResponse(BeamMode beam_mode, double* buffer,
- const std::vector<double>& time_array, double frequency,
- size_t field_id, size_t undersampling_factor,
- const std::vector<double>& baseline_weights) {
- std::vector<float> buffer_float(GetIntegratedBufferSize());
- IntegratedResponse(beam_mode, buffer_float.data(), time_array, frequency,
- field_id, undersampling_factor, baseline_weights);
- std::copy_n(buffer_float.begin(), GetIntegratedBufferSize(), buffer);
- }
-
std::size_t GetStationBufferSize(std::size_t nstations) const {
return nstations * width_ * height_ * 4u;
}
diff --git a/cpp/griddedresponse/mwagrid.cc b/cpp/griddedresponse/mwagrid.cc
index 08ef6852d3a591b5fd2e439b4ca187a9c460893c..a9c3bc6e7d446c8f642af70bbb978532fa352a5b 100644
--- a/cpp/griddedresponse/mwagrid.cc
+++ b/cpp/griddedresponse/mwagrid.cc
@@ -23,8 +23,7 @@ void MWAGrid::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
static_cast<const telescope::MWA&>(*telescope_);
casacore::MEpoch time_epoch(casacore::Quantity(time, "s"));
- casacore::MeasFrame frame(mwatelescope.ms_properties_.array_position,
- time_epoch);
+ casacore::MeasFrame frame(mwatelescope.GetArrayPosition(), time_epoch);
const casacore::MDirection::Ref hadec_ref(casacore::MDirection::HADEC, frame);
const casacore::MDirection::Ref azelgeo_ref(casacore::MDirection::AZELGEO,
@@ -33,12 +32,12 @@ void MWAGrid::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
casacore::MDirection::Convert j2000_to_hadecref(j2000_ref, hadec_ref),
j2000_to_azelgeoref(j2000_ref, azelgeo_ref);
casacore::MPosition wgs = casacore::MPosition::Convert(
- mwatelescope.ms_properties_.array_position, casacore::MPosition::WGS84)();
- double arr_latitude = wgs.getValue().getLat();
+ mwatelescope.GetArrayPosition(), casacore::MPosition::WGS84)();
+ const double arr_latitude = wgs.getValue().getLat();
if (!tile_beam_) {
tile_beam_.reset(
- new TileBeam2016(mwatelescope.ms_properties_.delays,
+ new TileBeam2016(mwatelescope.GetDelays().data(),
mwatelescope.GetOptions().frequency_interpolation,
mwatelescope.GetOptions().coeff_path));
}
diff --git a/cpp/griddedresponse/phasedarraygrid.cc b/cpp/griddedresponse/phasedarraygrid.cc
index 9044b732e48b5862fbbdac4ea071053c60c915dc..8fd0a1505fc89a6c0a24ccba6ed1f6416e110fa2 100644
--- a/cpp/griddedresponse/phasedarraygrid.cc
+++ b/cpp/griddedresponse/phasedarraygrid.cc
@@ -4,31 +4,22 @@
#include "phasedarraygrid.h"
#include "../telescope/phasedarray.h"
+#include <aocommon/lane.h>
#include <aocommon/imagecoordinates.h>
#include <aocommon/threadpool.h>
#include <cmath>
#include <iostream>
+namespace everybeam {
-using everybeam::griddedresponse::PhasedArrayGrid;
+using telescope::PhasedArray;
+
+namespace griddedresponse {
PhasedArrayGrid::PhasedArrayGrid(
const telescope::Telescope* telescope_ptr,
const coords::CoordinateSystem& coordinate_system)
: GriddedResponse(telescope_ptr, coordinate_system),
- beam_normalisation_mode_(
- telescope_->GetOptions().beam_normalisation_mode) {
- // Extract PhasedArrayPoint specific options from ms_properties_ and
- // telescope::Options
- const telescope::PhasedArray& phasedarray =
- dynamic_cast<const telescope::PhasedArray&>(*telescope_ptr);
- delay_dir_ = phasedarray.GetMSProperties().delay_dir;
- tile_beam_dir_ = phasedarray.GetMSProperties().tile_beam_dir;
- preapplied_beam_dir_ = phasedarray.GetMSProperties().preapplied_beam_dir;
- preapplied_correction_mode_ =
- phasedarray.GetMSProperties().preapplied_correction_mode;
- subband_frequency_ = phasedarray.GetMSProperties().subband_freq;
- use_channel_frequency_ = phasedarray.GetOptions().use_channel_frequency;
-
+ PhasedArrayResponse(static_cast<const PhasedArray*>(telescope_ptr)) {
// Compute and set number of threads
const size_t ncpus = aocommon::ThreadPool::NCPUs();
const size_t nthreads = std::min(ncpus, telescope_->GetNrStations());
@@ -66,7 +57,7 @@ void PhasedArrayGrid::ResponseAllStations(BeamMode beam_mode,
std::complex<float>* buffer,
double time, double frequency,
size_t) {
- const telescope::PhasedArray& phasedarraytelescope =
+ const telescope::PhasedArray& phased_array =
static_cast<const telescope::PhasedArray&>(*telescope_);
aocommon::Lane<Job> lane(threads_.size());
lane_ = &lane;
@@ -74,8 +65,8 @@ void PhasedArrayGrid::ResponseAllStations(BeamMode beam_mode,
SetITRFVectors(time);
bool apply_normalisation;
- inverse_central_gain_.resize(phasedarraytelescope.GetNrStations());
- for (size_t i = 0; i != phasedarraytelescope.GetNrStations(); ++i) {
+ inverse_central_gain_.resize(phased_array.GetNrStations());
+ for (size_t i = 0; i != phased_array.GetNrStations(); ++i) {
apply_normalisation = CalculateBeamNormalisation(
beam_mode, time, frequency, i, inverse_central_gain_[i]);
}
@@ -87,8 +78,8 @@ void PhasedArrayGrid::ResponseAllStations(BeamMode beam_mode,
}
for (size_t y = 0; y != height_; ++y) {
- for (size_t antenna_idx = 0;
- antenna_idx != phasedarraytelescope.GetNrStations(); ++antenna_idx) {
+ for (size_t antenna_idx = 0; antenna_idx != phased_array.GetNrStations();
+ ++antenna_idx) {
lane.write(Job(y, antenna_idx, antenna_idx));
}
}
@@ -126,72 +117,10 @@ void PhasedArrayGrid::SetITRFVectors(double time) {
diff_beam_centre_);
}
-bool PhasedArrayGrid::CalculateBeamNormalisation(
- BeamMode beam_mode, double time, double frequency, size_t station_idx,
- aocommon::MC2x2F& inverse_gain) const {
- const telescope::PhasedArray& phasedarraytelescope =
- static_cast<const telescope::PhasedArray&>(*telescope_);
- if (beam_normalisation_mode_ == BeamNormalisationMode::kNone) {
- return false;
- }
-
- const double sb_freq =
- use_channel_frequency_ ? frequency : subband_frequency_;
-
- // if the normalisation mode is kPreApplied, but no beam correction was pre
- // applied then there is nothing to do
- if (beam_normalisation_mode_ == BeamNormalisationMode::kPreApplied &&
- preapplied_correction_mode_ == BeamMode::kNone) {
- return false;
- }
-
- // If the normalisation mode is kPreApplied, or kPreAppliedOrFull and the
- // fallback to Full is not needed then the response for the diff_beam_centre_
- // with preapplied_correction_mode_ needs to be computed
- if (beam_normalisation_mode_ == BeamNormalisationMode::kPreApplied ||
- (beam_normalisation_mode_ == BeamNormalisationMode::kPreAppliedOrFull &&
- preapplied_correction_mode_ != BeamMode::kNone)) {
- inverse_gain = aocommon::MC2x2F(
- phasedarraytelescope.GetStation(station_idx)
- ->Response(preapplied_correction_mode_, time, frequency,
- diff_beam_centre_, sb_freq, station0_, tile0_)
- .Data());
- } else {
- // in all other cases the response for the reference direction with
- // beam_mode is needed
- inverse_gain = aocommon::MC2x2F(
- phasedarraytelescope.GetStation(station_idx)
- ->Response(beam_mode, time, frequency, diff_beam_centre_, sb_freq,
- station0_, tile0_)
- .Data());
- }
-
- switch (beam_normalisation_mode_) {
- case BeamNormalisationMode::kFull:
- case BeamNormalisationMode::kPreApplied:
- case BeamNormalisationMode::kPreAppliedOrFull:
- if (!inverse_gain.Invert()) {
- inverse_gain = aocommon::MC2x2F::Zero();
- }
- break;
- case BeamNormalisationMode::kAmplitude: {
- const float amplitude_inv = 1.0 / std::sqrt(0.5 * Norm(inverse_gain));
- inverse_gain[0] = std::isfinite(amplitude_inv) ? amplitude_inv : 0.0;
- inverse_gain[1] = 0.0;
- inverse_gain[2] = 0.0;
- inverse_gain[3] = std::isfinite(amplitude_inv) ? amplitude_inv : 0.0;
- break;
- }
- default:
- throw std::runtime_error("Invalid beam normalisation mode here");
- }
- return true;
-}
-
void PhasedArrayGrid::CalcThread(BeamMode beam_mode, bool apply_normalisation,
std::complex<float>* buffer, double time,
double frequency) {
- const telescope::PhasedArray& phasedarraytelescope =
+ const telescope::PhasedArray& phased_array =
static_cast<const telescope::PhasedArray&>(*telescope_);
const size_t values_per_ant = width_ * height_ * 4;
const double sb_freq =
@@ -216,7 +145,7 @@ void PhasedArrayGrid::CalcThread(BeamMode beam_mode, bool apply_normalisation,
base_buffer + job.buffer_offset * values_per_ant;
const aocommon::MC2x2F gain_matrix = aocommon::MC2x2F(
- phasedarraytelescope.GetStation(job.antenna_idx)
+ phased_array.GetStation(job.antenna_idx)
->Response(beam_mode, time, frequency, itrf_direction, sb_freq,
station0_, tile0_)
.Data());
@@ -231,3 +160,5 @@ void PhasedArrayGrid::CalcThread(BeamMode beam_mode, bool apply_normalisation,
}
}
}
+} // namespace griddedresponse
+} // namespace everybeam
\ No newline at end of file
diff --git a/cpp/griddedresponse/phasedarraygrid.h b/cpp/griddedresponse/phasedarraygrid.h
index 5ceea6fb1e1f77b355661f662191a06e9debe319..bff2507fed07406645d89be04b4a2806b9e034df 100644
--- a/cpp/griddedresponse/phasedarraygrid.h
+++ b/cpp/griddedresponse/phasedarraygrid.h
@@ -9,10 +9,16 @@
#include "griddedresponse.h"
#include "../beamnormalisationmode.h"
+#include "../phasedarrayresponse.h"
+
+namespace aocommon {
+template <typename Tp>
+class Lane;
+}
namespace everybeam {
namespace griddedresponse {
-class PhasedArrayGrid : public GriddedResponse {
+class PhasedArrayGrid : public GriddedResponse, protected PhasedArrayResponse {
public:
PhasedArrayGrid(const telescope::Telescope* telescope_ptr,
const coords::CoordinateSystem& coordinate_system);
@@ -25,17 +31,10 @@ class PhasedArrayGrid : public GriddedResponse {
double time, double frequency,
size_t field_id) override;
- protected:
- casacore::MDirection delay_dir_, tile_beam_dir_, preapplied_beam_dir_;
- vector3r_t station0_, tile0_, l_vector_itrf_, m_vector_itrf_, n_vector_itrf_,
- diff_beam_centre_;
-
- CorrectionMode preapplied_correction_mode_;
- BeamNormalisationMode beam_normalisation_mode_;
- bool use_channel_frequency_;
- double subband_frequency_;
-
private:
+ vector3r_t l_vector_itrf_;
+ vector3r_t m_vector_itrf_;
+ vector3r_t n_vector_itrf_;
std::vector<aocommon::MC2x2F> inverse_central_gain_;
std::vector<std::thread> threads_;
@@ -48,10 +47,6 @@ class PhasedArrayGrid : public GriddedResponse {
};
aocommon::Lane<Job>* lane_;
- bool CalculateBeamNormalisation(BeamMode beam_mode, double time,
- double frequency, size_t station_idx,
- aocommon::MC2x2F& inverse_gain) const;
-
/**
* @brief Method for computing the ITRF-vectors.
* NOTE: method not thread-safe due to casacore dependencies.
diff --git a/cpp/phasedarrayresponse.cc b/cpp/phasedarrayresponse.cc
new file mode 100644
index 0000000000000000000000000000000000000000..dc67cad88ce7fe65c38bcb523ed4cd3fae171aec
--- /dev/null
+++ b/cpp/phasedarrayresponse.cc
@@ -0,0 +1,91 @@
+// Copyright (C) 2020 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#include "telescope/phasedarray.h"
+#include "phasedarrayresponse.h"
+
+#include <cassert>
+
+namespace everybeam {
+
+PhasedArrayResponse::PhasedArrayResponse(
+ const telescope::PhasedArray* phased_array)
+ : delay_dir_(phased_array->GetMSProperties().delay_dir),
+ tile_beam_dir_(phased_array->GetMSProperties().tile_beam_dir),
+ preapplied_beam_dir_(phased_array->GetMSProperties().preapplied_beam_dir),
+ preapplied_correction_mode_(
+ phased_array->GetMSProperties().preapplied_correction_mode),
+ beam_normalisation_mode_(
+ phased_array->GetOptions().beam_normalisation_mode),
+ use_channel_frequency_(phased_array->GetOptions().use_channel_frequency),
+ subband_frequency_(phased_array->GetMSProperties().subband_freq),
+ phased_array_(phased_array) {
+ assert(phased_array_);
+}
+
+bool PhasedArrayResponse::CalculateBeamNormalisation(
+ BeamMode beam_mode, double time, double frequency, size_t station_index,
+ aocommon::MC2x2F& inverse_gain) const {
+ const telescope::PhasedArray& phased_array =
+ static_cast<const telescope::PhasedArray&>(*phased_array_);
+ if (beam_normalisation_mode_ == BeamNormalisationMode::kNone) {
+ return false;
+ }
+
+ const double subband_frequency =
+ use_channel_frequency_ ? frequency : subband_frequency_;
+
+ // if the normalisation mode is kPreApplied, but no beam correction was pre
+ // applied then there is nothing to do
+ if (beam_normalisation_mode_ == BeamNormalisationMode::kPreApplied &&
+ preapplied_correction_mode_ == BeamMode::kNone) {
+ return false;
+ }
+
+ // If the normalisation mode is kPreApplied, or kPreAppliedOrFull and the
+ // fallback to Full is not needed then the response for the diff_beam_centre_
+ // with preapplied_correction_mode_ needs to be computed
+ if (beam_normalisation_mode_ == BeamNormalisationMode::kPreApplied ||
+ (beam_normalisation_mode_ == BeamNormalisationMode::kPreAppliedOrFull &&
+ preapplied_correction_mode_ != BeamMode::kNone)) {
+ inverse_gain = aocommon::MC2x2F(
+ phased_array.GetStation(station_index)
+ ->Response(preapplied_correction_mode_, time, frequency,
+ diff_beam_centre_, subband_frequency, station0_, tile0_)
+ .Data());
+ } else {
+ // in all other cases the response for the reference direction with
+ // beam_mode is needed
+ inverse_gain = aocommon::MC2x2F(
+ phased_array.GetStation(station_index)
+ ->Response(beam_mode, time, frequency, diff_beam_centre_,
+ subband_frequency, station0_, tile0_)
+ .Data());
+ }
+
+ switch (beam_normalisation_mode_) {
+ case BeamNormalisationMode::kFull:
+ case BeamNormalisationMode::kPreApplied:
+ case BeamNormalisationMode::kPreAppliedOrFull:
+ if (!inverse_gain.Invert()) {
+ inverse_gain = aocommon::MC2x2F::Zero();
+ }
+ break;
+ case BeamNormalisationMode::kAmplitude: {
+ const float norm_inverse_gain = Norm(inverse_gain);
+ const float amplitude_inv =
+ (norm_inverse_gain == 0.0) ? 0.0
+ : 1.0 / std::sqrt(0.5 * norm_inverse_gain);
+ inverse_gain[0] = amplitude_inv;
+ inverse_gain[1] = 0.0;
+ inverse_gain[2] = 0.0;
+ inverse_gain[3] = amplitude_inv;
+ break;
+ }
+ case BeamNormalisationMode::kNone:
+ throw std::runtime_error("Invalid beam normalisation mode here");
+ }
+ return true;
+}
+
+} // namespace everybeam
\ No newline at end of file
diff --git a/cpp/phasedarrayresponse.h b/cpp/phasedarrayresponse.h
new file mode 100644
index 0000000000000000000000000000000000000000..34899c73b08c6038d7fd3c3808cfc6f1a0cfa11f
--- /dev/null
+++ b/cpp/phasedarrayresponse.h
@@ -0,0 +1,72 @@
+// phasedarrayresponse.h: Common functionality for PhasedArrayPoint and
+// PhasedArrayGrid.
+//
+// Copyright (C) 2021 ASTRON (Netherlands Institute for Radio Astronomy)
+// SPDX-License-Identifier: GPL-3.0-or-later
+
+#ifndef EVERYBEAM_PHASEDARRAYRESPONSE_PHASEDARRAYRESPONSE_H_
+#define EVERYBEAM_PHASEDARRAYRESPONSE_PHASEDARRAYRESPONSE_H_
+
+#include "common/types.h"
+#include "correctionmode.h"
+#include "beammode.h"
+#include "beamnormalisationmode.h"
+
+#include <aocommon/matrix2x2.h>
+#include <casacore/measures/Measures/MDirection.h>
+
+namespace everybeam {
+namespace telescope {
+
+/**
+ * @brief Class containing shared functionality and member variables
+ * for phased array beamformer responses.
+ *
+ * TODO: this class might migrate to a better suited namespace in the future.
+ */
+class PhasedArray;
+} // namespace telescope
+
+class PhasedArrayResponse {
+ public:
+ ~PhasedArrayResponse() = default;
+
+ /**
+ * @brief Determine whether beam normalisation should be applied,
+ * and calculate inverse central gain if
+ *
+ * @param beam_mode Requested beam mode (kNone, kFull, kArrayFactor, kElement)
+ * @param time Time (J2000,s)
+ * @param frequency Frequency (Hz)
+ * @param station_index Station index
+ * @param inverse_gain Inverse central gain with which to normalize beam.
+ * @return Is beam normalisation required?
+ */
+ bool CalculateBeamNormalisation(BeamMode beam_mode, double time,
+ double frequency, size_t station_index,
+ aocommon::MC2x2F& inverse_gain) const;
+
+ protected:
+ PhasedArrayResponse(const telescope::PhasedArray* phasedarray);
+
+ vector3r_t station0_;
+ vector3r_t tile0_;
+ vector3r_t diff_beam_centre_;
+
+ const casacore::MDirection delay_dir_;
+ const casacore::MDirection tile_beam_dir_;
+ const casacore::MDirection preapplied_beam_dir_;
+ const CorrectionMode preapplied_correction_mode_;
+
+ const BeamNormalisationMode beam_normalisation_mode_;
+ // non-const, since OSKAR overwrites these member variables
+ // to correct default
+ bool use_channel_frequency_;
+ double subband_frequency_;
+
+ private:
+ const telescope::PhasedArray* phased_array_;
+};
+} // namespace everybeam
+
+#endif // EVERYBEAM_PHASEDARRAYRESPONSE_PHASEDARRAYRESPONSE_H_
\ No newline at end of file
diff --git a/cpp/pointresponse/dishpoint.cc b/cpp/pointresponse/dishpoint.cc
index 5b077806560f4e63d0967a6e2139e104bc4ca88b..29f7cad24849e06bbe9078d90f25b09245594890 100644
--- a/cpp/pointresponse/dishpoint.cc
+++ b/cpp/pointresponse/dishpoint.cc
@@ -16,20 +16,20 @@ void DishPoint::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
const telescope::Dish& dish_telescope =
static_cast<const telescope::Dish&>(*telescope_);
- const double pdir_ra =
- dish_telescope.ms_properties_.field_pointing[field_id].first;
- const double pdir_dec =
- dish_telescope.ms_properties_.field_pointing[field_id].second;
+ double pdir_ra;
+ double pdir_dec;
+ std::tie(pdir_ra, pdir_dec) = dish_telescope.GetFieldPointing()[field_id];
const double max_radius_arc_min =
- dish_telescope.coefficients_->MaxRadiusInArcMin();
+ dish_telescope.GetDishCoefficients()->MaxRadiusInArcMin();
const double reference_frequency =
- dish_telescope.coefficients_->ReferenceFrequency();
+ dish_telescope.GetDishCoefficients()->ReferenceFrequency();
circularsymmetric::VoltagePattern vp(
- dish_telescope.coefficients_->GetFrequencies(freq), max_radius_arc_min,
- reference_frequency);
+ dish_telescope.GetDishCoefficients()->GetFrequencies(freq),
+ max_radius_arc_min, reference_frequency);
const aocommon::UVector<double> coefs_vec =
- dish_telescope.coefficients_->GetCoefficients(freq);
- vp.EvaluatePolynomial(coefs_vec, dish_telescope.coefficients_->AreInverted());
+ dish_telescope.GetDishCoefficients()->GetCoefficients(freq);
+ vp.EvaluatePolynomial(coefs_vec,
+ dish_telescope.GetDishCoefficients()->AreInverted());
vp.Render(buffer, ra, dec, pdir_ra, pdir_dec, freq);
}
diff --git a/cpp/pointresponse/mwapoint.cc b/cpp/pointresponse/mwapoint.cc
index 364751621c7dd9113314e0f9d88b996ac2b9f05a..e1d515099061b5fb8c7b7e03dca0f8f00fa559fa 100644
--- a/cpp/pointresponse/mwapoint.cc
+++ b/cpp/pointresponse/mwapoint.cc
@@ -23,17 +23,19 @@ void MWAPoint::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
if (!tile_beam_) {
tile_beam_.reset(
- new TileBeam2016(mwatelescope.ms_properties_.delays,
+ new TileBeam2016(mwatelescope.GetDelays().data(),
mwatelescope.GetOptions().frequency_interpolation,
mwatelescope.GetOptions().coeff_path));
}
- std::complex<double> gain[4];
+ std::array<std::complex<double>, 4> gains;
+
tile_beam_->ArrayResponse(ra, dec, j2000_ref_, j2000_to_hadecref_,
- j2000_to_azelgeoref_, arr_latitude_, freq, gain);
+ j2000_to_azelgeoref_, arr_latitude_, freq,
+ gains.data());
- for (size_t i = 0; i != 4; ++i) {
- *buffer = gain[i];
+ for (const auto& gain : gains) {
+ *buffer = gain;
++buffer;
}
}
@@ -57,8 +59,7 @@ void MWAPoint::SetJ200Vectors() {
std::unique_lock<std::mutex> lock(mutex_);
casacore::MEpoch time_epoch(
casacore::Quantity(time_ + 0.5 * update_interval_, "s"));
- casacore::MeasFrame frame(mwatelescope.ms_properties_.array_position,
- time_epoch);
+ casacore::MeasFrame frame(mwatelescope.GetArrayPosition(), time_epoch);
const casacore::MDirection::Ref hadec_ref(casacore::MDirection::HADEC, frame);
const casacore::MDirection::Ref azelgeo_ref(casacore::MDirection::AZELGEO,
@@ -67,7 +68,7 @@ void MWAPoint::SetJ200Vectors() {
j2000_to_hadecref_(j2000_ref_, hadec_ref);
j2000_to_azelgeoref_(j2000_ref_, azelgeo_ref);
casacore::MPosition wgs = casacore::MPosition::Convert(
- mwatelescope.ms_properties_.array_position, casacore::MPosition::WGS84)();
+ mwatelescope.GetArrayPosition(), casacore::MPosition::WGS84)();
arr_latitude_ = wgs.getValue().getLat();
}
} // namespace pointresponse
diff --git a/cpp/pointresponse/phasedarraypoint.cc b/cpp/pointresponse/phasedarraypoint.cc
index 1501d0ee7d557d70ee23dadbe23445624b290536..2fabb6dae115df397ed71afe01c476cad868daf4 100644
--- a/cpp/pointresponse/phasedarraypoint.cc
+++ b/cpp/pointresponse/phasedarraypoint.cc
@@ -11,32 +11,20 @@
#include <limits>
namespace everybeam {
+
+using telescope::PhasedArray;
+
namespace pointresponse {
PhasedArrayPoint::PhasedArrayPoint(const telescope::Telescope *telescope_ptr,
double time)
: PointResponse(telescope_ptr, time),
- beam_normalisation_mode_(
- telescope_->GetOptions().beam_normalisation_mode),
+ PhasedArrayResponse(static_cast<const PhasedArray *>(telescope_ptr)),
ra_(std::numeric_limits<double>::min()),
dec_(std::numeric_limits<double>::min()),
has_partial_itrf_update_(false),
is_local_(false),
- rotate_(true) {
- // Extract phased array specific options from ms_properties_ and
- // telescope::Options
- // TODO: code is largely a duplicate of PhasedArrayGrid
- // constructor
- const telescope::PhasedArray &phasedarray =
- dynamic_cast<const telescope::PhasedArray &>(*telescope_ptr);
- delay_dir_ = phasedarray.GetMSProperties().delay_dir;
- tile_beam_dir_ = phasedarray.GetMSProperties().tile_beam_dir;
- preapplied_beam_dir_ = phasedarray.GetMSProperties().preapplied_beam_dir;
- preapplied_correction_mode_ =
- phasedarray.GetMSProperties().preapplied_correction_mode;
- subband_frequency_ = phasedarray.GetMSProperties().subband_freq;
- use_channel_frequency_ = phasedarray.GetOptions().use_channel_frequency;
-}
+ rotate_(true) {}
void PhasedArrayPoint::Response(BeamMode beam_mode, std::complex<float> *buffer,
double ra, double dec, double freq,
@@ -50,17 +38,17 @@ void PhasedArrayPoint::Response(BeamMode beam_mode, std::complex<float> *buffer,
has_partial_itrf_update_ = false;
}
- const telescope::PhasedArray &phasedarraytelescope =
- static_cast<const telescope::PhasedArray &>(*telescope_);
+ const PhasedArray &phased_array =
+ static_cast<const PhasedArray &>(*telescope_);
aocommon::MC2x2F inverse_central_gain;
- bool apply_normalisation = CalculateBeamNormalisation(
+ const bool apply_normalisation = CalculateBeamNormalisation(
beam_mode, time_, freq, station_idx, inverse_central_gain);
const double sb_freq = use_channel_frequency_ ? freq : subband_frequency_;
const aocommon::MC2x2F gain_matrix =
- aocommon::MC2x2F(phasedarraytelescope.GetStation(station_idx)
- ->Response(beam_mode, time_, freq, dir_itrf_,
+ aocommon::MC2x2F(phased_array.GetStation(station_idx)
+ ->Response(beam_mode, time_, freq, itrf_direction_,
sb_freq, station0_, tile0_)
.Data());
@@ -71,9 +59,10 @@ void PhasedArrayPoint::Response(BeamMode beam_mode, std::complex<float> *buffer,
}
}
-aocommon::MC2x2 PhasedArrayPoint::FullResponse(size_t station_idx, double freq,
- const vector3r_t &direction,
- std::mutex *mutex) {
+aocommon::MC2x2 PhasedArrayPoint::Response(BeamMode beam_mode,
+ size_t station_idx, double freq,
+ const vector3r_t &direction,
+ std::mutex *mutex) {
if (has_time_update_) {
if (mutex != nullptr) {
// Caller takes over responsibility to be thread-safe
@@ -85,67 +74,40 @@ aocommon::MC2x2 PhasedArrayPoint::FullResponse(size_t station_idx, double freq,
has_time_update_ = false;
has_partial_itrf_update_ = true;
}
- return FullResponse(station_idx, freq, direction, station0_, tile0_);
-}
-aocommon::MC2x2 PhasedArrayPoint::FullResponse(size_t station_idx, double freq,
- const vector3r_t &direction,
- const vector3r_t &station0,
- const vector3r_t &tile0) {
- const telescope::PhasedArray &phasedarraytelescope =
- static_cast<const telescope::PhasedArray &>(*telescope_);
- const double sb_freq = use_channel_frequency_ ? freq : subband_frequency_;
- return phasedarraytelescope.GetStation(station_idx)
- ->Response(time_, freq, direction, sb_freq, station0, tile0, rotate_);
-}
+ aocommon::MC2x2F inverse_central_gain;
+ const bool apply_normalisation = CalculateBeamNormalisation(
+ beam_mode, time_, freq, station_idx, inverse_central_gain);
-aocommon::MC2x2Diag PhasedArrayPoint::ArrayFactor(size_t station_idx,
- double freq,
- const vector3r_t &direction,
- std::mutex *mutex) {
- if (has_time_update_) {
- if (mutex != nullptr) {
- // Caller takes over responsibility to be thread-safe
- UpdateITRFVectors(*mutex);
- } else {
- // Calllee assumes that caller is thread-safe
- UpdateITRFVectors(mutex_);
- }
- has_time_update_ = false;
- has_partial_itrf_update_ = true;
- }
- return ArrayFactor(station_idx, freq, direction, station0_, tile0_);
-}
+ aocommon::MC2x2 gain_matrix = UnnormalisedResponse(
+ beam_mode, station_idx, freq, direction, station0_, tile0_);
-aocommon::MC2x2Diag PhasedArrayPoint::ArrayFactor(size_t station_idx,
- double freq,
- const vector3r_t &direction,
- const vector3r_t &station0,
- const vector3r_t &tile0) {
- const telescope::PhasedArray &phasedarraytelescope =
- static_cast<const telescope::PhasedArray &>(*telescope_);
+ // Conversion to MC2x2 (double) for inverse_central_gain needed
+ return apply_normalisation
+ ? aocommon::MC2x2(inverse_central_gain.Data()) * gain_matrix
+ : gain_matrix;
+}
+aocommon::MC2x2 PhasedArrayPoint::UnnormalisedResponse(
+ BeamMode beam_mode, size_t station_idx, double freq,
+ const vector3r_t &direction, const vector3r_t &station0,
+ const vector3r_t &tile0) const {
+ const PhasedArray &phased_array =
+ static_cast<const PhasedArray &>(*telescope_);
const double sb_freq = use_channel_frequency_ ? freq : subband_frequency_;
- return phasedarraytelescope.GetStation(station_idx)
- ->ArrayFactor(time_, freq, direction, sb_freq, station0, tile0);
-}
-
-aocommon::MC2x2 PhasedArrayPoint::ElementResponse(
- size_t station_idx, double freq, const vector3r_t &direction) const {
- const telescope::PhasedArray &phasedarraytelescope =
- static_cast<const telescope::PhasedArray &>(*telescope_);
- return phasedarraytelescope.GetStation(station_idx)
- ->ComputeElementResponse(time_, freq, direction, is_local_, rotate_);
+ return phased_array.GetStation(station_idx)
+ ->Response(beam_mode, time_, freq, direction, sb_freq, station0, tile0,
+ rotate_);
}
aocommon::MC2x2 PhasedArrayPoint::ElementResponse(size_t station_idx,
double freq,
const vector3r_t &direction,
size_t element_idx) const {
- const telescope::PhasedArray &phasedarraytelescope =
- static_cast<const telescope::PhasedArray &>(*telescope_);
- return phasedarraytelescope.GetStation(station_idx)
+ const PhasedArray &phased_array =
+ static_cast<const PhasedArray &>(*telescope_);
+ return phased_array.GetStation(station_idx)
->ComputeElementResponse(time_, freq, direction, element_idx, is_local_,
rotate_);
}
@@ -168,7 +130,7 @@ void PhasedArrayPoint::UpdateITRFVectors(double ra, double dec) {
casacore::MVDirection(casacore::Quantity(ra, rad_unit),
casacore::Quantity(dec, rad_unit)),
casacore::MDirection::J2000);
- coords::SetITRFVector(itrf_converter.ToDirection(n_dir), dir_itrf_);
+ coords::SetITRFVector(itrf_converter.ToDirection(n_dir), itrf_direction_);
coords::SetITRFVector(itrf_converter.ToDirection(preapplied_beam_dir_),
diff_beam_centre_);
@@ -181,67 +143,5 @@ void PhasedArrayPoint::UpdateITRFVectors(std::mutex &mutex) {
coords::SetITRFVector(itrf_converter.ToDirection(tile_beam_dir_), tile0_);
}
-bool PhasedArrayPoint::CalculateBeamNormalisation(
- BeamMode beam_mode, double time, double frequency, size_t station_idx,
- aocommon::MC2x2F &inverse_gain) const {
- const telescope::PhasedArray &phasedarraytelescope =
- static_cast<const telescope::PhasedArray &>(*telescope_);
- if (beam_normalisation_mode_ == BeamNormalisationMode::kNone) {
- return false;
- }
-
- const double sb_freq =
- use_channel_frequency_ ? frequency : subband_frequency_;
-
- // if the normalisation mode is kPreApplied, but no beam correction was pre
- // applied then there is nothing to do
- if (beam_normalisation_mode_ == BeamNormalisationMode::kPreApplied &&
- preapplied_correction_mode_ == BeamMode::kNone) {
- return false;
- }
- // If the normalisation mode is kPreApplied, or kPreAppliedOrFull and the
- // fallback to Full is not needed then the response for the diff_beam_centre_
- // with preapplied_correction_mode_ needs to be computed
- if (beam_normalisation_mode_ == BeamNormalisationMode::kPreApplied ||
- (beam_normalisation_mode_ == BeamNormalisationMode::kPreAppliedOrFull &&
- preapplied_correction_mode_ != BeamMode::kNone)) {
- inverse_gain = aocommon::MC2x2F(
- phasedarraytelescope.GetStation(station_idx)
- ->Response(preapplied_correction_mode_, time, frequency,
- diff_beam_centre_, sb_freq, station0_, tile0_)
- .Data());
- } else {
- // in all other cases the response in for the reference direction with
- // beam_mode is needed
- inverse_gain = aocommon::MC2x2F(
- phasedarraytelescope.GetStation(station_idx)
- ->Response(beam_mode, time, frequency, diff_beam_centre_, sb_freq,
- station0_, tile0_)
- .Data());
- }
-
- switch (beam_normalisation_mode_) {
- case BeamNormalisationMode::kFull:
- case BeamNormalisationMode::kPreApplied:
- case BeamNormalisationMode::kPreAppliedOrFull:
- if (!inverse_gain.Invert()) {
- inverse_gain = aocommon::MC2x2F::Zero();
- }
- break;
- case BeamNormalisationMode::kAmplitude: {
- const float amplitude_inv = 1.0 / std::sqrt(0.5 * Norm(inverse_gain));
- inverse_gain[0] = std::isfinite(amplitude_inv) ? amplitude_inv : 0.0;
- inverse_gain[1] = 0.0;
- inverse_gain[2] = 0.0;
- inverse_gain[3] = std::isfinite(amplitude_inv) ? amplitude_inv : 0.0;
- break;
- }
- default:
- throw std::runtime_error("Invalid beam normalisation mode here");
- }
-
- return true;
-}
-
} // namespace pointresponse
} // namespace everybeam
diff --git a/cpp/pointresponse/phasedarraypoint.h b/cpp/pointresponse/phasedarraypoint.h
index 93d118c0ed701dbdfda42972da42216c68ce3718..f9c972cd49980dc9e3dbd34e94c7a42466d7be4a 100644
--- a/cpp/pointresponse/phasedarraypoint.h
+++ b/cpp/pointresponse/phasedarraypoint.h
@@ -11,6 +11,7 @@
#include "../common/types.h"
#include "../correctionmode.h"
#include "../beamnormalisationmode.h"
+#include "../phasedarrayresponse.h"
#include <aocommon/matrix2x2.h>
#include <casacore/measures/Measures/MDirection.h>
@@ -18,7 +19,7 @@
namespace everybeam {
namespace pointresponse {
-class PhasedArrayPoint : public PointResponse {
+class PhasedArrayPoint : public PointResponse, protected PhasedArrayResponse {
public:
PhasedArrayPoint(const telescope::Telescope *telescope_ptr, double time);
@@ -44,30 +45,31 @@ class PhasedArrayPoint : public PointResponse {
double ra, double dec, double freq, size_t station_idx,
size_t field_id) final override;
- aocommon::MC2x2 FullResponse(size_t station_idx, double freq,
- const vector3r_t &direction,
- std::mutex *mutex) final override;
-
/**
- * @brief Convenience function for the python bindings
+ * @brief Compute beam response. Optional beam normalisation is
+ * done in this function
+ *
+ * @param beam_mode BeamMode, can be any of kNone, kFull, kArrayFactor or
+ * kElement
+ * @param station_idx Station index for which to compute the beam response.
+ * @param freq Freq [Hz]
+ * @param direction Direction in ITRF
+ * @param mutex mutex. When provided, the caller keeps control over
+ * thread-safety. If not provided, the internal mutex will be used and the
+ * caller is assumed to be thread-safe.
+ * @return aocommon::MC2x2
*/
- aocommon::MC2x2 FullResponse(size_t station_idx, double freq,
- const vector3r_t &direction,
- const vector3r_t &station0,
- const vector3r_t &tile0);
-
- aocommon::MC2x2Diag ArrayFactor(size_t station_idx, double freq,
- const vector3r_t &direction,
- std::mutex *mutex) final override;
+ aocommon::MC2x2 Response(BeamMode beam_mode, size_t station_idx, double freq,
+ const vector3r_t &direction,
+ std::mutex *mutex) final override;
/**
- * @brief Compute the array factor given a direction, station0 direction
- * and tile0 direction. Method is used in the python bindings.
+ * @brief Compute the unnormalised response.
*/
- aocommon::MC2x2Diag ArrayFactor(size_t station_idx, double freq,
- const vector3r_t &direction,
- const vector3r_t &station0,
- const vector3r_t &tile0);
+ aocommon::MC2x2 UnnormalisedResponse(BeamMode beam_mode, size_t station_idx,
+ double freq, const vector3r_t &direction,
+ const vector3r_t &station0,
+ const vector3r_t &tile0) const;
/**
* @brief Convenience method for computing the element response, for a
@@ -86,10 +88,6 @@ class PhasedArrayPoint : public PointResponse {
const vector3r_t &direction,
size_t element_idx) const;
- aocommon::MC2x2 ElementResponse(
- size_t station_idx, double freq,
- const vector3r_t &direction) const final override;
-
/**
* @brief Method for computing the ITRF-vectors, given ra, dec position in
* radians and using the cached \param time ((MJD(UTC), s))
@@ -110,16 +108,6 @@ class PhasedArrayPoint : public PointResponse {
void SetParalacticRotation(bool rotate) { rotate_ = rotate; }
bool GetParalacticRotation() const { return rotate_; };
- protected:
- casacore::MDirection delay_dir_, tile_beam_dir_;
- vector3r_t station0_, tile0_, dir_itrf_, diff_beam_centre_;
- bool use_channel_frequency_;
- bool use_differential_beam_;
- casacore::MDirection preapplied_beam_dir_;
- CorrectionMode preapplied_correction_mode_;
- BeamNormalisationMode beam_normalisation_mode_;
- double subband_frequency_;
-
private:
/**
* @brief Update ITRF coordinates for reference station and reference tile
@@ -128,10 +116,7 @@ class PhasedArrayPoint : public PointResponse {
*/
void UpdateITRFVectors(std::mutex &mutex);
- bool CalculateBeamNormalisation(BeamMode beam_mode, double time,
- double frequency, size_t station_idx,
- aocommon::MC2x2F &inverse_gain) const;
-
+ vector3r_t itrf_direction_;
double ra_, dec_;
std::mutex mutex_;
diff --git a/cpp/pointresponse/pointresponse.h b/cpp/pointresponse/pointresponse.h
index d2e3a7045218b0302174cba128840360c1a7d558..8e302f723b282123fdb89c9a887fb8528c9318b1 100644
--- a/cpp/pointresponse/pointresponse.h
+++ b/cpp/pointresponse/pointresponse.h
@@ -20,12 +20,11 @@ namespace everybeam {
namespace pointresponse {
/**
- * @brief Virtual base class to compute the point response
- *
+ * @brief Virtual base class to compute the point response.
*/
class PointResponse {
public:
- virtual ~PointResponse() {}
+ virtual ~PointResponse() = default;
/**
* @brief Update the (cached) time
@@ -58,44 +57,6 @@ class PointResponse {
*/
bool HasTimeUpdate() const { return has_time_update_; }
- /**
- * @brief See FullResponse.
- */
- [[deprecated("Use FullResponse() instead.")]] void CalculateStation(
- std::complex<float>* response_matrix, double ra, double dec, double freq,
- size_t station_idx, size_t field_id) {
- FullResponse(response_matrix, ra, dec, freq, station_idx, field_id);
- };
-
- /**
- * @brief See FullResponseAllStations.
- */
- [[deprecated("Use FullResponseAllStations() instead.")]] void
- CalculateAllStations(std::complex<float>* response_matrices, double ra,
- double dec, double freq, size_t field_id) {
- FullResponseAllStations(response_matrices, ra, dec, freq, field_id);
- };
-
- /**
- * @brief Get beam response for a given station at a prescribed ra, dec
- * position.
- *
- * @param response_matrix Buffer with a size of 4 complex floats to receive
- * the beam response
- * @param ra Right ascension (rad)
- * @param dec Declination (rad)
- * @param freq Frequency (Hz)
- * @param station_id Station index, corresponding to measurement set antenna
- * index.
- * @param field_id Field index as used in the measurement set
- */
- virtual void FullResponse(std::complex<float>* response_matrix, double ra,
- double dec, double freq, size_t station_id,
- size_t field_id) {
- Response(BeamMode::kFull, response_matrix, ra, dec, freq, station_id,
- field_id);
- }
-
/**
* @brief Get beam response for a given station at a prescribed ra, dec
* position.
@@ -117,74 +78,8 @@ class PointResponse {
size_t field_id) = 0;
/**
- * @brief Get the full beam response for a station, given a pointing direction
- * in ITRF coordinates
- *
- * @param beam_mode Selects beam mode (element, array factor or full)
- * @param station_idx Station index
- * @param freq Frequency (Hz)
- * @param direction Direction in ITRF
- * @param mutex Optional mutex. When provided, the caller keeps control over
- * thread-safety. If not provided, the internal mutex will be used and the
- * caller is assumed to be thread-safe.
- * @return aocommon::MC2x2
- */
- aocommon::MC2x2 Response(BeamMode beam_mode, size_t station_idx, double freq,
- const vector3r_t& direction,
- std::mutex* mutex = nullptr) {
- switch (beam_mode) {
- case BeamMode::kNone:
- return aocommon::MC2x2::Unity();
- case BeamMode::kFull:
- return FullResponse(station_idx, freq, direction, mutex);
- case BeamMode::kElement:
- return ElementResponse(station_idx, freq, direction);
- case BeamMode::kArrayFactor:
- return ArrayFactor(station_idx, freq, direction, mutex);
- default:
- throw std::runtime_error("Invalid beam mode.");
- }
- }
-
- /**
- * @brief Get the full beam response for a station, given a pointing direction
- * in ITRF coordinates
- *
- * @param station_idx Station index
- * @param freq Frequency (Hz)
- * @param direction Direction in ITRF
- * @param mutex Optional mutex. When provided, the caller keeps control over
- * thread-safety. If not provided, the internal mutex will be used and the
- * caller is assumed to be thread-safe.
- * @return aocommon::MC2x2
- */
- virtual aocommon::MC2x2 FullResponse(size_t station_idx, double freq,
- const vector3r_t& direction,
- std::mutex* mutex = nullptr) {
- throw std::runtime_error("Not yet implemented");
- }
-
- /**
- * @brief Same as FullResponse, but now iterate over all stations in
- * MS.
- *
- * @param response_matrices Buffer with a size of 4 * nr_stations complex
- * floats to receive the beam response
- * @param ra Right ascension (rad)
- * @param dec Declination (rad)
- * @param freq Frequency (Hz)
- * @param field_id Field index as used in the measurement set
- */
- virtual void FullResponseAllStations(std::complex<float>* response_matrices,
- double ra, double dec, double freq,
- size_t field_id) {
- ResponseAllStations(BeamMode::kFull, response_matrices, ra, dec, freq,
- field_id);
- }
-
- /**
- * @brief Same as FullResponse, but now iterate over all stations in
- * MS.
+ * @brief Same as Response, but now iterate over all stations in
+ * measurement set.
*
* @param beam_mode Selects beam mode (element, array factor or full)
* @param response_matrices Buffer with a size of 4 * nr_stations complex
@@ -205,7 +100,7 @@ class PointResponse {
}
/**
- * @brief Get the array factor for a station, given a pointing direction
+ * @brief Get the beam response for a station, given a pointing direction
* in ITRF coordinates
*
* @param station_idx Station index
@@ -214,16 +109,11 @@ class PointResponse {
* @param mutex Optional mutex. When provided, the caller keeps control over
* thread-safety. If not provided, the internal mutex will be used and the
* caller is assumed to be thread-safe.
- * @return aocommon::MC2x2Diag
+ * @return aocommon::MC2x2
*/
- virtual aocommon::MC2x2Diag ArrayFactor(size_t station_idx, double freq,
- const vector3r_t& direction,
- std::mutex* mutex = nullptr) {
- throw std::runtime_error("Not yet implemented");
- };
-
- virtual aocommon::MC2x2 ElementResponse(size_t station_idx, double freq,
- const vector3r_t& direction) const {
+ virtual aocommon::MC2x2 Response(BeamMode beam_mode, size_t station_idx,
+ double freq, const vector3r_t& direction,
+ std::mutex* mutex = nullptr) {
throw std::runtime_error("Not yet implemented");
}
diff --git a/cpp/telescope/dish.cc b/cpp/telescope/dish.cc
index 2efa0f76e02e371ca48606520dfe687dd97bbcc1..532bad464b3592367c5429376c524d01948e1f12 100644
--- a/cpp/telescope/dish.cc
+++ b/cpp/telescope/dish.cc
@@ -16,7 +16,7 @@ using everybeam::telescope::Dish;
Dish::Dish(const casacore::MeasurementSet &ms,
std::unique_ptr<circularsymmetric::Coefficients> coefficients,
const Options &options)
- : Telescope(ms, options), coefficients_(std::move(coefficients)) {
+ : Telescope(ms, options), dish_coefficients_(std::move(coefficients)) {
casacore::MSField field_table = ms.field();
casacore::ArrayColumn<double> pointing_dir_col(
field_table, casacore::MSField::columnName(casacore::MSField::DELAY_DIR));
@@ -25,7 +25,7 @@ Dish::Dish(const casacore::MeasurementSet &ms,
casacore::Array<double> pdir = pointing_dir_col(field_id);
double pdir_ra = *pdir.cbegin();
double pdir_dec = *(pdir.cbegin() + 1);
- ms_properties_.field_pointing.emplace_back(pdir_ra, pdir_dec);
+ field_pointing_.emplace_back(pdir_ra, pdir_dec);
}
SetIsTimeRelevant(false);
}
diff --git a/cpp/telescope/dish.h b/cpp/telescope/dish.h
index c91880dfc54a2f6570fa5c06162adc90c64aaadd..3e93bdc69e72f1f2a6bc76d62029185ee2b308f2 100644
--- a/cpp/telescope/dish.h
+++ b/cpp/telescope/dish.h
@@ -17,12 +17,10 @@ namespace everybeam {
namespace griddedresponse {
class DishGrid;
-class GriddedResponse;
} // namespace griddedresponse
namespace pointresponse {
class PointResponse;
-class DishPoint;
} // namespace pointresponse
namespace telescope {
@@ -32,9 +30,6 @@ namespace telescope {
* response.
*/
class Dish final : public Telescope {
- friend class griddedresponse::DishGrid;
- friend class pointresponse::DishPoint;
-
public:
Dish(const casacore::MeasurementSet &ms,
std::unique_ptr<circularsymmetric::Coefficients> coefficients,
@@ -46,12 +41,24 @@ class Dish final : public Telescope {
std::unique_ptr<pointresponse::PointResponse> GetPointResponse(
double time) const override;
+ /**
+ * @brief Get (ra, dec) pointings of fields.
+ *
+ * @return std::vector<std::pair<double, double>> Vector of size number of
+ * fields, and (ra, dec) pointings as entries.
+ */
+ const std::vector<std::pair<double, double>> &GetFieldPointing() const {
+ return field_pointing_;
+ }
+
+ const circularsymmetric::Coefficients *GetDishCoefficients() const {
+ return dish_coefficients_.get();
+ }
+
private:
- struct MSProperties {
- std::vector<std::pair<double, double>> field_pointing;
- };
- MSProperties ms_properties_;
- std::unique_ptr<circularsymmetric::Coefficients> coefficients_;
+ std::unique_ptr<circularsymmetric::Coefficients> dish_coefficients_;
+ /// Store ra, dec pointing per field id from measurement set
+ std::vector<std::pair<double, double>> field_pointing_;
};
} // namespace telescope
} // namespace everybeam
diff --git a/cpp/telescope/mwa.cc b/cpp/telescope/mwa.cc
index e3f50d6a2727f6f409d44e05909e6b9614c935dd..b6f12e38d6c77eabd007d22a02622c59902f1aeb 100644
--- a/cpp/telescope/mwa.cc
+++ b/cpp/telescope/mwa.cc
@@ -16,19 +16,24 @@ using everybeam::telescope::MWA;
MWA::MWA(const casacore::MeasurementSet &ms, const Options &options)
: Telescope(ms, options) {
- if (GetNrStations() == 0) throw std::runtime_error("No antennae in set");
+ if (GetNrStations() == 0) {
+ throw std::runtime_error("No antennae in set");
+ }
casacore::MSAntenna antenna(ms.antenna());
casacore::MPosition::ScalarColumn antenna_pos_col(
antenna, antenna.columnName(casacore::MSAntennaEnums::POSITION));
- ms_properties_.array_position = antenna_pos_col(0);
+ array_position_ = antenna_pos_col(0);
casacore::Table mwa_tile_pointing =
ms.keywordSet().asTable("MWA_TILE_POINTING");
casacore::ArrayColumn<int> delays_col(mwa_tile_pointing, "DELAYS");
casacore::Array<int> delays_arr = delays_col(0);
casacore::Array<int>::contiter delays_arr_ptr = delays_arr.cbegin();
- for (int i = 0; i != 16; ++i) ms_properties_.delays[i] = delays_arr_ptr[i];
+ for (auto &delay : delays_) {
+ delay = *delays_arr_ptr;
+ ++delays_arr_ptr;
+ }
}
std::unique_ptr<GriddedResponse> MWA::GetGriddedResponse(
@@ -38,7 +43,6 @@ std::unique_ptr<GriddedResponse> MWA::GetGriddedResponse(
}
std::unique_ptr<PointResponse> MWA::GetPointResponse(double time) const {
- // Get and return PointResponse ptr
std::unique_ptr<PointResponse> point_response(new MWAPoint(this, time));
return point_response;
}
\ No newline at end of file
diff --git a/cpp/telescope/mwa.h b/cpp/telescope/mwa.h
index dadb50e2f902e800a0c57967368f25c7e355da52..9e9182d7dcc2b2c5f0f817c002ce52016f9c5ffd 100644
--- a/cpp/telescope/mwa.h
+++ b/cpp/telescope/mwa.h
@@ -14,21 +14,16 @@
namespace everybeam {
namespace griddedresponse {
-class MWAGrid;
class GriddedResponse;
} // namespace griddedresponse
namespace pointresponse {
-class MWAPoint;
class PointResponse;
} // namespace pointresponse
namespace telescope {
class MWA final : public Telescope {
- friend class griddedresponse::MWAGrid;
- friend class pointresponse::MWAPoint;
-
public:
/**
* @brief Construct a new MWA object
@@ -45,12 +40,12 @@ class MWA final : public Telescope {
std::unique_ptr<pointresponse::PointResponse> GetPointResponse(
double time) const override;
+ casacore::MPosition GetArrayPosition() const { return array_position_; }
+ const std::array<double, 16> &GetDelays() const { return delays_; }
+
private:
- struct MSProperties {
- double delays[16];
- casacore::MPosition array_position;
- };
- MSProperties ms_properties_;
+ casacore::MPosition array_position_;
+ std::array<double, 16> delays_;
};
} // namespace telescope
diff --git a/python/test/test_pybindings.py b/python/test/test_pybindings.py
index aefd6cceed257d688a42bd4d1aec3db22ca7b17c..8b9665b379d268060f17a1052f588e6dfe33119b 100644
--- a/python/test/test_pybindings.py
+++ b/python/test/test_pybindings.py
@@ -170,13 +170,6 @@ def test_coordinate_system():
def test_lofar(ref, differential_beam):
ms_path = os.path.join(DATADIR, ref["filename"])
telescope = load_telescope(ms_path, use_differential_beam=differential_beam)
- a = telescope.gridded_response(
- ref["coordinate_system"],
- ref["time"],
- ref["freq"],
- ref["station_id"],
- field_index=10,
- )
assert isinstance(telescope, LOFAR)
time = ref["time"]
diff --git a/python/wrappers/pytelescope.cc b/python/wrappers/pytelescope.cc
index e8d9cf750b4c65f21f9c24270552a21c5b0efb53..ad888be6bf4ea685a25cd78126896242075c0b2c 100644
--- a/python/wrappers/pytelescope.cc
+++ b/python/wrappers/pytelescope.cc
@@ -438,8 +438,8 @@ void init_telescope(py::module &m) {
(self.GetOptions().beam_normalisation_mode ==
BeamNormalisationMode::kPreApplied)
? aocommon::MC2x2::Unity()
- : phased_array_point.FullResponse(idx, freq, direction,
- &mutex);
+ : phased_array_point.Response(BeamMode::kFull, idx, freq,
+ direction, &mutex);
return cast_matrix(response);
},
R"pbdoc(
@@ -597,9 +597,12 @@ void init_telescope(py::module &m) {
PhasedArrayPoint &phased_array_point =
static_cast<PhasedArrayPoint &>(*point_response);
phased_array_point.SetParalacticRotation(rotate);
- const aocommon::MC2x2 response = phased_array_point.FullResponse(
- idx, freq, direction, station0, tile0);
+ aocommon::MC2x2 response = phased_array_point.UnnormalisedResponse(
+ BeamMode::kFull, idx, freq, direction, station0, tile0);
+ // We can't delegate the normalisation to Response, in this case.
+ // Since the diff_beam_centre direction would not be computed
+ // correctly.
if (self.GetOptions().beam_normalisation_mode ==
BeamNormalisationMode::kPreApplied) {
vector3r_t diff_beam_centre;
@@ -609,13 +612,15 @@ void init_telescope(py::module &m) {
diff_beam_centre);
aocommon::MC2x2 response_diff_beam =
- phased_array_point.FullResponse(idx, freq, diff_beam_centre,
- station0, tile0);
+ phased_array_point.UnnormalisedResponse(
+ BeamMode::kFull, idx, freq, diff_beam_centre, station0,
+ tile0);
apply_differential_beam(response_diff_beam, response);
return cast_matrix(response_diff_beam);
} else {
return cast_matrix(response);
}
+ return cast_matrix(response);
},
R"pbdoc(
Get station response in user-specified direction
@@ -809,20 +814,30 @@ void init_telescope(py::module &m) {
phased_array_point.SetUseLocalCoordinateSystem(is_local);
phased_array_point.SetParalacticRotation(rotate);
+ // Avoid any beam normalisations, so compute station0 and tile0
+ // manually
+ ITRFConverter itrf_converter(time);
+ vector3r_t station0;
+ vector3r_t tile0;
+ SetITRFVector(itrf_converter.ToDirection(self.GetDelayDirection()),
+ station0);
+ SetITRFVector(
+ itrf_converter.ToDirection(self.GetTileBeamDirection()), tile0);
const aocommon::MC2x2 response =
- phased_array_point.ElementResponse(idx, freq, direction);
+ phased_array_point.UnnormalisedResponse(
+ BeamMode::kElement, idx, freq, direction, station0, tile0);
if (self.GetOptions().beam_normalisation_mode ==
BeamNormalisationMode::kPreApplied) {
vector3r_t diff_beam_centre;
- ITRFConverter itrf_converter(time);
SetITRFVector(
itrf_converter.ToDirection(self.GetPreappliedBeamDirection()),
diff_beam_centre);
aocommon::MC2x2 response_diff_beam =
- phased_array_point.ElementResponse(idx, freq,
- diff_beam_centre);
+ phased_array_point.UnnormalisedResponse(
+ BeamMode::kElement, idx, freq, diff_beam_centre, station0,
+ tile0);
apply_differential_beam(response_diff_beam, response);
return cast_matrix(response_diff_beam);
} else {
@@ -875,8 +890,10 @@ void init_telescope(py::module &m) {
static_cast<PhasedArrayPoint &>(*point_response);
// Diagonal to 2x2 matrix
- const aocommon::MC2x2 response(phased_array_point.ArrayFactor(
- idx, freq, direction, station0, tile0));
+ const aocommon::MC2x2 response(
+ phased_array_point.UnnormalisedResponse(BeamMode::kArrayFactor,
+ idx, freq, direction,
+ station0, tile0));
return cast_matrix(response);
},
R"pbdoc(