From 32d16e46449cfe26b6bfb89ffabe545fc49cf2fd Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Fri, 17 Jul 2020 12:21:00 +0200
Subject: [PATCH] fixed interpolation
---
Environment/ShowerAxis.cc | 6 ++++--
Environment/ShowerAxis.h | 20 +++++++++++++-------
Environment/testShowerAxis.cc | 10 +++++-----
3 files changed, 22 insertions(+), 14 deletions(-)
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
index 5b3117873..c13fdd297 100644
--- a/Environment/ShowerAxis.cc
+++ b/Environment/ShowerAxis.cc
@@ -34,15 +34,17 @@ GrammageType ShowerAxis::X(LengthType l) const {
assert(0 <= lambda && lambda <= 1.);
+ std::cout << l << ": " << lower << " " << lambda << " " << upper << std::endl;
+
// linear interpolation between X[lower] and X[upper]
- return X_[lower] * lambda + X_[upper] * (1 - lambda);
+ return X_[upper] * lambda + X_[lower] * (1 - lambda);
}
LengthType ShowerAxis::steplength() const { return steplength_; }
GrammageType ShowerAxis::maximumX() const { return *X_.rbegin(); }
-GrammageType ShowerAxis::minimumX() const { return *X_.cbegin(); }
+GrammageType ShowerAxis::minimumX() const { return GrammageType::zero(); }
GrammageType ShowerAxis::projectedX(geometry::Point const& p) const {
auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
diff --git a/Environment/ShowerAxis.h b/Environment/ShowerAxis.h
index 3e713dd46..413455a64 100644
--- a/Environment/ShowerAxis.h
+++ b/Environment/ShowerAxis.h
@@ -23,6 +23,8 @@
#include <stdexcept>
#include <vector>
+#include <iostream>
+
#include <boost/math/quadrature/gauss_kronrod.hpp>
namespace corsika::environment {
@@ -37,7 +39,7 @@ namespace corsika::environment {
, max_length_(length_.norm())
, steplength_(max_length_ / steps)
, axis_normalized_(length / max_length_)
- , X_(steps) {
+ , X_(steps + 1) {
auto const* const universe = env.GetUniverse().get();
auto rho = [pStart, length, universe](double x) {
@@ -48,19 +50,23 @@ namespace corsika::environment {
double error;
int k = 0;
- for (int i = 1; i < steps; ++i) {
- auto const x_prev = (i - 1) / (steps - 1.);
+ X_[0] = units::si::GrammageType::zero();
+ auto sum = units::si::GrammageType::zero();
+
+ for (int i = 1; i <= steps; ++i) {
+ auto const x_prev = (i - 1.) / steps;
auto const d_prev = max_length_ * x_prev;
- auto const x = i / (steps - 1.);
+ auto const x = double(i) / steps;
auto const r = boost::math::quadrature::gauss_kronrod<double, 15>::integrate(
rho, x_prev, x, 15, 1e-9, &error);
auto const result =
units::si::MassDensityType(phys::units::detail::magnitude_tag, r) *
max_length_;
- auto const resultTotal = result + X_[i - 1];
- X_[i] = resultTotal;
- for (; resultTotal > k * X_binning_; ++k) {
+ sum += result;
+ X_[i] = sum;
+
+ for (; sum > k * X_binning_; ++k) {
d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
}
}
diff --git a/Environment/testShowerAxis.cc b/Environment/testShowerAxis.cc
index b541e63f0..8c2441b98 100644
--- a/Environment/testShowerAxis.cc
+++ b/Environment/testShowerAxis.cc
@@ -66,18 +66,18 @@ TEST_CASE("Homogeneous Density") {
Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
environment::ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos),
- *env, 10000};
+ *env, 20};
- CHECK(showerAxis.steplength() == 1_m);
+ CHECK(showerAxis.steplength() == 500_m);
CHECK(showerAxis.maximumX() / (10_km * density) == Approx(1).epsilon(1e-8));
CHECK(showerAxis.minimumX() == 0_g / square(1_cm));
- const Point p{cs, 10_km, 20_km, 9_km};
- CHECK(showerAxis.projectedX(p) / (1_km * density) == Approx(1).epsilon(1e-8));
+ const Point p{cs, 10_km, 20_km, 8.3_km};
+ CHECK(showerAxis.projectedX(p) / (1.7_km * density) == Approx(1).epsilon(1e-8));
- const units::si::LengthType d = 1_km;
+ const units::si::LengthType d = 6.789_km;
CHECK(showerAxis.X(d) / (d * density) == Approx(1).epsilon(1e-8));
const Vector<dimensionless_d> dir{cs, {0, 0, -1}};
--
GitLab