Skip to content
Snippets Groups Projects
Commit 424660c9 authored by Andre Offringa's avatar Andre Offringa
Browse files

Average over squared Mueller matrices

parent b2bc2b93
No related branches found
Tags v0.6.0
1 merge request!351Average over squared Mueller matrices
Pipeline #82068 passed with warnings
Showing with 183 additions and 139 deletions
......@@ -9,7 +9,7 @@ include(CheckFunctionExists)
#------------------------------------------------------------------------------
# Set version name and project number
set(EVERYBEAM_VERSION 0.5.8) # Keep in sync with `pyproject.toml` file
set(EVERYBEAM_VERSION 0.6.0) # Keep in sync with `pyproject.toml` file
if(EVERYBEAM_VERSION MATCHES "^([0-9]+)\\.([0-9]+)\\.([0-9]+)")
set(EVERYBEAM_VERSION_MAJOR "${CMAKE_MATCH_1}")
set(EVERYBEAM_VERSION_MINOR "${CMAKE_MATCH_2}")
......
......@@ -28,10 +28,10 @@ void AartfaacGrid::ResponseAllStations(BeamMode beam_mode,
}
}
void AartfaacGrid::MakeIntegratedSnapshot(
void AartfaacGrid::MakeIntegratedCorrectionSnapshot(
BeamMode beam_mode, std::vector<aocommon::HMC4x4>& matrices, double time,
double frequency, size_t field_id,
const double* baseline_weights_interval) {
double frequency, size_t field_id, const double* baseline_weights_interval,
bool square_mueller) {
const size_t n_stations = telescope_->GetNrStations();
aocommon::UVector<std::complex<float>> buffer_undersampled(
GetStationBufferSize(n_stations));
......@@ -53,9 +53,15 @@ void AartfaacGrid::MakeIntegratedSnapshot(
// Mueller matrix constant for all baselines, so just compute once for
// each individual pixel
matrices[y * width_ + x] =
aocommon::HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A) *
snapshot_weight;
if (square_mueller) {
matrices[y * width_ + x] =
aocommon::HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A)
.Square() *
snapshot_weight;
} else {
aocommon::HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A) *
snapshot_weight;
}
}
}
}
......
......@@ -39,10 +39,10 @@ class [[gnu::visibility("default")]] AartfaacGrid final
private:
// Override MakeIntegratedSnapshot for efficiency.
// Implementation is similar to MWAGrid::MakeIntegratedSnapshot
void MakeIntegratedSnapshot(BeamMode beam_mode,
std::vector<aocommon::HMC4x4> & matrices,
double time, double frequency, size_t field_id,
const double* baseline_weights_interval) override;
void MakeIntegratedCorrectionSnapshot(
BeamMode beam_mode, std::vector<aocommon::HMC4x4> & matrices, double time,
double frequency, size_t field_id,
const double* baseline_weights_interval, bool square_mueller) override;
};
} // namespace griddedresponse
} // namespace everybeam
......
......@@ -26,15 +26,17 @@ class [[gnu::visibility("default")]] AiryGrid final : public GriddedResponse {
double time, double frequency, size_t field_id)
override;
void IntegratedResponse(
void IntegratedCorrection(
BeamMode beam_mode, float* buffer,
[[maybe_unused]] const std::vector<double>& time_array, double frequency,
size_t field_id, size_t undersampling_factor,
[[maybe_unused]] const std::vector<double>& baseline_weights) override {
[[maybe_unused]] const std::vector<double>& baseline_weights,
bool square_mueller) override {
// The integrated response of an Airy disk telescope is independent of time,
// so call IntegratedResponse as if it were one time step
GriddedResponse::IntegratedResponse(beam_mode, buffer, 0.0, frequency,
field_id, undersampling_factor, {0.0});
GriddedResponse::IntegratedCorrection(beam_mode, buffer, 0.0, frequency,
field_id, undersampling_factor, {0.0},
square_mueller);
}
bool PerformUndersampling() const override { return false; }
......
......@@ -41,10 +41,10 @@ void DishGrid::Response([[maybe_unused]] BeamMode beam_mode,
l_shift_, m_shift_, frequency);
}
void DishGrid::IntegratedResponse(
void DishGrid::IntegratedCorrection(
BeamMode /* beam_mode */, float* buffer, double, double frequency,
size_t field_id, size_t undersampling_factor,
const std::vector<double>& /*baseline_weights*/) {
const std::vector<double>& /*baseline_weights*/, bool square_mueller) {
// Copy coordinate members
const size_t width_original = width_;
const size_t height_original = height_;
......@@ -60,6 +60,11 @@ void DishGrid::IntegratedResponse(
const size_t npixels = width_ * height_;
std::vector<HMC4x4> matrices(npixels, HMC4x4::Zero());
MakeIntegratedDishSnapshot(matrices, frequency, field_id);
if (square_mueller) {
for (HMC4x4& m : matrices) {
m = m.Square();
}
}
DoFFTResampling(buffer, width_, height_, width_original, height_original,
matrices);
......
......@@ -30,21 +30,22 @@ class [[gnu::visibility("default")]] DishGrid : public GriddedResponse {
HomogeneousAllStationResponse(beam_mode, buffer, time, frequency, field_id);
}
void IntegratedResponse(
BeamMode beam_mode, float* buffer, double time, double frequency,
size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights) final override;
void IntegratedCorrection(BeamMode beam_mode, float* buffer, double time,
double frequency, size_t field_id,
size_t undersampling_factor,
const std::vector<double>& baseline_weights,
bool square_mueller) final override;
void IntegratedResponse(
void IntegratedCorrection(
BeamMode beam_mode, float* buffer,
[[maybe_unused]] const std::vector<double>& time_array, double frequency,
size_t field_id, size_t undersampling_factor,
[[maybe_unused]] const std::vector<double>& baseline_weights)
final override {
[[maybe_unused]] const std::vector<double>& baseline_weights,
bool square_mueller) final override {
// Time does not play a role in the integrated response of a dish telescope,
// so call IntegratedResponse as if it were one time step
IntegratedResponse(beam_mode, buffer, 0.0, frequency, field_id,
undersampling_factor, {0.0});
IntegratedCorrection(beam_mode, buffer, 0.0, frequency, field_id,
undersampling_factor, {0.0}, square_mueller);
}
bool PerformUndersampling() const final override { return false; }
......
......@@ -23,10 +23,10 @@ using aocommon::UVector;
namespace everybeam {
namespace griddedresponse {
void GriddedResponse::IntegratedResponse(
void GriddedResponse::IntegratedCorrection(
BeamMode beam_mode, float* buffer, double time, double frequency,
size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights) {
const std::vector<double>& baseline_weights, bool square_mueller) {
const size_t nstations = telescope_->GetNrStations();
const size_t nbaselines = nstations * (nstations + 1) / 2;
if (baseline_weights.size() != nbaselines) {
......@@ -52,8 +52,9 @@ void GriddedResponse::IntegratedResponse(
// Init (Hermitian) Mueller matrix for every pixel in the coarse grid
size_t npixels = width_ * height_;
std::vector<HMC4x4> matrices(npixels, HMC4x4::Zero());
MakeIntegratedSnapshot(beam_mode, matrices, time, frequency, field_id,
baseline_weights.data());
MakeIntegratedCorrectionSnapshot(beam_mode, matrices, time, frequency,
field_id, baseline_weights.data(),
square_mueller);
for (HMC4x4& matrix : matrices) {
matrix /= baseline_total_weight;
......@@ -69,10 +70,10 @@ void GriddedResponse::IntegratedResponse(
dm_ = dm_original;
}
std::vector<HMC4x4> GriddedResponse::UndersampledIntegratedResponse(
std::vector<HMC4x4> GriddedResponse::UndersampledIntegratedCorrection(
BeamMode beam_mode, const std::vector<double>& time_array, double frequency,
size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights) {
const std::vector<double>& baseline_weights, bool square_mueller) {
const size_t nstations = telescope_->GetNrStations();
const size_t nbaselines = nstations * (nstations + 1) / 2;
if (baseline_weights.size() != time_array.size() * nbaselines) {
......@@ -97,9 +98,9 @@ std::vector<HMC4x4> GriddedResponse::UndersampledIntegratedResponse(
std::vector<HMC4x4> matrices(npixels, HMC4x4::Zero());
for (std::size_t tstep = 0; tstep != time_array.size(); ++tstep) {
MakeIntegratedSnapshot(beam_mode, matrices, time_array[tstep], frequency,
field_id,
baseline_weights.data() + tstep * nbaselines);
MakeIntegratedCorrectionSnapshot(
beam_mode, matrices, time_array[tstep], frequency, field_id,
baseline_weights.data() + tstep * nbaselines, square_mueller);
}
for (HMC4x4& matrix : matrices) {
......@@ -115,13 +116,14 @@ std::vector<HMC4x4> GriddedResponse::UndersampledIntegratedResponse(
return matrices;
}
void GriddedResponse::IntegratedResponse(
void GriddedResponse::IntegratedCorrection(
BeamMode beam_mode, float* destination,
const std::vector<double>& time_array, double frequency, size_t field_id,
size_t undersampling_factor, const std::vector<double>& baseline_weights) {
const std::vector<HMC4x4> matrices =
UndersampledIntegratedResponse(beam_mode, time_array, frequency, field_id,
undersampling_factor, baseline_weights);
size_t undersampling_factor, const std::vector<double>& baseline_weights,
bool square_mueller) {
const std::vector<HMC4x4> matrices = UndersampledIntegratedCorrection(
beam_mode, time_array, frequency, field_id, undersampling_factor,
baseline_weights, square_mueller);
const size_t undersampled_width = width_ / undersampling_factor;
const size_t undersampled_height = height_ / undersampling_factor;
......@@ -130,7 +132,7 @@ void GriddedResponse::IntegratedResponse(
height_, matrices);
}
void GriddedResponse::UpsampleResponse(
void GriddedResponse::UpsampleCorrection(
float* destination, size_t element_index, size_t width, size_t height,
const std::vector<aocommon::HMC4x4>& undersampled_beam,
size_t undersampling_factor) {
......@@ -157,10 +159,10 @@ void GriddedResponse::UpsampleResponse(
}
}
void GriddedResponse::MakeIntegratedSnapshot(
void GriddedResponse::MakeIntegratedCorrectionSnapshot(
BeamMode beam_mode, std::vector<HMC4x4>& matrices, double time,
double frequency, size_t field_id,
const double* baseline_weights_interval) {
double frequency, size_t field_id, const double* baseline_weights_interval,
bool square_mueller) {
const size_t nstations = telescope_->GetNrStations();
UVector<std::complex<float>> buffer_undersampled(
GetStationBufferSize(nstations));
......@@ -174,28 +176,41 @@ void GriddedResponse::MakeIntegratedSnapshot(
HMC4x4 gain = HMC4x4::Zero();
for (size_t s1 = 0; s1 != nstations; ++s1) {
size_t offset_s1 = (s1 * npixels + y * width_ + x) * 4;
MC2x2 A(buffer_undersampled[offset_s1],
buffer_undersampled[offset_s1 + 1],
buffer_undersampled[offset_s1 + 2],
buffer_undersampled[offset_s1 + 3]);
const MC2x2 A(buffer_undersampled[offset_s1],
buffer_undersampled[offset_s1 + 1],
buffer_undersampled[offset_s1 + 2],
buffer_undersampled[offset_s1 + 3]);
for (size_t s2 = s1; s2 != nstations; ++s2) {
size_t offset_s2 = offset_s1 + (s2 - s1) * npixels * 4;
MC2x2 B(buffer_undersampled[offset_s2],
buffer_undersampled[offset_s2 + 1],
buffer_undersampled[offset_s2 + 2],
buffer_undersampled[offset_s2 + 3]);
// Compute Mueller matrix and apply vec trick, see
// https://en.wikipedia.org/wiki/Kronecker_product#Matrix_equations
MC4x4 baseline_mueller =
MC4x4::KroneckerProduct(B.HermTranspose().Transpose(), A) +
MC4x4::KroneckerProduct(A.HermTranspose().Transpose(), B);
const MC2x2 B(buffer_undersampled[offset_s2],
buffer_undersampled[offset_s2 + 1],
buffer_undersampled[offset_s2 + 2],
buffer_undersampled[offset_s2 + 3]);
HMC4x4 baseline_mueller;
if (square_mueller) {
// Compute squared Mueller matrix by applying vec trick, see
// https://en.wikipedia.org/wiki/Kronecker_product#Matrix_equations
// Also uses the rule (B* (x) A)^H (B* (x) A) = (B^H B)* (x) (A^H A)
// where (x) is the Kronecker product.
baseline_mueller = HMC4x4(MC4x4::KroneckerProduct(
(B.HermTranspose() * B).Conjugate(),
A.HermTranspose() * A)) +
HMC4x4(MC4x4::KroneckerProduct(
(A.HermTranspose() * A).Conjugate(),
B.HermTranspose() * B));
} else {
// This is a simplification: this matrix is not necessarily
// Hermitian. For homogeneous arrays it is accurate.
baseline_mueller =
HMC4x4(MC4x4::KroneckerProduct(B.Conjugate(), A)) +
HMC4x4(MC4x4::KroneckerProduct(A.Conjugate(), B));
}
// Note: keep scalar factor in brackets to avoid a redundant
// element-wise multiplication of the HMC4x4 matrix. Factor 0.5 is to
// account for Kronecker delta additions.
gain += HMC4x4(baseline_mueller) *
(0.5 * baseline_weights_interval[index]);
gain += baseline_mueller * (0.5 * baseline_weights_interval[index]);
++index;
}
}
......
......@@ -33,7 +33,7 @@ namespace griddedresponse {
*/
class [[gnu::visibility("default")]] GriddedResponse {
public:
virtual ~GriddedResponse() {}
virtual ~GriddedResponse() = default;
/**
* @brief Compute the beam for a single station, given a prescribed beam mode.
......@@ -69,10 +69,11 @@ class [[gnu::visibility("default")]] GriddedResponse {
double frequency, size_t field_id) = 0;
/**
* @brief Calculate baseline-integrated (undersampled) beam for a single time
* step. This function makes use of @ref MakeIntegratedSnapshot(). Subclasses
* may override MakeIntegratedSnapshot() to implement a more efficient
* version.
* @brief Calculate a grid of baseline-integrated (undersampled) linear or
* squared beam correction Mueller matrix for a single
* time step. This function makes use of @ref MakeIntegratedSnapshot().
* Subclasses may override MakeIntegratedSnapshot() to implement a more
* efficient version.
*
* @param beam_mode Selects beam mode (element, array factor or full)
* @param destination Buffer for storing the result, should have size width *
......@@ -84,18 +85,19 @@ class [[gnu::visibility("default")]] GriddedResponse {
* @param baseline_weights Baseline weights, size should equal
* Telescope::GetNrStations() * (Telescope::GetNrStations() + 1)/2
*/
virtual void IntegratedResponse(BeamMode beam_mode, float* destination,
double time, double frequency,
size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights);
virtual void IntegratedCorrection(
BeamMode beam_mode, float* destination, double time, double frequency,
size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights, bool square_mueller);
/**
* @brief Calculate baseline-integrated beam over multiple time steps.
* @brief Calculate a grid of baseline-integrated linear or squared beam
* Mueller matrix correction over multiple time steps.
* This function makes use of @ref MakeIntegratedSnapshot(). Subclasses
* may override MakeIntegratedSnapshot() to implement a more efficient
* version. This function stores all Mueller matrices for the full size
* image in memory. If this is undesirable, the functions
* @ref UndersampledIntegratedResponse() and @ref UpsampleResponse() may
* @ref UndersampledIntegratedCorrection() and @ref UpsampleCorrection() may
* be used.
*
* @param beam_mode Selects beam mode (element, array factor or full)
......@@ -110,35 +112,36 @@ class [[gnu::visibility("default")]] GriddedResponse {
* (Telescope::GetNrStations() * (Telescope::GetNrStations() + 1)/2) *
* time_array.size()
*/
virtual void IntegratedResponse(
virtual void IntegratedCorrection(
BeamMode beam_mode, float* destination,
const std::vector<double>& time_array, double frequency, size_t field_id,
size_t undersampling_factor, const std::vector<double>& baseline_weights);
size_t undersampling_factor, const std::vector<double>& baseline_weights,
bool square_mueller);
/**
* Same as @ref IntegratedResponse(), but without performing the upsampling.
* Same as @ref IntegratedCorrection(), but without performing the upsampling.
* This function therefore returns the undersampled data.
* Function @ref UpsampleResponse() can be used to upsample the returned
* Function @ref UpsampleCorrection() can be used to upsample the returned
* data. This route is useful for big images, for which it is undesirable to
* hold all 16 elements of the Mueller matrix in memory at the same time.
* @returns The undersampled data, input to @ref UpsampleResponse().
* @returns The undersampled data, input to @ref UpsampleCorrection().
*/
virtual std::vector<aocommon::HMC4x4> UndersampledIntegratedResponse(
virtual std::vector<aocommon::HMC4x4> UndersampledIntegratedCorrection(
BeamMode beam_mode, const std::vector<double>& time_array,
double frequency, size_t field_id, size_t undersampling_factor,
const std::vector<double>& baseline_weights);
const std::vector<double>& baseline_weights, bool square_mueller);
/**
* Upsample a single element from an undersampled response.
* @param destination Result buffer for one Mueller matrix element of size
* Width() x Height().
* Upsample a grid with a single element from an undersampled response.
* @param destination Result buffer for one Mueller matrix element for a
* grid of size Width() x Height().
* @param element_index Value from 0 to 15 indicating which Mueller matrix to
* upsample.
* @param undersampled_beam The previously calculated undersampled response
* @param undersampling_factor Undersampling factor that was used in the
* response calculation call.
*/
static void UpsampleResponse(
static void UpsampleCorrection(
float* destination, size_t element_index, size_t width, size_t height,
const std::vector<aocommon::HMC4x4>& undersampled_beam,
size_t undersampling_factor);
......@@ -208,15 +211,17 @@ class [[gnu::visibility("default")]] GriddedResponse {
private:
/**
* @brief Calculate a baseline-integrated snapshot.
* By default, this function will request the response for all antennas and
* perform a weighted average. (Partly) homogeneous arrays may implement a
* faster implementation by overriding this method.
* @brief Calculate a Mueller matrix grid with a baseline-integrated snapshot
* of beam correction.
* By default, this function will request the response for all antennas, form
* the Mueller matrices, square these and perform a weighted average. (Partly)
* homogeneous arrays may implement a faster implementation by overriding this
* method.
*/
virtual void MakeIntegratedSnapshot(
virtual void MakeIntegratedCorrectionSnapshot(
BeamMode beam_mode, std::vector<aocommon::HMC4x4> & matrices, double time,
double frequency, size_t field_id,
const double* baseline_weights_interval);
const double* baseline_weights_interval, bool square_mueller);
};
} // namespace griddedresponse
} // namespace everybeam
......
......@@ -60,11 +60,10 @@ void MWAGrid::Response(BeamMode /* beam_mode */, std::complex<float>* buffer,
}
}
void MWAGrid::MakeIntegratedSnapshot(BeamMode beam_mode,
std::vector<aocommon::HMC4x4>& matrices,
double time, double frequency,
size_t field_id,
const double* baseline_weights_interval) {
void MWAGrid::MakeIntegratedCorrectionSnapshot(
BeamMode beam_mode, std::vector<aocommon::HMC4x4>& matrices, double time,
double frequency, size_t field_id, const double* baseline_weights_interval,
bool square_mueller) {
const size_t n_stations = telescope_->GetNrStations();
aocommon::UVector<std::complex<float>> buffer_undersampled(
GetStationBufferSize(n_stations));
......@@ -86,9 +85,16 @@ void MWAGrid::MakeIntegratedSnapshot(BeamMode beam_mode,
// Mueller matrix constant for all baselines, so just compute once for
// each individual pixel
matrices[y * width_ + x] =
aocommon::HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A) *
snapshot_weight;
if (square_mueller) {
matrices[y * width_ + x] =
aocommon::HMC4x4::KroneckerProduct(A.HermTranspose().Transpose(), A)
.Square() *
snapshot_weight;
} else {
matrices[y * width_ + x] = aocommon::HMC4x4::KroneckerProduct(
A.HermTranspose().Transpose(), A) *
snapshot_weight;
}
}
}
}
......
......@@ -32,10 +32,10 @@ class [[gnu::visibility("default")]] MWAGrid final : public GriddedResponse {
std::unique_ptr<everybeam::mwabeam::TileBeam2016> tile_beam_;
// Override MakeIntegrated snapshot for efficiency
void MakeIntegratedSnapshot(BeamMode beam_mode,
std::vector<aocommon::HMC4x4> & matrices,
double time, double frequency, size_t field_id,
const double* baseline_weights_interval) override;
void MakeIntegratedCorrectionSnapshot(
BeamMode beam_mode, std::vector<aocommon::HMC4x4> & matrices, double time,
double frequency, size_t field_id,
const double* baseline_weights_interval, bool square_mueller) override;
};
} // namespace griddedresponse
} // namespace everybeam
......
......@@ -89,9 +89,9 @@ void CheckUndersampledIntegratedResponse(size_t undersampling_factor) {
1.0);
const std::vector<aocommon::HMC4x4> result =
gmrt.grid_response->UndersampledIntegratedResponse(
gmrt.grid_response->UndersampledIntegratedCorrection(
everybeam::CorrectionMode::kFull, time_array, kFrequency, 0,
undersampling_factor, baseline_weights);
undersampling_factor, baseline_weights, false);
const size_t undersampled_width = kWidth / undersampling_factor;
const size_t undersampled_height = kHeight / undersampling_factor;
BOOST_REQUIRE_EQUAL(result.size(), undersampled_width * undersampled_height);
......@@ -110,8 +110,9 @@ void CheckUndersampledIntegratedResponse(size_t undersampling_factor) {
result[undersampled_width / 2].Data(element_index);
std::vector<float> buffer(kWidth * kHeight);
gmrt.grid_response->UpsampleResponse(buffer.data(), element_index, kWidth,
kHeight, result, undersampling_factor);
gmrt.grid_response->UpsampleCorrection(buffer.data(), element_index, kWidth,
kHeight, result,
undersampling_factor);
const float upsampled_centre_pixel =
buffer[kWidth * (kHeight / 2) + kWidth / 2];
......
......@@ -59,24 +59,27 @@ BOOST_AUTO_TEST_CASE(undersampled) {
const std::vector<double> baseline_weights(
telescope->GetNrStations() * (telescope->GetNrStations() + 1) / 2, 1.0);
// Use the simple interface to calculate the reference values
grid_response->IntegratedResponse(everybeam::BeamMode::kFull,
reference_data.data(), times, frequency, 0,
sampling_factor, baseline_weights);
// Now recalculate in separate steps
const std::vector<aocommon::HMC4x4> undersampled =
grid_response->UndersampledIntegratedResponse(
everybeam::BeamMode::kFull, times, frequency, 0, sampling_factor,
baseline_weights);
for (size_t element_index = 0; element_index != 16; ++element_index) {
std::vector<float> result(width * height);
GriddedResponse::UpsampleResponse(result.data(), element_index, width,
height, undersampled, sampling_factor);
for (size_t i = 0; i != result.size(); ++i) {
BOOST_CHECK_EQUAL(result[i],
reference_data[i + element_index * width * height]);
for (bool square_mueller : {false, true}) {
// Use the simple interface to calculate the reference values
grid_response->IntegratedCorrection(
everybeam::BeamMode::kFull, reference_data.data(), times, frequency, 0,
sampling_factor, baseline_weights, square_mueller);
// Now recalculate in separate steps
const std::vector<aocommon::HMC4x4> undersampled =
grid_response->UndersampledIntegratedCorrection(
everybeam::BeamMode::kFull, times, frequency, 0, sampling_factor,
baseline_weights, square_mueller);
for (size_t element_index = 0; element_index != 16; ++element_index) {
std::vector<float> result(width * height);
GriddedResponse::UpsampleCorrection(result.data(), element_index, width,
height, undersampled,
sampling_factor);
for (size_t i = 0; i != result.size(); ++i) {
BOOST_CHECK_EQUAL(result[i],
reference_data[i + element_index * width * height]);
}
}
}
}
......
......@@ -400,9 +400,9 @@ BOOST_AUTO_TEST_CASE(integrated_beam) {
grid_response->GetIntegratedBufferSize());
std::vector<double> baseline_weights(
telescope->GetNrStations() * (telescope->GetNrStations() + 1) / 2, 1.);
grid_response->IntegratedResponse(BeamMode::kFull,
antenna_buffer_integrated.data(), time,
frequency, 0, 2, baseline_weights);
grid_response->IntegratedCorrection(BeamMode::kFull,
antenna_buffer_integrated.data(), time,
frequency, 0, 2, baseline_weights, false);
// Just check whether some (rather arbitrary) numbers are reproduced
BOOST_CHECK_LT(std::abs(antenna_buffer_integrated[10] - 0.0309436), 1e-6);
......@@ -415,9 +415,9 @@ BOOST_AUTO_TEST_CASE(integrated_beam) {
std::vector<double> tarray = {time, time};
baseline_weights.resize(baseline_weights.size() * tarray.size());
std::fill(baseline_weights.begin(), baseline_weights.end(), 1.);
grid_response->IntegratedResponse(BeamMode::kFull,
antenna_buffer_integrated.data(), tarray,
frequency, 0, 2, baseline_weights);
grid_response->IntegratedCorrection(BeamMode::kFull,
antenna_buffer_integrated.data(), tarray,
frequency, 0, 2, baseline_weights, false);
BOOST_CHECK_LT(std::abs(antenna_buffer_integrated[10] - 0.0309436), 1e-6);
BOOST_CHECK_LT(std::abs(antenna_buffer_integrated[20] - 0.156267), 1e-6);
......@@ -461,9 +461,9 @@ BOOST_AUTO_TEST_CASE(integrated_beam) {
grid_response_pb->GetIntegratedBufferSize());
std::vector<double> baseline_weights_pb(
telescope->GetNrStations() * (telescope->GetNrStations() + 1) / 2, 1.);
grid_response_pb->IntegratedResponse(BeamMode::kFull,
antenna_buffer_pb.data(), time,
frequency, 0, 8, baseline_weights_pb);
grid_response_pb->IntegratedCorrection(
BeamMode::kFull, antenna_buffer_pb.data(), time, frequency, 0, 8,
baseline_weights_pb, false);
// Check diagonal and off-diagonal term in component 0 and 5 of HMC4x4
// representation of Mueller matrix
std::size_t offset_01616 = 16 * width_pb + 16,
......
Subproject commit ea534da2fda4ec24329c5ca2bbee344f42806546
Subproject commit 9485d91d660513519416d63ad2ce10e8a49c557d
......@@ -15,7 +15,7 @@ build-backend = "scikit_build_core.build"
[project]
name = "everybeam"
version = "0.5.8" # Keep in sync with top-level `CMakeLists.txt` file
version = "0.6.0" # Keep in sync with top-level `CMakeLists.txt` file
description = "EveryBeam"
readme = {file = "README.md", content-type = "text/markdown"}
requires-python = ">=3.7"
......
......@@ -297,13 +297,13 @@ void init_telescope(py::module& m) {
// Call overload depending on the size of the time_vec
if (time_vec.size() == 1) {
grid_response->IntegratedResponse(BeamMode::kFull, buffer.data(),
time_vec[0], freq, field_id,
undersampling_factor, bw_vec);
grid_response->IntegratedCorrection(
BeamMode::kFull, buffer.data(), time_vec[0], freq, field_id,
undersampling_factor, bw_vec, false);
} else {
grid_response->IntegratedResponse(BeamMode::kFull, buffer.data(),
time_vec, freq, field_id,
undersampling_factor, bw_vec);
grid_response->IntegratedCorrection(
BeamMode::kFull, buffer.data(), time_vec, freq, field_id,
undersampling_factor, bw_vec, false);
}
// The Hermitian matrices need to be stored as a Matrix4x4
const size_t npixels =
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment