diff --git a/benchmark/predict.cpp b/benchmark/predict.cpp
index 978bb579fca81bfc4896bd789a6d3f8cfa59b1d1..2b175ef17d23dd2b015f75d08c1f15b9588869df 100644
--- a/benchmark/predict.cpp
+++ b/benchmark/predict.cpp
@@ -62,7 +62,6 @@ BENCHMARK_REGISTER_F(PredictBenchmark, PointSource)
nstations,
nsources,
{false, true},
- {false, true},
})
->ArgNames({"stat", "chan", "fullstokes", "freqsmear"})
->Unit(benchmark::kMillisecond);
diff --git a/benchmark/spectrum_cross.cpp b/benchmark/spectrum_cross.cpp
index 13def0baaa8bd11c2244f643fc8eae57bbff9813..92cc24a547b4cc8aca34e90c658534058111980b 100644
--- a/benchmark/spectrum_cross.cpp
+++ b/benchmark/spectrum_cross.cpp
@@ -29,9 +29,8 @@ public:
spec.SetReferenceFlux({1.0, 0.0, 0.0, 0.0});
spec.SetPolarization(0.3, 0.2);
- computed_spectrum =
- std::make_unique<xt::xtensor<std::complex<double>, 2>>();
- computed_spectrum->resize({4, n_frequencies});
+ computed_spectrum = std::make_unique<xt::xtensor<double, 3>>();
+ computed_spectrum->resize({2, 4, n_frequencies});
}
void TearDown(benchmark::State &) override {
spectrum.reset();
@@ -40,7 +39,7 @@ public:
protected:
std::unique_ptr<Spectrum> spectrum;
- std::unique_ptr<xt::xtensor<std::complex<double>, 2>> computed_spectrum;
+ std::unique_ptr<xt::xtensor<double, 3>> computed_spectrum;
size_t n_frequencies;
xt::xtensor<double, 1> frequencies;
};
diff --git a/include/PointSource.h b/include/PointSource.h
index bc09b0b8e120191c1b4f693d8fac812e11ce9e6d..fb2bf3ec4f5afc1158721dfcff88b1825b5b7eba 100644
--- a/include/PointSource.h
+++ b/include/PointSource.h
@@ -31,7 +31,7 @@ public:
void SetDirection(const Direction &position);
void ComputeSpectrum(const xt::xtensor<double, 1> &frequencies,
- xt::xtensor<std::complex<double>, 2> &result) const;
+ xt::xtensor<double, 3> &result) const;
const Spectrum &GetSpectrum() const { return spectrum_; }
Stokes GetStokes(double freq) const;
diff --git a/include/PredictPlan.h b/include/PredictPlan.h
index d4ac95a3359b57fbe94a68c0bf2cb6eda1cc364b..7b103dbcda5df8b14fa5801c9f3ee06579370cd1 100644
--- a/include/PredictPlan.h
+++ b/include/PredictPlan.h
@@ -43,13 +43,13 @@ struct PredictPlan {
// standalone PredictPlan object
virtual void Precompute(const PointSourceCollection &) {}
- virtual void Compute(const PointSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3,
- xt::layout_type::row_major> &buffer) {}
+ virtual void
+ Compute(const PointSourceCollection &sources,
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major> &buffer) {}
- virtual void Compute(const GaussianSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3,
- xt::layout_type::row_major> &buffer) {}
+ virtual void
+ Compute(const GaussianSourceCollection &sources,
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major> &buffer) {}
};
#endif // PREDICT_PLAN_H_
diff --git a/include/PredictPlanExec.h b/include/PredictPlanExec.h
index f088bbc5614596183bd9e848f3ad3442dc5036dc..a9106ad2948e0cdb919fabee2fa955ed387b9d76 100644
--- a/include/PredictPlanExec.h
+++ b/include/PredictPlanExec.h
@@ -27,19 +27,17 @@ public:
virtual void Precompute(const PointSourceCollection &sources) final override;
- virtual void
- Compute(const PointSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3, xt::layout_type::row_major>
- &buffer) final override;
+ virtual void Compute(const PointSourceCollection &sources,
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major>
+ &buffer) final override;
- virtual void
- Compute(const GaussianSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3, xt::layout_type::row_major>
- &buffer) final override;
+ virtual void Compute(const GaussianSourceCollection &sources,
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major>
+ &buffer) final override;
void ComputeWithTarget(
const PointSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3, xt::layout_type::row_major> &buffer,
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major> &buffer,
const computation_strategy strat);
void ComputeStationPhases();
@@ -52,6 +50,10 @@ public:
return shift_data;
}
+ const xt::xtensor<double, 3> &GetShiftDataComplexSeperated() const {
+ return shift_data_complex_seperated;
+ }
+
void FillEulerMatrix(xt::xtensor<double, 2> &mat, const double ra,
const double dec);
@@ -71,6 +73,7 @@ public:
xt::xtensor<double, 2> lmn;
xt::xtensor<double, 1> station_phases;
xt::xtensor<std::complex<double>, 2> shift_data;
+ xt::xtensor<double, 3> shift_data_complex_seperated;
xt::xtensor<double, 1> smear_terms;
private:
diff --git a/include/Spectrum.h b/include/Spectrum.h
index a9d7306c22413d5a67a6caee495b610d5b453d0e..7cf2800eec8f91769bbeed0f0b7483e5b8a6923b 100644
--- a/include/Spectrum.h
+++ b/include/Spectrum.h
@@ -117,7 +117,7 @@ public:
const Stokes &GetReferenceFlux() const { return reference_flux_; }
void EvaluateCrossCorrelations(const xt::xtensor<double, 1> &frequencies,
- xt::xtensor<std::complex<double>, 2> &result,
+ xt::xtensor<double, 3> &result,
const bool need_reshape = true) const {
xt::xtensor<double, 1>::shape_type kSpectralCorrectionShape{
frequencies.size()};
@@ -129,7 +129,7 @@ public:
xt::zeros<double>(kRotationCoefficientsShape);
if (need_reshape) {
- result.resize({4, frequencies.size()});
+ result.resize({2, 4, frequencies.size()});
}
if (spectral_terms_.size() > 0 && has_logarithmic_spectral_index_) {
@@ -168,14 +168,14 @@ public:
auto Q = I * xt::view(rotation_coefficients, 0, xt::all());
auto U = I * xt::view(rotation_coefficients, 1, xt::all());
- xt::real(xt::view(result, CrossCorrelation::XX, xt::all())) = I + Q;
- xt::imag(xt::view(result, CrossCorrelation::XX, xt::all())) = 0.0;
- xt::real(xt::view(result, CrossCorrelation::XY, xt::all())) = U;
- xt::imag(xt::view(result, CrossCorrelation::XY, xt::all())) = V;
- xt::real(xt::view(result, CrossCorrelation::YX, xt::all())) = U;
- xt::imag(xt::view(result, CrossCorrelation::YX, xt::all())) = -V;
- xt::real(xt::view(result, CrossCorrelation::YY, xt::all())) = I - Q;
- xt::imag(xt::view(result, CrossCorrelation::YY, xt::all())) = 0.0;
+ xt::view(result, 0, CrossCorrelation::XX, xt::all()) = I + Q;
+ xt::view(result, 1, CrossCorrelation::XX, xt::all()) = 0.0;
+ xt::view(result, 0, CrossCorrelation::XY, xt::all()) = U;
+ xt::view(result, 1, CrossCorrelation::XY, xt::all()) = V;
+ xt::view(result, 0, CrossCorrelation::YX, xt::all()) = U;
+ xt::view(result, 1, CrossCorrelation::YX, xt::all()) = -V;
+ xt::view(result, 0, CrossCorrelation::YY, xt::all()) = I - Q;
+ xt::view(result, 1, CrossCorrelation::YY, xt::all()) = 0.0;
}
private:
diff --git a/src/PointSource.cpp b/src/PointSource.cpp
index 4b1be08f09026ad15287420e4edf913de71ed20b..7b71cf1de851da8fbe15d932ac6c29fa999731ab 100644
--- a/src/PointSource.cpp
+++ b/src/PointSource.cpp
@@ -28,9 +28,8 @@ void PointSource::SetDirection(const Direction &direction) {
direction_ = direction;
}
-void PointSource::ComputeSpectrum(
- const xt::xtensor<double, 1> &frequencies,
- xt::xtensor<std::complex<double>, 2> &result) const {
+void PointSource::ComputeSpectrum(const xt::xtensor<double, 1> &frequencies,
+ xt::xtensor<double, 3> &result) const {
spectrum_.EvaluateCrossCorrelations(frequencies, result);
}
diff --git a/src/PredictPlanExec.cpp b/src/PredictPlanExec.cpp
index 9126c4ef16800c43f57b8d18bfa3e0587831994a..85d15dfc2f48339230b865444ce628c92b6510ca 100644
--- a/src/PredictPlanExec.cpp
+++ b/src/PredictPlanExec.cpp
@@ -6,10 +6,12 @@
#include "PredictPlan.h"
#include "Spectrum.h"
+#include <cstdint>
#include <immintrin.h>
#include <iomanip>
#include <iostream>
+#include <thread>
#include <xtensor-blas/xlinalg.hpp>
#include <xtensor/xlayout.hpp>
#include <xtensor/xtensor.hpp>
@@ -19,6 +21,7 @@ void PredictPlanExec::Initialize() {
// Initialize the station phases and shifts
station_phases.resize({nstations});
shift_data.resize({nstations, nchannels});
+ shift_data_complex_seperated.resize({2, nstations, nchannels});
lmn.resize({nchannels, 3});
uvw.resize({nstations, 3});
frequencies.resize({nchannels});
@@ -83,6 +86,8 @@ void PredictPlanExec::ComputeStationPhases() {
// Assign to shifts as std::complex
for (size_t ch = 0; ch < nchannels; ++ch) {
shift_data(st, ch) = std::complex<double>(cos_phase(ch), sin_phase(ch));
+ shift_data_complex_seperated(0, st, ch) = cos_phase(ch);
+ shift_data_complex_seperated(1, st, ch) = sin_phase(ch);
}
}
}
@@ -149,8 +154,8 @@ void PredictPlanExec::ProcessPolarizationComponentSingle(
const double y_p = (GetShiftData()(p, ch).imag());
const double x_q = (GetShiftData()(q, ch).real());
const double y_q = (GetShiftData()(q, ch).imag());
- const double x_c = stoke_spectrum(ch).real();
- const double y_c = stoke_spectrum(ch).imag();
+ const double x_c = stoke_spectrum(0, ch);
+ const double y_c = stoke_spectrum(1, ch);
// Compute baseline phase shift.
// Compute visibilities.
@@ -200,44 +205,48 @@ void PredictPlanExec::ProcessPolarizationComponentMultiSimd(E &buffer_expr,
const S &spectrum) {
// Use preallocated storage (outside for loops)
- const size_t simd_channels = (nchannels / 4) * 4;
- const size_t total_vector_size = simd_channels * 4;
+ const uint_fast32_t simd_channels = (nchannels / 4) * 4;
+ const uint_fast32_t total_vector_size = simd_channels * 4;
- // std::cout << "simd_channels: " << simd_channels << std::endl;
- // std::cout << "remainder: " << remainder << std::endl;
-
- for (size_t bl = 0; bl < baselines.size(); ++bl) {
+ for (uint_fast32_t bl = 0; bl < baselines.size(); ++bl) {
const size_t p = baselines[bl].first;
const size_t q = baselines[bl].second;
- const auto &shift_data = GetShiftData(); // Cache reference
+ const auto &shift_data = GetShiftDataComplexSeperated(); // Cache reference
- size_t ch = 0;
+ uint_fast32_t ch = 0;
for (; ch < simd_channels; ch += 4) {
// Load 4 values from shift_data for p and q
- __m256d x_p = _mm256_set_pd(
- shift_data(p, ch + 3).real(), shift_data(p, ch + 2).real(),
- shift_data(p, ch + 1).real(), shift_data(p, ch + 0).real());
+ // __m256d x_p = _mm256_set_pd(
+ // shift_data(p, ch + 3).real(), shift_data(p, ch + 2).real(),
+ // shift_data(p, ch + 1).real(), shift_data(p, ch + 0).real());
- __m256d y_p = _mm256_set_pd(
- shift_data(p, ch + 3).imag(), shift_data(p, ch + 2).imag(),
- shift_data(p, ch + 1).imag(), shift_data(p, ch + 0).imag());
+ // __m256d y_p = _mm256_set_pd(
+ // shift_data(p, ch + 3).imag(), shift_data(p, ch + 2).imag(),
+ // shift_data(p, ch + 1).imag(), shift_data(p, ch + 0).imag());
- __m256d x_q = _mm256_set_pd(
- shift_data(q, ch + 3).real(), shift_data(q, ch + 2).real(),
- shift_data(q, ch + 1).real(), shift_data(q, ch + 0).real());
+ // __m256d x_q = _mm256_set_pd(
+ // shift_data(q, ch + 3).real(), shift_data(q, ch + 2).real(),
+ // shift_data(q, ch + 1).real(), shift_data(q, ch + 0).real());
- __m256d y_q = _mm256_set_pd(
- shift_data(q, ch + 3).imag(), shift_data(q, ch + 2).imag(),
- shift_data(q, ch + 1).imag(), shift_data(q, ch + 0).imag());
+ // __m256d y_q = _mm256_set_pd(
+ // shift_data(q, ch + 3).imag(), shift_data(q, ch + 2).imag(),
+ // shift_data(q, ch + 1).imag(), shift_data(q, ch + 0).imag());
- __m256d x_c =
- _mm256_set_pd(spectrum(ch + 3).real(), spectrum(ch + 2).real(),
- spectrum(ch + 1).real(), spectrum(ch + 0).real());
+ // __m256d x_c =
+ // _mm256_set_pd(spectrum(ch + 3).real(), spectrum(ch + 2).real(),
+ // spectrum(ch + 1).real(), spectrum(ch + 0).real());
- __m256d y_c =
- _mm256_set_pd(spectrum(ch + 3).imag(), spectrum(ch + 2).imag(),
- spectrum(ch + 1).imag(), spectrum(ch + 0).imag());
+ // __m256d y_c =
+ // _mm256_set_pd(spectrum(ch + 3).imag(), spectrum(ch + 2).imag(),
+ // spectrum(ch + 1).imag(), spectrum(ch + 0).imag());
+
+ __m256d x_p = _mm256_loadu_pd(&shift_data(0, p, ch + 0));
+ __m256d y_p = _mm256_loadu_pd(&shift_data(1, p, ch + 0));
+ __m256d x_q = _mm256_loadu_pd(&shift_data(0, q, ch + 0));
+ __m256d y_q = _mm256_loadu_pd(&shift_data(1, q, ch + 0));
+ __m256d x_c = _mm256_loadu_pd(&spectrum(0, ch + 0));
+ __m256d y_c = _mm256_loadu_pd(&spectrum(1, ch + 0));
// Complex multiply (q_conj * p) * stokes
__m256d real_part =
@@ -250,29 +259,27 @@ void PredictPlanExec::ProcessPolarizationComponentMultiSimd(E &buffer_expr,
__m256d temp_imag = _mm256_add_pd(_mm256_mul_pd(real_part, y_c),
_mm256_mul_pd(imag_part, x_c));
+ if (correct_frequency_smearing) {
+ __m256d smear_terms_temp = _mm256_loadu_pd(&smear_terms[ch]);
+ temp_real = _mm256_mul_pd(temp_real, smear_terms_temp);
+ temp_imag = _mm256_mul_pd(temp_imag, smear_terms_temp);
+ }
+
alignas(32) double r[4], i[4];
_mm256_store_pd(r, temp_real);
_mm256_store_pd(i, temp_imag);
- if (correct_frequency_smearing) {
- for (int k = 0; k < 4; ++k) {
- const double smear = smear_terms[ch + k];
- buffer_expr(bl, ch + k) += dcomplex{r[k] * smear, i[k] * smear};
- }
- } else {
- for (int k = 0; k < 4; ++k) {
- buffer_expr(bl, ch + k) += dcomplex{r[k], i[k]};
- }
+ for (uint_fast8_t k = 0; k < 4; ++k) {
+ buffer_expr(bl, ch + k) += dcomplex{r[k], i[k]};
}
}
-
for (; ch < nchannels; ++ch) {
- const double x_p = (shift_data(p, ch).real());
- const double y_p = (shift_data(p, ch).imag());
- const double x_q = (shift_data(q, ch).real());
- const double y_q = (shift_data(q, ch).imag());
- const double x_c = spectrum(ch).real();
- const double y_c = spectrum(ch).imag();
+ const double x_p = (shift_data(0, p, ch));
+ const double y_p = (shift_data(1, p, ch));
+ const double x_q = (shift_data(0, q, ch));
+ const double y_q = (shift_data(1, q, ch));
+ const double x_c = spectrum(0, ch);
+ const double y_c = spectrum(1, ch);
// Compute baseline phase shift.
// Compute visibilities.
@@ -296,14 +303,13 @@ void PredictPlanExec::ProcessPolarizationComponentMultiSimd(E &buffer_expr,
void PredictPlanExec::Compute(
const PointSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3, xt::layout_type::row_major> &buffer) {
- xt::xtensor<std::complex<double>, 2> computed_spectrum(
- {4, frequencies.size()});
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major> &buffer) {
+ xt::xtensor<double, 3> computed_spectrum({2, 4, frequencies.size()});
- const auto spectrum_XX = xt::view(computed_spectrum, 0, xt::all());
- const auto spectrum_XY = xt::view(computed_spectrum, 1, xt::all());
- const auto spectrum_YX = xt::view(computed_spectrum, 2, xt::all());
- const auto spectrum_YY = xt::view(computed_spectrum, 3, xt::all());
+ const auto spectrum_XX = xt::view(computed_spectrum, xt::all(), 0, xt::all());
+ const auto spectrum_XY = xt::view(computed_spectrum, xt::all(), 1, xt::all());
+ const auto spectrum_YX = xt::view(computed_spectrum, xt::all(), 2, xt::all());
+ const auto spectrum_YY = xt::view(computed_spectrum, xt::all(), 3, xt::all());
for (size_t source_index = 0; source_index < sources.Size(); ++source_index) {
const Spectrum &spectrum = sources.spectrums[source_index];
@@ -330,15 +336,14 @@ void PredictPlanExec::Compute(
void PredictPlanExec::ComputeWithTarget(
const PointSourceCollection &sources,
- xt::xtensor<std::complex<double>, 3, xt::layout_type::row_major> &buffer,
+ xt::xtensor<dcomplex, 3, xt::layout_type::row_major> &buffer,
const computation_strategy strat) {
- xt::xtensor<std::complex<double>, 2> computed_spectrum(
- {4, frequencies.size()});
+ xt::xtensor<double, 3> computed_spectrum({2, 4, frequencies.size()});
- const auto spectrum_XX = xt::view(computed_spectrum, 0, xt::all());
- const auto spectrum_XY = xt::view(computed_spectrum, 1, xt::all());
- const auto spectrum_YX = xt::view(computed_spectrum, 2, xt::all());
- const auto spectrum_YY = xt::view(computed_spectrum, 3, xt::all());
+ const auto spectrum_XX = xt::view(computed_spectrum, xt::all(), 0, xt::all());
+ const auto spectrum_XY = xt::view(computed_spectrum, xt::all(), 1, xt::all());
+ const auto spectrum_YX = xt::view(computed_spectrum, xt::all(), 2, xt::all());
+ const auto spectrum_YY = xt::view(computed_spectrum, xt::all(), 3, xt::all());
for (size_t source_index = 0; source_index < sources.Size(); ++source_index) {
const Spectrum &spectrum = sources.spectrums[source_index];
@@ -368,7 +373,7 @@ void PredictPlanExec::Compute(
xt::xtensor<std::complex<double>, 3, xt::layout_type::row_major> &buffer) {
const std::array<size_t, 2> kStationUvwShape = {get_uvw().shape(1),
get_uvw().shape(0)};
- xt::xtensor<std::complex<double>, 2> computed_spectrum;
+ xt::xtensor<double, 3> computed_spectrum;
xt::xtensor<double, 2> xtStationUvw;
xt::xtensor<double, 2> uvwShifted;
@@ -468,10 +473,8 @@ void PredictPlanExec::Compute(
const double y_p = (GetShiftData()(p, ch).imag());
const double x_q = (GetShiftData()(q, ch).real());
const double y_q = (GetShiftData()(q, ch).imag());
- const double x_c =
- (computed_spectrum(StokesComponent::I, ch).real());
- const double y_c =
- (computed_spectrum(StokesComponent::I, ch).imag());
+ const double x_c = (computed_spectrum(0, StokesComponent::I, ch));
+ const double y_c = (computed_spectrum(1, StokesComponent::I, ch));
// Compute baseline phase shift.
// Compute visibilities.
@@ -492,22 +495,20 @@ void PredictPlanExec::Compute(
const double y_p = (GetShiftData()(p, ch).imag());
const double x_q = (GetShiftData()(q, ch).real());
const double y_q = (GetShiftData()(q, ch).imag());
- const double x_c =
- (computed_spectrum(CrossCorrelation::XX, ch).real());
- const double y_c =
- (computed_spectrum(CrossCorrelation::XX, ch).imag());
+ const double x_c = (computed_spectrum(0, CrossCorrelation::XX, ch));
+ const double y_c = (computed_spectrum(1, CrossCorrelation::XX, ch));
const double x_c2 =
- (computed_spectrum(CrossCorrelation::XY, ch).real());
+ (computed_spectrum(0, CrossCorrelation::XY, ch));
const double y_c2 =
- (computed_spectrum(CrossCorrelation::XY, ch).imag());
+ (computed_spectrum(1, CrossCorrelation::XY, ch));
const double x_c3 =
- (computed_spectrum(CrossCorrelation::YX, ch).real());
+ (computed_spectrum(0, CrossCorrelation::YX, ch));
const double y_c3 =
- (computed_spectrum(CrossCorrelation::YX, ch).imag());
+ (computed_spectrum(1, CrossCorrelation::YX, ch));
const double x_c4 =
- (computed_spectrum(CrossCorrelation::YY, ch).real());
+ (computed_spectrum(0, CrossCorrelation::YY, ch));
const double y_c4 =
- (computed_spectrum(CrossCorrelation::YY, ch).imag());
+ (computed_spectrum(1, CrossCorrelation::YY, ch));
// Compute baseline phase shift.
// Compute visibilities.
diff --git a/tests/test_simulator.cpp b/tests/test_simulator.cpp
index 2c6758ea2402c13bfcd10e1c0f11c7df030199ca..43d158525723ef59c34c9d880f9a4cbdf6b5410a 100644
--- a/tests/test_simulator.cpp
+++ b/tests/test_simulator.cpp
@@ -380,6 +380,7 @@ BOOST_AUTO_TEST_CASE(test_calc_phase) {
plan.frequencies = {1e8, 2e8};
PredictPlanExec plan_exec{plan};
+ plan_exec.correct_frequency_smearing = false;
plan_exec.lmn = {{0.1, 0.2, 0.3}};
plan_exec.ComputeStationPhases();
diff --git a/tests/test_spectrum.cpp b/tests/test_spectrum.cpp
index 10badfb3df0f39e90b90baf4ee09876255802d3a..6176cc8ca1c7b10f40c179a1198ec5ccf334684b 100644
--- a/tests/test_spectrum.cpp
+++ b/tests/test_spectrum.cpp
@@ -227,14 +227,14 @@ BOOST_FIXTURE_TEST_CASE(normal_log_with_rotation_cross, SpectrumFixture) {
stokes_nu1 = point_source.GetStokes(frequencies[1]);
stokes_nu2 = point_source.GetStokes(frequencies[2]);
- xt::xtensor<std::complex<double>, 2> result_complex;
+ xt::xtensor<double, 3> result_complex;
spectrum.EvaluateCrossCorrelations(frequencies, result_complex);
- BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0).real(),
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0, 0),
stokes_nu0.I + stokes_nu0.Q, 1.e-3);
- BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 1).real(),
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0, 1),
stokes_nu1.I + stokes_nu1.Q, 1.e-3);
- BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 2).real(),
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0, 2),
stokes_nu2.I + stokes_nu2.Q, 1.e-3);
}
@@ -261,14 +261,14 @@ BOOST_FIXTURE_TEST_CASE(normal_with_rotation_cross, SpectrumFixture) {
stokes_nu1 = point_source.GetStokes(frequencies[1]);
stokes_nu2 = point_source.GetStokes(frequencies[2]);
- xt::xtensor<std::complex<double>, 2> result_complex;
+ xt::xtensor<double, 3> result_complex;
spectrum.EvaluateCrossCorrelations(frequencies, result_complex);
- BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0).real(),
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0, 0),
stokes_nu0.I + stokes_nu0.Q, 1.e-3);
- BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 1).real(),
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0, 1),
stokes_nu1.I + stokes_nu1.Q, 1.e-3);
- BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 2).real(),
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0, 2),
stokes_nu2.I + stokes_nu2.Q, 1.e-3);
}