fixed interpolation

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_);

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);
         }
       }

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}};