implemented ShowerAxis

Environment/CMakeLists.txt

 set (
   ENVIRONMENT_SOURCES
   LayeredSphericalAtmosphereBuilder.cc
+  ShowerAxis.cc
 )
 
 set (
@@ -19,6 +20,7 @@ set (
   FlatExponential.h
   SlidingPlanarExponential.h
   LayeredSphericalAtmosphereBuilder.h
+  ShowerAxis.h
   )
 
 set (

Environment/ShowerAxis.cc

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/environment/ShowerAxis.h>
+
+using namespace corsika::environment;
+using namespace corsika::units::si;
+
+GrammageType ShowerAxis::X(LengthType l) const {
+  auto const fractionalBin = l / steplength_;
+  int const lower = fractionalBin; // indices of nearest X support points
+  auto const lambda = fractionalBin - lower;
+  int const upper = lower + 1;
+
+  if (upper >= steps) {
+    throw std::runtime_error("shower axis too short, cannot extrapolate");
+  } else if (lower < 0) {
+    throw std::runtime_error("cannot extrapolate to points behind point of injection");
+  }
+
+  assert(0 <= lambda && lambda <= 1.);
+  return (*X_)[lower] * lambda + (*X_)[upper] * (1 - lambda);
+}
+
+LengthType ShowerAxis::steplength() const { return steplength_; }
+
+GrammageType ShowerAxis::maximumX() const { return *(X_->rbegin()); }
+
+GrammageType ShowerAxis::minimumX() const { return *(X_->cbegin()); }
+
+GrammageType ShowerAxis::projectedX(geometry::Point const& p) const {
+  auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
+  return X(projectedLength);
+}

Environment/ShowerAxis.h

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/environment/Environment.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <cassert>
+#include <cstdlib>
+#include <fstream>
+#include <functional>
+#include <iterator>
+#include <memory>
+#include <stdexcept>
+#include <vector>
+
+#include <boost/math/quadrature/gauss_kronrod.hpp>
+
+namespace corsika::environment {
+  class ShowerAxis {
+  public:
+    template <typename TEnvModel>
+    ShowerAxis(geometry::Point const& pStart,
+               geometry::Vector<units::si::length_d> length,
+               environment::Environment<TEnvModel> const& env)
+        : pointStart_(pStart)
+        , length_(length)
+        , max_length_(length_.norm())
+        , steplength_(max_length_ / steps)
+        , axis_normalized_(length / max_length_) {
+      auto const* const universe = env.GetUniverse().get();
+
+      auto rho = [pStart, length, universe](double x) {
+        auto const p = pStart + length * x;
+        auto const* node = universe->GetContainingNode(p);
+        return node->GetModelProperties().GetMassDensity(p).magnitude();
+      };
+
+      double error;
+      int k = 0;
+      for (int i = 1; i < steps; ++i) {
+        auto const x_prev = (i - 1) / (steps - 1.);
+        auto const d_prev = max_length_ * x_prev;
+        auto const x = i / (steps - 1.);
+        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) {
+          d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
+        }
+      }
+
+      assert(std::is_sorted(X_->cbegin(), X_->cend()));
+      assert(std::is_sorted(d_.cbegin(), d_.cend()));
+    }
+
+    units::si::LengthType steplength() const;
+
+    units::si::GrammageType maximumX() const;
+
+    units::si::GrammageType minimumX() const;
+
+    units::si::GrammageType projectedX(geometry::Point const& p) const;
+
+    units::si::GrammageType X(units::si::LengthType) const;
+
+  private:
+    static constexpr int steps = 10'000;
+    units::si::GrammageType const X_binning_ = std::invoke([]() {
+      using namespace units::si;
+      return 1_g / 1_cm / 1_cm;
+    });
+    std::unique_ptr<std::array<units::si::GrammageType, steps>> const X_{
+        new std::array<units::si::GrammageType, steps>};
+    geometry::Point const pointStart_;
+    geometry::Vector<units::si::length_d> const length_;
+    units::si::LengthType const max_length_, steplength_;
+    geometry::Vector<units::si::dimensionless_d> const axis_normalized_;
+    std::vector<units::si::LengthType> d_;
+  };
+} // namespace corsika::environment