diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index 27847c8f4643b5990fdbef202df6118d1afb67da..165dcb7b30321c20fd2185cadbd8c2d029ef5226 100644
--- a/benchmark/CMakeLists.txt
+++ b/benchmark/CMakeLists.txt
@@ -14,7 +14,8 @@ FetchContent_MakeAvailable(googlebenchmark)
set(BENCHMARKS simulator)
-add_executable(microbenchmarks main.cpp predict.cpp directions.cpp spectrum.cpp)
+add_executable(microbenchmarks main.cpp predict.cpp directions.cpp spectrum.cpp
+ spectrum_cross.cpp)
target_compile_options(
microbenchmarks PRIVATE $<$<CONFIG:Debug>:-g> $<$<CONFIG:Release>:-O3>
diff --git a/benchmark/spectrum.cpp b/benchmark/spectrum.cpp
index 295b2601c166bae5144563a2c9e7646e39d1d586..4e8007decf90f3a5ebe3384fafc0e5ce563b8dc9 100644
--- a/benchmark/spectrum.cpp
+++ b/benchmark/spectrum.cpp
@@ -61,7 +61,7 @@ protected:
BENCHMARK_DEFINE_F(SpectrumBenchmark, SpectrumBenchmark)
(benchmark::State &state) {
for (auto _ : state) {
- spectrum->evaluate(frequencies, *computed_spectrum);
+ spectrum->Evaluate(frequencies, *computed_spectrum);
}
}
diff --git a/benchmark/spectrum_cross.cpp b/benchmark/spectrum_cross.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..23c43f177353b3f2579c3f0e67db1d9919ea7e7a
--- /dev/null
+++ b/benchmark/spectrum_cross.cpp
@@ -0,0 +1,67 @@
+#include <Simulator.h>
+#include <benchmark/benchmark.h>
+#include <random>
+#include <xtensor/xtensor.hpp>
+
+#include <Direction.h>
+#include <PointSource.h>
+#include <Spectrum.h>
+const std::vector<int64_t> n_frequencies = {128, 256, 512};
+
+const double kHz = 1.e3;
+const double MHz = 1.e6;
+namespace {} // namespace
+
+class SpectrumCrossBenchmark : public benchmark::Fixture {
+public:
+ void SetUp(benchmark::State &state) override {
+ n_frequencies = state.range(0);
+ const bool is_log = state.range(1);
+ const bool compute_rotation = state.range(2);
+ frequencies.resize({n_frequencies});
+ for (int i = 0; i < n_frequencies; i++) {
+ frequencies(i) = i * 100. * kHz + MHz;
+ }
+
+ spectrum = std::make_unique<Spectrum>();
+ auto &spec = *spectrum;
+ spec.has_logarithmic_spectral_index = is_log;
+ spec.reference_flux = {1.0, 0.0, 0.0, 0.0};
+ spec.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
+ spec.reference_frequency = MHz;
+ spec.polarization_angle = 0.3;
+ spec.polarization_factor = 0.2;
+ spec.rotation_measure = 0.1;
+ spec.has_rotation_measure = compute_rotation;
+ computed_spectrum =
+ std::make_unique<xt::xtensor<std::complex<double>, 2>>();
+ computed_spectrum->resize({4, n_frequencies});
+ }
+ void TearDown(benchmark::State &) override {
+ spectrum.reset();
+ computed_spectrum.reset();
+ }
+
+protected:
+ std::unique_ptr<Spectrum> spectrum;
+ std::unique_ptr<xt::xtensor<std::complex<double>, 2>> computed_spectrum;
+ size_t n_frequencies;
+ xt::xtensor<double, 1> frequencies;
+};
+
+BENCHMARK_DEFINE_F(SpectrumCrossBenchmark, SpectrumBenchmark)
+(benchmark::State &state) {
+ for (auto _ : state) {
+ spectrum->EvaluateCrossCorrelations(frequencies, *computed_spectrum);
+ }
+}
+
+BENCHMARK_REGISTER_F(SpectrumCrossBenchmark, SpectrumBenchmark)
+ ->ArgsProduct({
+ n_frequencies,
+ {false, true},
+ {false, true},
+
+ })
+ ->ArgNames({"#frequencies", "is_log", "compute_rotation"})
+ ->Unit(benchmark::kNanosecond);
diff --git a/include/Spectrum.h b/include/Spectrum.h
index 94f2154f019d8f836dedaee610570549afe2052b..194bce45413786f20fc01cece8285ddb8fee60d7 100644
--- a/include/Spectrum.h
+++ b/include/Spectrum.h
@@ -13,9 +13,13 @@
#include <StokesVector.h>
#include <casacore/casa/BasicSL/Constants.h>
+#include <complex>
+#include <xtensor/xcomplex.hpp>
+#include <xtensor/xshape.hpp>
#include <xtensor/xtensor.hpp>
-
+#include <xtensor/xview.hpp>
using Stokes = dp3::base::Stokes;
+enum CrossCorrelation { XX = 0, XY = 1, YX = 2, YY = 3 };
class Spectrum {
public:
@@ -29,7 +33,7 @@ public:
bool has_rotation_measure;
bool has_logarithmic_spectral_index;
- void evaluate(const xt::xtensor<double, 1> &frequencies,
+ void Evaluate(const xt::xtensor<double, 1> &frequencies,
StokesVector &result) const {
result.I.resize(frequencies.size());
result.Q.resize(frequencies.size());
@@ -37,20 +41,26 @@ public:
result.V.resize(frequencies.size());
double I0 = reference_flux.I;
+ double V0 = reference_flux.V;
if (spectral_terms.size() > 0 && has_logarithmic_spectral_index) {
for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
- result.I[f_idx] = evaluate_log_spectrum(
- I0, frequencies[f_idx], reference_frequency, spectral_terms);
+ const double spectral_term = EvaluateLogSpectrum(
+ frequencies[f_idx], reference_frequency, spectral_terms);
+ result.I[f_idx] = I0 * spectral_term;
+ result.V[f_idx] = V0 * spectral_term;
}
} else if (spectral_terms.size()) {
for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
- result.I[f_idx] = evaluate_spectrum(
- I0, frequencies[f_idx], reference_frequency, spectral_terms);
+ const double spectral_term = EvaluateSpectrum(
+ frequencies[f_idx], reference_frequency, spectral_terms);
+ result.I[f_idx] = I0 + spectral_term;
+ result.V[f_idx] = V0 + spectral_term;
}
} else {
for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
result.I[f_idx] = I0;
+ result.V[f_idx] = V0;
}
}
@@ -67,10 +77,70 @@ public:
}
}
+ void EvaluateCrossCorrelations(
+ const xt::xtensor<double, 1> &frequencies,
+ xt::xtensor<std::complex<double>, 2> &result) const {
+ xt::xtensor<double, 1>::shape_type kSpectralCorrectionShape{
+ frequencies.size()};
+ xt::xtensor<double, 2>::shape_type kRotationCoefficientsShape{
+ 2, frequencies.size()};
+ xt::xtensor<double, 1> spectral_correction =
+ xt::zeros<double>(kSpectralCorrectionShape);
+ xt::xtensor<double, 2> rotation_coefficients =
+ xt::zeros<double>(kRotationCoefficientsShape);
+
+ result.resize({4, frequencies.size()});
+
+ if (spectral_terms.size() > 0 && has_logarithmic_spectral_index) {
+ for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+ spectral_correction[f_idx] = EvaluateLogSpectrum(
+ frequencies[f_idx], reference_frequency, spectral_terms);
+ }
+ } else if (spectral_terms.size() > 0) {
+ for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+ spectral_correction[f_idx] = EvaluateSpectrum(
+ frequencies[f_idx], reference_frequency, spectral_terms);
+ }
+ }
+
+ if (has_rotation_measure) {
+ for (size_t f_idx = 0; f_idx < frequencies.size(); f_idx++) {
+ const double kLambda = casacore::C::c / frequencies[f_idx];
+ const double kChi =
+ 2.0 * (polarization_angle + (rotation_measure * kLambda * kLambda));
+ const double stokesQU = polarization_factor;
+ rotation_coefficients(0, f_idx) = stokesQU * cos(kChi);
+ rotation_coefficients(1, f_idx) = stokesQU * sin(kChi);
+ }
+ }
+
+ xt::xtensor<double, 1> I;
+ xt::xtensor<double, 1> V;
+
+ if (has_logarithmic_spectral_index) {
+ I = spectral_correction * reference_flux.I;
+ V = spectral_correction * reference_flux.V;
+
+ } else {
+ I = spectral_correction + reference_flux.I;
+ V = spectral_correction + reference_flux.V;
+ }
+ 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;
+ }
+
private:
inline double
- evaluate_polynomial(double x,
- const xt::xtensor<double, 1> ¶meters) const {
+ EvaluatePolynomial(double x, const xt::xtensor<double, 1> ¶meters) const {
double partial = 0.0;
for (auto it = spectral_terms.rbegin(), end = spectral_terms.rend();
it != end; ++it) {
@@ -79,26 +149,24 @@ private:
return partial;
}
inline double
- evaluate_log_polynomial(double x,
- const xt::xtensor<double, 1> ¶meters) const {
- double base = log10(x);
- return evaluate_polynomial(base, parameters);
+ EvaluateLogPolynomial(double x,
+ const xt::xtensor<double, 1> ¶meters) const {
+ const double base = log10(x);
+ return EvaluatePolynomial(base, parameters);
}
inline double
- evaluate_log_spectrum(double I0, const double frequency,
- const double reference_frequency,
- const xt::xtensor<double, 1> ¶meters) const {
+ EvaluateLogSpectrum(const double frequency, const double reference_frequency,
+ const xt::xtensor<double, 1> ¶meters) const {
double x = frequency / reference_frequency;
- return I0 * pow(10, x * evaluate_log_polynomial(x, parameters));
+ return pow(10, x * EvaluateLogPolynomial(x, parameters));
}
inline double
- evaluate_spectrum(double I0, const double frequency,
- const double reference_frequency,
- const xt::xtensor<double, 1> ¶meters) const {
- double x = frequency / reference_frequency - 1.0;
+ EvaluateSpectrum(const double frequency, const double reference_frequency,
+ const xt::xtensor<double, 1> ¶meters) const {
+ const double x = frequency / reference_frequency - 1.0;
- return I0 + evaluate_polynomial(x, parameters) * x;
+ return EvaluatePolynomial(x, parameters) * x;
}
};
diff --git a/src/PointSource.cpp b/src/PointSource.cpp
index 4fab35065a00557e6819a670571b6b0c296702ce..c024db982c0d4f1f2f087322ae1db116450ea73a 100644
--- a/src/PointSource.cpp
+++ b/src/PointSource.cpp
@@ -65,6 +65,7 @@ Stokes PointSource::stokes(double freq) const {
// Compute I * (v / v0) ^ exponent, where I is the value of Stokes
// I at the reference frequency.
stokes.I *= pow(10., base * exponent);
+ stokes.V *= pow(10., base * exponent);
} else {
double x = freq / itsRefFreq - 1.0;
typedef std::vector<double>::const_reverse_iterator iterator_type;
@@ -75,6 +76,7 @@ Stokes PointSource::stokes(double freq) const {
val = val * x + *it;
}
stokes.I += val * x;
+ stokes.V += val * x;
}
}
diff --git a/tests/test_spectrum.cpp b/tests/test_spectrum.cpp
index 74bdaeb95946af02140ec90317d7a021ae3a4362..ad1013f5aaf1fed1e0bcfd917cde8fbd58eff16e 100644
--- a/tests/test_spectrum.cpp
+++ b/tests/test_spectrum.cpp
@@ -45,7 +45,7 @@ BOOST_FIXTURE_TEST_CASE(normal_no_rotation, SpectrumFixture) {
stokes_nu2 = point_source.stokes(frequencies[2]);
StokesVector result;
- spectrum.evaluate(frequencies, result);
+ spectrum.Evaluate(frequencies, result);
BOOST_CHECK_CLOSE(result.I[0], stokes_nu0.I, 1.e-6);
BOOST_CHECK_CLOSE(result.I[1], stokes_nu1.I, 1.e-6);
@@ -82,7 +82,7 @@ BOOST_FIXTURE_TEST_CASE(normal_with_rotation, SpectrumFixture) {
stokes_nu2 = point_source.stokes(frequencies[2]);
StokesVector result;
- spectrum.evaluate(frequencies, result);
+ spectrum.Evaluate(frequencies, result);
BOOST_CHECK_CLOSE(result.I[0], stokes_nu0.I, 1.e-6);
BOOST_CHECK_CLOSE(result.I[1], stokes_nu1.I, 1.e-6);
@@ -117,7 +117,7 @@ BOOST_FIXTURE_TEST_CASE(normal_log_no_rotation, SpectrumFixture) {
stokes_nu2 = point_source.stokes(frequencies[2]);
StokesVector result;
- spectrum.evaluate(frequencies, result);
+ spectrum.Evaluate(frequencies, result);
BOOST_CHECK_CLOSE_FRACTION(result.I[0], stokes_nu0.I, 1.e-3);
BOOST_CHECK_CLOSE_FRACTION(result.I[1], stokes_nu1.I, 1.e-3);
@@ -158,7 +158,7 @@ BOOST_FIXTURE_TEST_CASE(normal_log_with_rotation, SpectrumFixture) {
stokes_nu2 = point_source.stokes(frequencies[2]);
StokesVector result;
- spectrum.evaluate(frequencies, result);
+ spectrum.Evaluate(frequencies, result);
BOOST_CHECK_CLOSE_FRACTION(result.I[0], stokes_nu0.I, 1.e-3);
BOOST_CHECK_CLOSE_FRACTION(result.I[1], stokes_nu1.I, 1.e-3);
@@ -174,4 +174,74 @@ BOOST_FIXTURE_TEST_CASE(normal_log_with_rotation, SpectrumFixture) {
BOOST_CHECK_CLOSE_FRACTION(result.V[2], stokes_nu2.V, 1.e-3);
}
+BOOST_FIXTURE_TEST_CASE(normal_log_with_rotation_cross, SpectrumFixture) {
+ spectrum.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
+ spectrum.has_logarithmic_spectral_index = true;
+ spectrum.has_rotation_measure = true;
+ spectrum.has_logarithmic_spectral_index = true;
+
+ dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
+ spectrum.reference_flux);
+
+ point_source.setSpectralTerms(
+ spectrum.reference_frequency, spectrum.has_logarithmic_spectral_index,
+ spectrum.spectral_terms.begin(), spectrum.spectral_terms.end());
+
+ point_source.setRotationMeasure(spectrum.polarization_factor,
+ spectrum.polarization_angle,
+ spectrum.rotation_measure);
+ Stokes stokes_nu0(0, 0, 0, 0);
+ Stokes stokes_nu1(0, 0, 0, 0);
+ Stokes stokes_nu2(0, 0, 0, 0);
+
+ stokes_nu0 = point_source.stokes(frequencies[0]);
+ stokes_nu1 = point_source.stokes(frequencies[1]);
+ stokes_nu2 = point_source.stokes(frequencies[2]);
+
+ xt::xtensor<std::complex<double>, 2> result_complex;
+ spectrum.EvaluateCrossCorrelations(frequencies, result_complex);
+
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0).real(),
+ stokes_nu0.I + stokes_nu0.Q, 1.e-3);
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 1).real(),
+ stokes_nu1.I + stokes_nu1.Q, 1.e-3);
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 2).real(),
+ stokes_nu2.I + stokes_nu2.Q, 1.e-3);
+}
+
+BOOST_FIXTURE_TEST_CASE(normal_with_rotation_cross, SpectrumFixture) {
+ spectrum.spectral_terms = {1.0e-6, 2.0e-6, 3.0e-6};
+ spectrum.has_logarithmic_spectral_index = true;
+ spectrum.has_rotation_measure = true;
+ spectrum.has_logarithmic_spectral_index = false;
+
+ dp3::base::PointSource point_source(dp3::base::Direction(0.0, 0.0),
+ spectrum.reference_flux);
+
+ point_source.setSpectralTerms(
+ spectrum.reference_frequency, spectrum.has_logarithmic_spectral_index,
+ spectrum.spectral_terms.begin(), spectrum.spectral_terms.end());
+
+ point_source.setRotationMeasure(spectrum.polarization_factor,
+ spectrum.polarization_angle,
+ spectrum.rotation_measure);
+ Stokes stokes_nu0(0, 0, 0, 0);
+ Stokes stokes_nu1(0, 0, 0, 0);
+ Stokes stokes_nu2(0, 0, 0, 0);
+
+ stokes_nu0 = point_source.stokes(frequencies[0]);
+ stokes_nu1 = point_source.stokes(frequencies[1]);
+ stokes_nu2 = point_source.stokes(frequencies[2]);
+
+ xt::xtensor<std::complex<double>, 2> result_complex;
+ spectrum.EvaluateCrossCorrelations(frequencies, result_complex);
+
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 0).real(),
+ stokes_nu0.I + stokes_nu0.Q, 1.e-3);
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 1).real(),
+ stokes_nu1.I + stokes_nu1.Q, 1.e-3);
+ BOOST_CHECK_CLOSE_FRACTION(result_complex(0, 2).real(),
+ stokes_nu2.I + stokes_nu2.Q, 1.e-3);
+}
+
BOOST_AUTO_TEST_SUITE_END()
\ No newline at end of file