corsika/detail/framework/utility/LinearSolver.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+  inline std::vector<double> solve_linear_real(double a, double b) {
+
+    if (a == 0) {
+      return {}; // no (b!=0), or infinite number (b==0) of solutions....
+    }
+
+    return {-b / a};
+  }
+
+  inline std::vector<std::complex<double>> solve_linear(double a, double b) {
+
+    if (std::abs(a) == 0) {
+      return {}; // no (b!=0), or infinite number (b==0) of solutions....
+    }
+
+    return {std::complex<double>(-b / a, 0)};
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/QuadraticSolver.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+                                                           long double c,
+                                                           double const epsilon) {
+    if (std::abs(a) < epsilon) { return solve_linear(b, c); }
+
+    if (std::abs(c) < epsilon) {
+      std::vector<std::complex<double>> lin_result = solve_linear(a, b);
+      lin_result.push_back({0.});
+      return lin_result;
+    }
+
+    long double const radicant = static_pow<2>(b) - a * c * 4;
+
+    if (radicant < -epsilon) { // two complex solutions
+      double const rpart = -b / 2 * a;
+      double const ipart = std::sqrt(-radicant);
+      return {{rpart, ipart}, {rpart, -ipart}};
+    }
+
+    if (radicant < epsilon) { // just one real solution
+      return {{double(-b / 2 * a), 0}};
+    }
+
+    // two real solutions, use Citardauq formula and actively avoid cancellation in the
+    // denominator
+
+    const long double x1 =
+        2 * c / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+    return {{double(x1), 0}, {double(c / (a * x1)), 0}};
+  }
+
+  inline std::vector<double> solve_quadratic_real(long double a, long double b,
+                                                  long double c, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quadratic: a={} b={} c={}", a, b, c);
+
+    if (std::abs(a) < epsilon) { return solve_linear_real(b, c); }
+    if (std::abs(c) < epsilon) {
+      std::vector<double> lin_result = solve_linear_real(a, b);
+      lin_result.push_back(0.);
+      return lin_result;
+    }
+
+    // long double const radicant = std::pow(b, 2) - a * c * 4;
+    // Thanks!
+    // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+    long double w = a * 4 * c;
+    long double e = std::fma(a * 4, c, -w);
+    long double f = std::fma(b, b, -w);
+    long double radicant = f + e;
+
+    CORSIKA_LOG_TRACE("radicant={} {} ", radicant, b * b - a * c * 4);
+
+    if (std::abs(radicant) < epsilon) { // just one real solution
+      return {double(-b / (2 * a))};
+    }
+
+    if (radicant < 0) { return {}; }
+
+    // two real solutions, use Citardauq formula and actively avoid cancellation in the
+    // denominator
+
+    long double const x1 =
+        c * 2 / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+    CORSIKA_LOG_TRACE("quadratic x1={} x2={}", double(x1), double(c / (a * x1)));
+
+    return {double(x1), double(c / (a * x1))};
+  }
+} // namespace corsika

corsika/detail/framework/utility/QuarticSolver.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <cmath>
+
+namespace corsika {
+
+  namespace andre {
+
+    inline std::vector<double> solve_quartic_real(long double a, long double b,
+                                                  long double c, long double d,
+                                                  long double e, double const epsilon) {
+
+      if (std::abs(a) < epsilon) { return solve_cubic_real(b, c, d, e, epsilon); }
+
+      b /= a;
+      c /= a;
+      d /= a;
+      e /= a;
+
+      long double a3 = -c;
+      long double b3 = b * d - 4. * e;
+      long double c3 = -b * b * e - d * d + 4. * c * e;
+
+      // cubic resolvent
+      // y^3 − c*y^2 + (bd−4e)*y − b^2*e−d^2+4*c*e = 0
+
+      std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
+      if (!x3.size()) {
+        return {}; // no solution, numeric problem (LCOV_EXCL_LINE)
+      }
+      long double y = x3[0]; // there is always at least one solution
+      // The essence - choosing Y with maximal absolute value.
+      if (x3.size() == 3) {
+        if (std::abs(x3[1]) > std::abs(y)) y = x3[1];
+        if (std::abs(x3[2]) > std::abs(y)) y = x3[2];
+      }
+
+      long double q1, q2, p1, p2;
+      // h1+h2 = y && h1*h2 = e  <=>  h^2 -y*h + e = 0    (h === q)
+
+      long double Det = y * y - 4 * e;
+      CORSIKA_LOG_TRACE("Det={}", Det);
+      if (std::abs(Det) < epsilon) // in other words - D==0
+      {
+        q1 = q2 = y * 0.5;
+        // g1+g2 = b && g1+g2 = c-y   <=>   g^2 - b*g + c-y = 0    (p === g)
+        Det = b * b - 4 * (c - y);
+        CORSIKA_LOG_TRACE("Det={}", Det);
+        if (std::abs(Det) < epsilon) { // in other words - D==0
+          p1 = p2 = b * 0.5;
+        } else {
+          if (Det < 0) return {};
+          long double sqDet = sqrt(Det);
+          p1 = (b + sqDet) * 0.5;
+          p2 = (b - sqDet) * 0.5;
+        }
+      } else {
+        if (Det < 0) return {};
+        long double sqDet1 = sqrt(Det);
+        q1 = (y + sqDet1) * 0.5;
+        q2 = (y - sqDet1) * 0.5;
+        // g1+g2 = b && g1*h2 + g2*h1 = c       ( && g === p )  Krammer
+        p1 = (b * q1 - d) / (q1 - q2);
+        p2 = (d - b * q2) / (q1 - q2);
+      }
+
+      // solving quadratic eqs.  x^2 + p1*x + q1 = 0
+      //                         x^2 + p2*x + q2 = 0
+
+      std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
+      std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
+      if (quad2.size() > 0) {
+        for (auto const val : quad2) quad1.push_back(val);
+      }
+      return quad1;
+    }
+  } // namespace andre
+
+  inline std::vector<double> solve_quartic_depressed_real(long double p, long double q,
+                                                          long double r,
+                                                          double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f},  epsilon={}", p, q, r,
+                      epsilon);
+
+    long double const p2 = static_pow<2>(p);
+    long double const q2 = static_pow<2>(q);
+
+    std::vector<double> const resolve_cubic =
+        solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+    CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(),
+                      fmt::join(resolve_cubic, ", "));
+
+    if (!resolve_cubic.size()) return {};
+
+    long double m = 0;
+    for (auto const& v : resolve_cubic) {
+      CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+                                             v * (p2 / 4 - r) - q2 / 8));
+      if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
+    }
+    CORSIKA_LOG_TRACE("check m={}", m);
+    if (m == 0) { return {0}; }
+    if (m < 0) {
+      // this is a rare numerical instability
+      // first: try analytical solution, second: discard (curved->straight tracking)
+      std::vector<double> const resolve_cubic_analytic =
+          andre::solve_cubic_real_analytic(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+      CORSIKA_LOG_TRACE("andre::resolve_cubic_analytic: N={}, m=[{}]",
+                        resolve_cubic_analytic.size(),
+                        fmt::join(resolve_cubic_analytic, ", "));
+
+      if (!resolve_cubic_analytic.size()) return {};
+
+      for (auto const& v : resolve_cubic_analytic) {
+        CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+                                               v * (p2 / 4 - r) - q2 / 8));
+        if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
+      }
+      CORSIKA_LOG_TRACE("check m={}", m);
+      if (m == 0) { return {0}; }
+      if (m < 0) {
+        return {}; // now we are out of options, cannot solve: curved->straight tracking
+      }
+    }
+    long double const quad_term1 = p / 2 + m;
+    long double const quad_term2 = std::sqrt(2 * m);
+    long double const quad_term3 = q / (2 * quad_term2);
+
+    std::vector<double> z_quad1 =
+        solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, 1e-5);
+    std::vector<double> z_quad2 =
+        solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, 1e-5);
+    for (auto const& z : z_quad2) z_quad1.push_back(z);
+    return z_quad1;
+  }
+
+  inline std::vector<double> solve_quartic_real(long double a, long double b,
+                                                long double c, long double d,
+                                                long double e, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quartic: a={:f}, b={:f}, c={:f}, d={:f}, e={:f}, epsilon={}", a, b,
+                      c, d, e, epsilon);
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_cubic_real(b, c, d, e, epsilon);
+    }
+
+    if ((std::abs(a - 1) < epsilon) &&
+        (std::abs(b) < epsilon)) { // this is a depressed quartic
+      return solve_quartic_depressed_real(c, d, e, epsilon);
+    }
+
+    long double const b2 = static_pow<2>(b);
+    long double const b3 = static_pow<3>(b);
+    long double const b4 = static_pow<4>(b);
+    long double const a2 = static_pow<2>(a);
+    long double const a3 = static_pow<3>(a);
+    long double const a4 = static_pow<4>(a);
+
+    long double const p = (c * a * 8 - b2 * 3) / (a4 * 8);
+    long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8);
+    long double const r =
+        (-b4 * 3 + e * a3 * 256 - b * d * a2 * 64 + b2 * c * a * 16) / (a4 * 256);
+
+    std::vector<double> zs = solve_quartic_depressed_real(p, q, r, epsilon);
+    CORSIKA_LOG_TRACE("quartic: solve_depressed={}, b/4a={}", fmt::join(zs, ", "),
+                      b / (4 * a));
+    for (auto& z : zs) { z -= b / (4 * a); }
+    CORSIKA_LOG_TRACE("quartic: solve_quartic_real returns={}", fmt::join(zs, ", "));
+    return zs;
+  }
+} // namespace corsika

Framework/Utilities/SaveBoostHistogram.hpp → corsika/detail/framework/utility/SaveBoostHistogram.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * 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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
-#include <corsika/third_party/cnpy/cnpy.hpp>
+#pragma once
+
+#include <cnpy.hpp>
 
 #include <boost/histogram.hpp>
+#include <boost/filesystem.hpp> // can be changed to std::filesystem if compiler supports it
 
 #include <functional>
 #include <memory>
@@ -17,27 +19,22 @@
 #include <vector>
 #include <string>
 
-#pragma once
-
-namespace corsika::utl {
-
-  enum class SaveMode { overwrite, append };
+namespace corsika {
 
-  /**
-   * This functions saves a boost::histogram into a numpy file. Only rather basic axis
-   * types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
-   * counts (i.e. a double or int) are supported, nothing fancy like profiles.
-   *
-   * Note that this function makes a temporary, dense copy of the histogram, which could
-   * be an issue for huge sizes (e.g. for high dimensions)
-   */
   template <class Axes, class Storage>
   inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
-                        std::string const& filename, SaveMode mode = SaveMode::append) {
-    unsigned const rank = h.rank();
+                        std::string const& filename, bool overwrite) {
+    if (boost::filesystem::exists(filename)) {
+      if (overwrite) {
+        boost::filesystem::remove(filename);
+      } else {
+        using namespace std::literals;
+        throw std::runtime_error(
+            ("save_hist(): "s + filename + " already exists"s).c_str());
+      }
+    }
 
-    // append vs. overwrite
-    const std::string mode_str = (mode == SaveMode::append ? "a" : "w");
+    unsigned const rank = h.rank();
 
     std::vector<size_t> axes_dims;
     axes_dims.reserve(rank);
@@ -68,7 +65,7 @@ namespace corsika::utl {
         ax_edges.push_back(ax.bin(ax.size() - 1).upper());
 
         cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
-                       ax_edges.data(), {ax_edges.size()}, mode_str);
+                       ax_edges.data(), {ax_edges.size()}, "a");
       } else {
         axis_types.push_back('d');
         std::vector<int64_t> bins; // we assume that discrete axes have integer bins
@@ -77,14 +74,13 @@ namespace corsika::utl {
         for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
 
         cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
-                       {bins.size()}, mode_str);
+                       {bins.size()}, "a");
       }
     }
 
-    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank},
-                   mode_str);
-    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, mode_str);
-    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, mode_str);
+    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, "a");
 
     auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
                                                 unsigned{1}, std::multiplies<>());
@@ -108,8 +104,9 @@ namespace corsika::utl {
       temp[p] = *x;
     }
 
-    cnpy::npz_save(filename, "data", temp.get(), axes_dims, mode_str);
+    cnpy::npz_save(filename, "data", temp.get(), axes_dims, "a");
     // In Python this array can directly be assigned to a histogram view if that
     // histogram has its axes correspondingly: hist.view(flow=True)[:] = file['data']
-  }
-} // namespace corsika::utl
+  } // namespace corsika
+
+} // namespace corsika

corsika/detail/media/BaseExponential.inl

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
 
-#include <corsika/media/BaseExponential.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
 
 namespace corsika {
 
+  template <typename TDerived>
+  inline BaseExponential<TDerived>::BaseExponential(Point const& point,
+                                                    LengthType const referenceHeight,
+                                                    MassDensityType const rho0,
+                                                    LengthType const lambda)
+      : rho0_(rho0)
+      , lambda_(lambda)
+      , invLambda_(1 / lambda)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {}
 
-  template <class TDerived>
-    auto const& BaseExponential<TDerived>::GetImplementation() const 
-{ return *static_cast<TDerived const*>(this); }
+  template <typename TDerived>
+  inline auto const& BaseExponential<TDerived>::getImplementation() const {
+    return *static_cast<TDerived const*>(this);
+  }
 
-  template <class TDerived>
-    units::si::GrammageType BaseExponential<TDerived>::IntegratedGrammage(
-        Trajectory<Line> const& vLine,
-		units::si::LengthType vL,
-        Vector<units::si::dimensionless_d> const& vAxis) const {
-      if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
+  template <typename TDerived>
+  inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+      LengthType const height) const {
+    return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+  }
 
-      auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
-      auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+  template <typename TDerived>
+  inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj, DirectionVector const& axis) const {
+    LengthType const length = traj.getLength();
+    if (length == LengthType::zero()) { return GrammageType::zero(); }
 
-      if (uDotA == 0) {
-        return vL * rhoStart;
-      } else {
-        return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1);
-      }
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis);
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
+
+    if (uDotA == 0) {
+      return length * rhoStart;
+    } else {
+      return rhoStart * (lambda_ / uDotA) * expm1(uDotA * length * invLambda_);
     }
+  }
 
-  template <class TDerived>
-    units::si::LengthType BaseExponential<TDerived>::ArclengthFromGrammage(
-        Trajectory<Line> const& vLine,
-        units::si::GrammageType vGrammage,
-        Vector<units::si::dimensionless_d> const& vAxis) const {
-      auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
-      auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+  template <typename TDerived>
+  inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage,
+      DirectionVector const& axis) const {
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis);
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
 
-      if (uDotA == 0) {
-        return vGrammage / rhoStart;
+    if (uDotA == 0) {
+      return grammage / rhoStart;
+    } else {
+      auto const logArg = grammage * invLambda_ * uDotA / rhoStart;
+      if (logArg > -1) {
+        return lambda_ / uDotA * log1p(logArg);
       } else {
-        auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1;
-        if (logArg > 0) {
-          return fLambda / uDotA * log(logArg);
-        } else {
-          return std::numeric_limits<typename decltype(
-                     vGrammage)::value_type>::infinity() *
-                 units::si::meter;
-        }
+        return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
+               meter;
       }
     }
-
+  }
 
 } // namespace corsika
-

corsika/detail/media/BaseTabular.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+
+#include <exception>
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline BaseTabular<TDerived>::BaseTabular(
+      Point const& point, LengthType const referenceHeight,
+      std::function<MassDensityType(LengthType)> const& rho, unsigned int const nBins,
+      LengthType const deltaHeight)
+      : nBins_(nBins)
+      , deltaHeight_(deltaHeight)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {
+    density_.resize(nBins_);
+    for (unsigned int bin = 0; bin < nBins; ++bin) {
+      density_[bin] = rho(deltaHeight_ * bin);
+      CORSIKA_LOG_DEBUG("new tabulated atm bin={} h={} rho={}", bin, deltaHeight_ * bin,
+                        density_[bin]);
+    }
+  }
+
+  template <typename TDerived>
+  inline auto const& BaseTabular<TDerived>::getImplementation() const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline MassDensityType BaseTabular<TDerived>::getMassDensity(
+      LengthType const height) const {
+    double const fbin = (height - referenceHeight_) / deltaHeight_;
+    int const bin = int(fbin);
+    if (bin < 0) return MassDensityType::zero();
+    if (bin >= int(nBins_ - 1)) {
+      CORSIKA_LOG_ERROR(
+          "invalid height {} (corrected {}) in BaseTabular atmosphere. Min 0, max {}. If "
+          "max is too low: increase!",
+          height, height - referenceHeight_, nBins_ * deltaHeight_);
+      throw std::runtime_error("invalid height");
+    }
+    return density_[bin] + (fbin - bin) * (density_[bin + 1] - density_[bin]);
+  }
+
+  template <typename TDerived>
+  inline GrammageType BaseTabular<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+
+    Point pCurr = traj.getPosition(0);
+    DirectionVector dCurr = traj.getDirection(0);
+    LengthType height1 = (traj.getPosition(0) - point_).getNorm() - referenceHeight_;
+    LengthType height2 = (traj.getPosition(1) - point_).getNorm() - referenceHeight_;
+
+    LengthType const fullLength = traj.getLength(1);
+
+    int sign = +1; // normal direction
+    if (height1 > height2) {
+      std::swap(height1, height2);
+      pCurr = traj.getPosition(1);
+      dCurr = traj.getDirection(1);
+      sign = -1; // inverted direction
+    }
+
+    double const fbin1 = height1 / deltaHeight_;
+    unsigned int const bin1 = int(fbin1);
+
+    double const fbin2 = height2 / deltaHeight_;
+    unsigned int const bin2 = int(fbin2);
+
+    if (fbin1 == fbin2) { return GrammageType::zero(); }
+
+    if (bin1 >= nBins_ - 1 || bin2 >= nBins_ - 1) {
+      CORSIKA_LOG_ERROR("invalid height {} {} in BaseTabular atmosphere integration",
+                        height1, height2);
+      throw std::runtime_error("invalid height");
+    }
+
+    // interpolated start/end densities
+    MassDensityType const rho1 = getMassDensity(height1 + referenceHeight_);
+    MassDensityType const rho2 = getMassDensity(height2 + referenceHeight_);
+
+    // within first bin
+    if (bin1 == bin2) { return fullLength * (rho2 + rho1) / 2; }
+
+    // inclination of trajectory (local)
+    DirectionVector axis((pCurr - point_).normalized()); // to gravity center
+    double cosTheta = axis.dot(dCurr);
+
+    // distance to next height bin
+    unsigned int bin = bin1;
+    LengthType dD = (deltaHeight_ * (bin + 1) - height1) / cosTheta * sign;
+    LengthType distance = dD;
+
+    GrammageType X = dD * (rho1 + density_[bin + 1]) / 2;
+    double frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength);
+    pCurr = traj.getPosition(frac);
+    dCurr = traj.getDirection(frac);
+
+    for (++bin; bin < bin2; ++bin) {
+      // inclination of trajectory
+      axis = (pCurr - point_).normalized();
+      cosTheta = axis.dot(dCurr);
+      // distance to next height bin
+      dD = deltaHeight_ / cosTheta * sign;
+      distance += dD;
+      GrammageType const dX = dD * (density_[bin] + density_[bin + 1]) / 2;
+      X += dX;
+      frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength);
+      pCurr = traj.getPosition(frac);
+      dCurr = traj.getDirection(frac);
+    }
+
+    // inclination of trajectory
+    axis = ((pCurr - point_).normalized());
+    cosTheta = axis.dot(dCurr);
+    // distance to next height bin
+    dD = (height2 - deltaHeight_ * bin2) / cosTheta * sign;
+    X += dD * (rho2 + density_[bin2]) / 2;
+    distance += dD;
+    return X;
+  }
+
+  template <typename TDerived>
+  inline LengthType BaseTabular<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+
+    if (grammage < GrammageType::zero()) {
+      CORSIKA_LOG_ERROR("cannot integrate negative grammage");
+      throw std::runtime_error("negative grammage error");
+    }
+    LengthType const height = (traj.getPosition(0) - point_).getNorm() - referenceHeight_;
+
+    double const fbin = height / deltaHeight_;
+    int bin = int(fbin);
+
+    if (bin >= int(nBins_ - 1)) {
+      CORSIKA_LOG_ERROR("invalid height {} in BaseTabular atmosphere integration",
+                        height);
+      throw std::runtime_error("invalid height");
+    }
+
+    // interpolated start/end densities
+    MassDensityType const rho = getMassDensity(height + referenceHeight_);
+
+    // inclination of trajectory
+    Point pCurr = traj.getPosition(0);
+    DirectionVector dCurr = traj.getDirection(0);
+    DirectionVector axis((pCurr - point_).normalized());
+    double cosTheta = axis.dot(dCurr);
+    int sign = +1; // height increasing along traj
+    if (cosTheta < 0) {
+      cosTheta = -cosTheta; // absolute value only
+      sign = -1;            // height decreasing along traj
+    }
+
+    // height -> distance
+    LengthType const deltaDistance = deltaHeight_ / cosTheta;
+
+    // start with 0 g/cm2
+    GrammageType X(GrammageType::zero());
+    LengthType distance(LengthType::zero());
+
+    // within first bin
+    distance =
+        (sign > 0 ? deltaDistance * (bin + 1 - fbin) : deltaDistance * (fbin - bin));
+    GrammageType binGrammage = (sign > 0 ? distance * (rho + density_[bin + 1]) / 2
+                                         : distance * (rho + density_[bin]) / 2);
+    if (X + binGrammage > grammage) {
+      double const binFraction = (grammage - X) / binGrammage;
+      return distance * binFraction;
+    }
+    X += binGrammage;
+
+    // the following bins (along trajectory)
+    for (bin += sign; bin < int(nBins_) && bin >= 0; bin += sign) {
+
+      binGrammage = deltaDistance * (density_[bin] + density_[bin + 1]) / 2;
+      if (X + binGrammage > grammage) {
+        double const binFraction = (grammage - X) / binGrammage;
+        return distance + deltaDistance * binFraction;
+      }
+      X += binGrammage;
+      distance += deltaDistance;
+    }
+    return std::numeric_limits<double>::infinity() * meter;
+  }
+
+} // namespace corsika

corsika/detail/media/CORSIKA7Atmospheres.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  template <typename TEnvironmentInterface, template <typename> typename TExtraEnv,
+            typename TEnvironment, typename... TArgs>
+  void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
+                                Point const& center, TArgs... args) {
+
+    // construct the atmosphere builder
+    auto builder = make_layered_spherical_atmosphere_builder<
+        TEnvironmentInterface, TExtraEnv>::create(center, constants::EarthRadius::Mean,
+                                                  std::forward<TArgs>(args)...);
+
+    builder.setNuclearComposition(standardAirComposition);
+
+    // add the standard atmosphere layers
+    auto const params = atmosphereParameterList[static_cast<uint8_t>(atmId)];
+    for (int i = 0; i < 4; ++i) {
+      builder.addExponentialLayer(params[i].offset, params[i].scaleHeight,
+                                  params[i].altitude);
+    }
+    builder.addLinearLayer(params[4].offset, params[4].scaleHeight, params[4].altitude);
+
+    // and assemble the environment
+    builder.assemble(env);
+  }
+
+} // namespace corsika

corsika/detail/media/Environment.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -14,26 +11,61 @@
 
 namespace corsika {
 
-    template <typename IEnvironmentModel>
-    auto& Environment<IEnvironmentModel>::GetUniverse() {
- return fUniverse;     }
+  template <typename IEnvironmentModel>
+  inline Environment<IEnvironmentModel>::Environment()
+      : coordinateSystem_{get_root_CoordinateSystem()}
+      , universe_(std::make_unique<BaseNodeType>(
+            std::make_unique<Universe>(coordinateSystem_))) {}
+
+  template <typename IEnvironmentModel>
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr&
+  Environment<IEnvironmentModel>::getUniverse() {
+    return universe_;
+  }
+
+  template <typename IEnvironmentModel>
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr const&
+  Environment<IEnvironmentModel>::getUniverse() const {
+    return universe_;
+  }
+
+  template <typename IEnvironmentModel>
+  inline CoordinateSystemPtr const& Environment<IEnvironmentModel>::getCoordinateSystem()
+      const {
+    return coordinateSystem_;
+  }
 
-    template <typename IEnvironmentModel>
-    auto const& Environment<IEnvironmentModel>::GetUniverse() const {  return fUniverse; }
+  // factory method for creation of VolumeTreeNodes
+  template <typename IEnvironmentModel>
+  template <class TVolumeType, typename... TVolumeArgs>
+  std::unique_ptr<VolumeTreeNode<IEnvironmentModel> > inline Environment<
+      IEnvironmentModel>::createNode(TVolumeArgs&&... args) {
+    static_assert(std::is_base_of_v<IVolume, TVolumeType>,
+                  "unusable type provided, needs to be derived from "
+                  "\"Volume\"");
 
-    template <typename IEnvironmentModel>
-    auto const& Environment<IEnvironmentModel>::GetCoordinateSystem() const { return fCoordinateSystem; }
+    return std::make_unique<BaseNodeType>(
+        std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
+  }
 
-    // factory method for creation of VolumeTreeNodes
-    template <typename IEnvironmentModel>
-    template <class TVolumeType, typename... TVolumeArgs>
-    auto Environment<IEnvironmentModel>::CreateNode(TVolumeArgs&&... args) {
-      static_assert(std::is_base_of_v<corsika::Volume, TVolumeType>,
-                    "unusable type provided, needs to be derived from "
-                    "\"corsika::Volume\"");
+  template <typename IEnvironmentModel>
+  std::set<Code> const get_all_elements_in_universe(
+      Environment<IEnvironmentModel> const& env) {
 
-      return std::make_unique<BaseNodeType>(
-          std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
-    }
+    auto const& universe = *(env.getUniverse());
+    auto const allElementsInUniverse = std::invoke([&]() {
+      std::set<Code> allElementsInUniverse;
+      auto collectElements = [&](auto& vtn) {
+        if (vtn.hasModelProperties()) {
+          auto const& comp =
+              vtn.getModelProperties().getNuclearComposition().getComponents();
+          for (auto const c : comp) allElementsInUniverse.insert(c);
+        }
+      };
+      universe.walk(collectElements);
+      return allElementsInUniverse;
+    });
+    return allElementsInUniverse;
+  }
 
-}
\ No newline at end of file
+} // namespace corsika

corsika/detail/media/ExponentialRefractiveIndex.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
+      double const n0, InverseLengthType const lambda, Point const center,
+      LengthType const radius, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , n0_(n0)
+      , lambda_(lambda)
+      , center_(center)
+      , radius_(radius) {}
+
+  template <typename T>
+  inline double ExponentialRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
+  }
+
+} // namespace corsika

corsika/detail/media/FlatExponential.inl

 /*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
 
-#include <corsika/media/FlatExponential.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/BaseExponential.hpp>
+#include <corsika/media/NuclearComposition.hpp>
 
 namespace corsika {
 
-  template <class T>
-  units::si::MassDensityType FlatExponential<T>::GetMassDensity(Point const& vP) const {
-    return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
+  template <typename T>
+  inline FlatExponential<T>::FlatExponential(Point const& point,
+                                             DirectionVector const& axis,
+                                             MassDensityType const rho,
+                                             LengthType const lambda,
+                                             NuclearComposition const& nuclComp)
+      : BaseExponential<FlatExponential<T>>(point, 0_m, rho, lambda)
+      , axis_(axis)
+      , nuclComp_(nuclComp) {}
+
+  template <typename T>
+  inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
+    return BaseExponential<FlatExponential<T>>::getMassDensity(
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
   }
 
-  template <class T>
-  NuclearComposition const& FlatExponential<T>::GetNuclearComposition() const {
-    return fNuclComp;
+  template <typename T>
+  inline NuclearComposition const& FlatExponential<T>::getNuclearComposition() const {
+    return nuclComp_;
   }
 
-  template <class T>
-  units::si::GrammageType FlatExponential<T>::IntegratedGrammage(
-      Trajectory<Line> const& vLine, units::si::LengthType vTo) const {
-    return Base::IntegratedGrammage(vLine, vTo, fAxis);
+  template <typename T>
+  inline GrammageType FlatExponential<T>::getIntegratedGrammage(
+      BaseTrajectory const& line) const {
+    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
   }
 
-  template <class T>
-  units::si::LengthType FlatExponential<T>::ArclengthFromGrammage(
-      Trajectory<Line> const& vLine, units::si::GrammageType vGrammage) const {
-    return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
+  template <typename T>
+  inline LengthType FlatExponential<T>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType const grammage) const {
+    return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
+                                                                         axis_);
   }
 
 } // namespace corsika

corsika/detail/media/GeomagneticModel.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/framework/core/Logging.hpp>
+#include <boost/math/tr1.hpp>
+#include <boost/filesystem.hpp>
+#include <boost/filesystem/fstream.hpp>
+
+#include <stdexcept>
+#include <string>
+#include <cmath>
+
+namespace corsika {
+
+  inline GeomagneticModel::GeomagneticModel(Point const& center,
+                                            boost::filesystem::path const path)
+      : center_(center) {
+
+    // Read in coefficients
+    boost::filesystem::ifstream file(path, std::ios::in);
+
+    // Exit if file opening failed
+    if (!file.is_open()) {
+      CORSIKA_LOG_ERROR("Failed opening data file {}", path);
+      throw std::runtime_error("Cannot load GeomagneticModel data.");
+    }
+
+    // GeomagneticModel supports two types of input data: WMM.COF and IGRF.COF
+    // They have only slightly different format and content and can be easily
+    // differentiated here.
+
+    std::string line;
+    while (getline(file >> std::ws, line)) {
+
+      double epoch;
+      std::string model_name;
+      std::string release_date; // just for WMM
+      int nMax = 12; // the spherical max n (l) shell (for IGRF), for WMM this is 12
+                     // Note that n=l=0 is the monopole and is not included in the model.
+      int dummyInt;
+      double dummyDouble;
+
+      std::istringstream firstLine(line);
+
+      // check comments and ignore:
+      if (firstLine.peek() == '#' ||                     // normal comment
+          line.size() == 0 ||                            // empty line
+          line.find("9999999999999999999999999") == 0) { // crazy WMM comment
+        continue;
+      }
+
+      // check IGRF format:
+      if (firstLine >> model_name >> epoch >> nMax >> dummyInt >> dummyInt >>
+          dummyDouble >> dummyDouble >> dummyDouble >> dummyDouble >> model_name >>
+          dummyInt) {
+        static bool info = false;
+        if (!info) {
+          CORSIKA_LOG_INFO("Reading IGRF input data format.");
+          info = true;
+        }
+      } else {
+        // check WMM format:
+        firstLine.clear();
+        firstLine.seekg(0, std::ios::beg);
+        if (firstLine >> epoch >> model_name >> release_date) {
+          CORSIKA_LOG_INFO("Reading WMM input data format.");
+        } else {
+          CORSIKA_LOG_ERROR("line: {}", line);
+          throw std::runtime_error("Incompatible input data for GeomagneticModel");
+        }
+      }
+
+      int nPar = 0;
+      for (int i = 0; i < nMax; ++i) { nPar += i + 2; }
+      int iEpoch = int(epoch);
+
+      if (parameters_.count(iEpoch) != 0) {
+        throw std::runtime_error("GeomagneticModel input file has duplicate Epoch. Fix.");
+      }
+      parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
+
+      for (int i = 0; i < nPar; i++) {
+        file >> parameters_[iEpoch][i].n >> parameters_[iEpoch][i].m >>
+            parameters_[iEpoch][i].g >> parameters_[iEpoch][i].h >>
+            parameters_[iEpoch][i].dg >> parameters_[iEpoch][i].dh;
+        file.ignore(9999999, '\n');
+      }
+    }
+    file.close();
+
+    if (parameters_.size() == 0) {
+      CORSIKA_LOG_ERROR("No input data read!");
+      throw std::runtime_error("No input data read");
+    }
+  }
+
+  inline MagneticFieldVector GeomagneticModel::getField(double const year,
+                                                        LengthType const altitude,
+                                                        double const latitude,
+                                                        double const longitude) {
+
+    int iYear = int(year);
+    auto iEpoch = parameters_.rbegin();
+    for (; iEpoch != parameters_.rend(); ++iEpoch) {
+      if (iEpoch->first <= iYear) { break; }
+    }
+    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch->first, year);
+    if (iEpoch == parameters_.rend()) {
+      CORSIKA_LOG_WARN("Year {} is before first EPOCH. Results unclear.", year);
+      iEpoch--; // move one epoch back
+    }
+    if (altitude < -1_km || altitude > 600_km) {
+      CORSIKA_LOG_WARN("Altitude should be between -1_km and 600_km.");
+    }
+    if (latitude <= -90 || latitude >= 90) {
+      CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
+      throw std::runtime_error("Latitude has to be between -90 and 90 degree.");
+    } else if (latitude < -89.992 || latitude > 89.992) {
+      CORSIKA_LOG_WARN("Latitude is close to the poles.");
+    }
+    if (longitude < -180 || longitude > 180) {
+      CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
+    }
+
+    double const epoch = double(iEpoch->first);
+    auto iNextEpoch = iEpoch; // next epoch for interpolation
+    --iNextEpoch;
+    bool const lastEpoch = (iEpoch == parameters_.rbegin());
+    auto const delta_t = year - epoch;
+    CORSIKA_LOG_DEBUG(
+        "identified: t_epoch={}, delta_t={}, lastEpoch={} (false->interpolate)", epoch,
+        delta_t, lastEpoch);
+
+    double const lat_geo = latitude * constants::pi / 180;
+    double const lon = longitude * constants::pi / 180;
+
+    // Transform into spherical coordinates
+    double constexpr f = 1 / 298.257223563;
+    double constexpr e_squared = f * (2 - f);
+    LengthType R_c =
+        constants::EarthRadius::Equatorial / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
+    LengthType p = (R_c + altitude) * cos(lat_geo);
+    LengthType z = sin(lat_geo) * (altitude + R_c * (1 - e_squared));
+    LengthType r = sqrt(p * p + z * z);
+    double lat_sph = asin(z / r);
+
+    double legendre, next_legendre, derivate_legendre;
+    double magneticfield[3] = {0, 0, 0};
+
+    for (size_t j = 0; j < iEpoch->second.size(); j++) {
+
+      ParameterLine p = iEpoch->second[j];
+
+      // Time interpolation
+      if (iEpoch == parameters_.rbegin()) {
+        // this is the latest epoch in time, or time-dependence (dg/dh) was specified
+        // we use the extrapolation factors dg/dh:
+        p.g = p.g + delta_t * p.dg;
+        p.h = p.h + delta_t * p.dh;
+      } else {
+        // we linearly interpolate between two epochs
+        ParameterLine const next_p = iNextEpoch->second[j];
+        double const length = iNextEpoch->first - epoch;
+        double p_g = p.g + (next_p.g - p.g) * delta_t / length;
+        double p_h = p.h + (next_p.h - p.h) * delta_t / length;
+        CORSIKA_LOG_TRACE(
+            "interpolation: delta-g={}, delta-h={}, delta-t={}, length={} g1={} g2={} "
+            "g={} h={} ",
+            next_p.g - p.g, next_p.h - p.h, year - epoch, length, next_p.g, p.g, p_g,
+            p_h);
+        p.g = p_g;
+        p.h = p_h;
+      }
+
+      legendre = boost::math::tr1::assoc_legendre(p.n, p.m, sin(lat_sph));
+      next_legendre = boost::math::tr1::assoc_legendre(p.n + 1, p.m, sin(lat_sph));
+
+      // Schmidt semi-normalization
+      if (p.m > 0) {
+        // Note: n! = tgamma(n+1)
+        legendre *= sqrt(2 * std::tgamma(p.n - p.m + 1) / std::tgamma(p.n + p.m + 1));
+        next_legendre *=
+            sqrt(2 * std::tgamma(p.n + 1 - p.m + 1) / std::tgamma(p.n + 1 + p.m + 1));
+      }
+      derivate_legendre =
+          (p.n + 1) * tan(lat_sph) * legendre -
+          sqrt(pow(p.n + 1, 2) - pow(p.m, 2)) / cos(lat_sph) * next_legendre;
+
+      magneticfield[0] +=
+          pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+          (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * derivate_legendre;
+      magneticfield[1] +=
+          pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) * p.m *
+          (p.g * sin(p.m * lon) - p.h * cos(p.m * lon)) * legendre;
+      magneticfield[2] +=
+          (p.n + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+          (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * legendre;
+    }
+    magneticfield[0] *= -1;
+    magneticfield[1] /= cos(lat_sph);
+    magneticfield[2] *= -1;
+
+    // Transform back into geodetic coordinates
+    double magneticfield_geo[3];
+    magneticfield_geo[0] = magneticfield[0] * cos(lat_sph - lat_geo) -
+                           magneticfield[2] * sin(lat_sph - lat_geo);
+    magneticfield_geo[1] = magneticfield[1];
+    magneticfield_geo[2] = magneticfield[0] * sin(lat_sph - lat_geo) +
+                           magneticfield[2] * cos(lat_sph - lat_geo);
+
+    return MagneticFieldVector{center_.getCoordinateSystem(), magneticfield_geo[0] * 1_nT,
+                               magneticfield_geo[1] * -1_nT,
+                               magneticfield_geo[2] * -1_nT};
+  }
+
+} // namespace corsika

corsika/detail/media/GladstoneDaleRefractiveIndex.inl

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  inline GladstoneDaleRefractiveIndex<T>::GladstoneDaleRefractiveIndex(
+      double const referenceRefractiveIndex, Point const point, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , referenceRefractivity_(referenceRefractiveIndex - 1)
+      , referenceInvDensity_(1 / this->getMassDensity(point)) {}
+
+  template <typename T>
+  inline double GladstoneDaleRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return referenceRefractivity_ * (this->getMassDensity(point) * referenceInvDensity_) +
+           1.;
+  }
+
+} // namespace corsika

corsika/detail/media/HomogeneousMedium.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  inline HomogeneousMedium<T>::HomogeneousMedium(MassDensityType density,
+                                                 NuclearComposition const& nuclComp)
+      : density_(density)
+      , nuclComp_(nuclComp) {}
+
+  template <typename T>
+  inline MassDensityType HomogeneousMedium<T>::getMassDensity(Point const&) const {
+    return density_;
+  }
+
+  template <typename T>
+  inline NuclearComposition const& HomogeneousMedium<T>::getNuclearComposition() const {
+    return nuclComp_;
+  }
+
+  template <typename T>
+  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+      BaseTrajectory const& track) const {
+    return track.getLength() * density_;
+  }
+
+  template <typename T>
+  inline LengthType HomogeneousMedium<T>::getArclengthFromGrammage(
+      BaseTrajectory const&, GrammageType grammage) const {
+    return grammage / density_;
+  }
+} // namespace corsika

corsika/detail/media/InhomogeneousMedium.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika {
+
+  template <typename T, typename TDensityFunction>
+  template <typename... TArgs>
+  inline InhomogeneousMedium<T, TDensityFunction>::InhomogeneousMedium(
+      NuclearComposition const& nuclComp, TArgs&&... rhoTArgs)
+      : nuclComp_(nuclComp)
+      , densityFunction_(rhoTArgs...) {}
+
+  template <typename T, typename TDensityFunction>
+  inline MassDensityType InhomogeneousMedium<T, TDensityFunction>::getMassDensity(
+      Point const& point) const {
+    return densityFunction_.evaluateAt(point);
+  }
+
+  template <typename T, typename TDensityFunction>
+  inline NuclearComposition const&
+  InhomogeneousMedium<T, TDensityFunction>::getNuclearComposition() const {
+    return nuclComp_;
+  }
+
+  template <typename T, typename TDensityFunction>
+  inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
+      BaseTrajectory const& line) const {
+    return densityFunction_.getIntegrateGrammage(line);
+  }
+
+  template <typename T, typename TDensityFunction>
+  inline LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType grammage) const {
+    return densityFunction_.getArclengthFromGrammage(line, grammage);
+  }
+
+} // namespace corsika

corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  namespace detail {
+
+    struct NoExtraModelInner {};
+
+    template <typename M>
+    struct NoExtraModel {};
+
+    template <template <typename> typename M>
+    struct has_extra_models : std::true_type {};
+
+    template <>
+    struct has_extra_models<NoExtraModel> : std::false_type {};
+
+  } // namespace detail
+
+} // namespace corsika

corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
-#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
 #include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
 
 namespace corsika {
 
-void LayeredSphericalAtmosphereBuilder::checkRadius(units::si::LengthType r) const {
-  if (r <= previousRadius_) {
-    throw std::runtime_error("radius must be greater than previous");
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::checkRadius(LengthType const r)
+      const {
+    if (r <= previousRadius_) {
+      throw std::runtime_error("radius must be greater than previous");
+    }
   }
-}
-
-void LayeredSphericalAtmosphereBuilder::setNuclearComposition(
-    NuclearComposition composition) {
-  composition_ = std::make_unique<NuclearComposition>(composition);
-}
-
-void LayeredSphericalAtmosphereBuilder::addExponentialLayer(
-    units::si::GrammageType b, units::si::LengthType c,
-    units::si::LengthType upperBoundary) {
-  auto const radius = seaLevel_ + upperBoundary;
-  checkRadius(radius);
-  previousRadius_ = radius;
-
-  auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
-      std::make_unique<Sphere>(center_, radius));
-
-  auto const rho0 = b / c;
-  std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
-
-  node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(
-      center_, rho0, -c, *composition_, seaLevel_);
-
-  layers_.push(std::move(node));
-}
 
-void LayeredSphericalAtmosphereBuilder::addLinearLayer(
-    units::si::LengthType c, units::si::LengthType upperBoundary) {
-  using namespace units::si;
-
-  auto const radius = seaLevel_ + upperBoundary;
-  checkRadius(radius);
-  previousRadius_ = radius;
-
-  units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
-  auto const rho0 = b / c;
-
-  std::cout << "rho0 = " << rho0;
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra,
+      TModelArgs...>::setNuclearComposition(NuclearComposition const& composition) {
+    composition_ = std::make_unique<NuclearComposition>(composition);
+  }
 
-  auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
-      std::make_unique<Sphere>(center_, radius));
-  node->SetModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline typename LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
+                                                    TModelArgs...>::volume_tree_node*
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addExponentialLayer(GrammageType const b, LengthType const scaleHeight,
+                          LengthType const upperBoundary) {
+
+    // outer radius
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    auto const rho0 = b / scaleHeight;
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 5 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<
+            TMediumModelExtra<SlidingPlanarExponential<TMediumInterface>>>(
+            argPack..., center_, rho0, -scaleHeight, *composition_, planetRadius_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
+          center_, rho0, -scaleHeight, *composition_, planetRadius_);
+    }
+
+    layers_.push(std::move(node));
+    return layers_.top().get();
+  }
 
-  layers_.push(std::move(node));
-}
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra,
+      TModelArgs...>::addLinearLayer(GrammageType const b, LengthType const scaleHeight,
+                                     LengthType const upperBoundary) {
+    // outer radius
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    auto const rho0 = b / scaleHeight;
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 2 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<TMediumModelExtra<HomogeneousMedium<TMediumInterface>>>(
+            argPack..., rho0, *composition_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
+          rho0, *composition_);
+    }
+
+    layers_.push(std::move(node));
+  }
 
-Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
-  Environment<IMediumModel> env;
-  assemble(env);
-  return env;
-}
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+                      unsigned int const nBins, LengthType const deltaHeight,
+                      LengthType const upperBoundary) {
+
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 5 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<
+            TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
+            argPack..., center_, funcRho, nBins, deltaHeight, *composition_,
+            planetRadius_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
+          center_, funcRho, nBins, deltaHeight, *composition_, planetRadius_);
+    }
+    layers_.push(std::move(node));
+  }
 
-void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
-  auto& universe = env.GetUniverse();
-  auto* outmost = universe.get();
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::assemble() {
+    Environment<TMediumInterface> env;
+    assemble(env);
+    return env;
+  }
 
-  while (!layers_.empty()) {
-    auto l = std::move(layers_.top());
-    auto* tmp = l.get();
-    outmost->AddChild(std::move(l));
-    layers_.pop();
-    outmost = tmp;
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra,
+      TModelArgs...>::assemble(Environment<TMediumInterface>& env) {
+    auto& universe = env.getUniverse();
+    auto* outmost = universe.get();
+
+    while (!layers_.empty()) {
+      auto l = std::move(layers_.top());
+      auto* tmp = l.get();
+      outmost->addChild(std::move(l));
+      layers_.pop();
+      outmost = tmp;
+    }
   }
-}
+
+  template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
+  struct make_layered_spherical_atmosphere_builder {
+    template <typename... TArgs>
+    static auto create(Point const& center, LengthType const planetRadius,
+                       TArgs... args) {
+      return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
+                                               TArgs...>{std::forward<TArgs>(args)...,
+                                                         center, planetRadius};
+    }
+  };
 
 } // namespace corsika

corsika/detail/media/LinearApproximationIntegrator.inl

 /*
  * (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.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -14,44 +11,43 @@
 
 namespace corsika {
 
-
-   template <class TDerived>
-     auto const& LinearApproximationIntegrator<TDerived>::GetImplementation() const {
-		   return *static_cast<TDerived const*>(this);
-	   }
-
-
-   template <class TDerived>
-    auto LinearApproximationIntegrator<TDerived>::IntegrateGrammage(
-        corsika::Trajectory<corsika::Line> const& line,
-        corsika::units::si::LengthType length) const {
-      auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-      return (c0 + 0.5 * c1 * length) * length;
-    }
-
-   template <class TDerived>
-    auto LinearApproximationIntegrator<TDerived>::ArclengthFromGrammage(
-        corsika::Trajectory<corsika::Line> const& line,
-        corsika::units::si::GrammageType grammage) const {
-      auto const c0 = GetImplementation().fRho(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
-    }
-
-   template <class TDerived>
-    auto LinearApproximationIntegrator<TDerived>::MaximumLength(corsika::Trajectory<corsika::Line> const& line,
-                       [[maybe_unused]] double relError) const {
-      using namespace corsika::units::si;
-      [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      // todo: provide a real, working implementation
-      return 1_m * std::numeric_limits<double>::infinity();
-    }
-
-
-}
\ No newline at end of file
+  template <typename TDerived>
+  inline TDerived const& LinearApproximationIntegrator<TDerived>::getImplementation()
+      const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline GrammageType LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
+      BaseTrajectory const& line) const {
+    LengthType const length = line.getLength();
+    auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
+    auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
+                                                                line.getDirection(0));
+    CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1,
+                     line.getPosition(0), line.getDirection(0),
+                     (c0 + 0.5 * c1 * length) * length);
+    return (c0 + 0.5 * c1 * length) * length;
+  }
+
+  template <typename TDerived>
+  inline LengthType LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType grammage) const {
+    auto const c0 = getImplementation().rho_(line.getPosition(0));
+    auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
+                                                                line.getDirection(0));
+
+    return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
+  }
+
+  template <typename TDerived>
+  inline LengthType LinearApproximationIntegrator<TDerived>::getMaximumLength(
+      BaseTrajectory const& line, [[maybe_unused]] double relError) const {
+    [[maybe_unused]] auto const c1 = getImplementation().rho_.getSecondDerivative(
+        line.getPosition(0), line.getDirection(0));
+
+    // todo: provide a real, working implementation
+    return 1_m * std::numeric_limits<double>::infinity();
+  }
+
+} // namespace corsika

corsika/detail/media/MediumPropertyModel.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IMediumPropertyModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  inline MediumPropertyModel<T>::MediumPropertyModel(Medium const medium, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , medium_(medium) {}
+
+  template <typename T>
+  inline Medium MediumPropertyModel<T>::getMedium() const {
+    return medium_;
+  }
+
+  template <typename T>
+  inline void MediumPropertyModel<T>::setMedium(Medium const medium) {
+    medium_ = medium;
+  }
+
+} // namespace corsika

corsika/detail/media/NuclearComposition.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+
+#include <cassert>
+#include <functional>
+#include <numeric>
+#include <random>
+#include <stdexcept>
+#include <vector>
+
+namespace corsika {
+
+  inline NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
+                                                std::vector<double> const& pFractions)
+      : numberFractions_(pFractions)
+      , components_(pComponents)
+      , avgMassNumber_(getWeightedSum([](Code const compID) -> double {
+        if (is_nucleus(compID)) {
+          return get_nucleus_A(compID);
+        } else {
+          return get_mass(compID) / convert_SI_to_HEP(constants::u);
+        }
+      })) {
+    if (pComponents.size() != pFractions.size()) {
+      throw std::runtime_error(
+          "Cannot construct NuclearComposition from vectors of different sizes.");
+    }
+    auto const sumFractions = std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.);
+
+    if (!(0.999 < sumFractions && sumFractions < 1.001)) {
+      throw std::runtime_error("element fractions do not add up to 1");
+    }
+    this->updateHash();
+  }
+
+  template <typename TFunction>
+  inline auto NuclearComposition::getWeighted(TFunction func) const {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
+    auto const product = [&](auto const compID, auto const fraction) {
+      return func(compID) * fraction;
+    };
+
+    if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity::zero());
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    } else {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity(0));
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    }
+  } // namespace corsika
+
+  template <typename TFunction>
+  inline auto NuclearComposition::getWeightedSum(TFunction func) const
+      -> decltype(func(std::declval<Code>())) {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
+
+    auto const prod = [&](auto const compID, auto const fraction) {
+      return func(compID) * fraction;
+    };
+
+    if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+      return std::inner_product(
+          components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+          ResultQuantity::zero(), // .zero() is defined for quantity types only
+          std::plus<ResultQuantity>(), prod);
+    } else {
+      return std::inner_product(
+          components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+          ResultQuantity(0), // in other cases we have to use a bare 0
+          std::plus<ResultQuantity>(), prod);
+    }
+  }
+
+  inline size_t NuclearComposition::getSize() const { return numberFractions_.size(); }
+
+  inline std::vector<double> const& NuclearComposition::getFractions() const {
+    return numberFractions_;
+  }
+
+  inline std::vector<Code> const& NuclearComposition::getComponents() const {
+    return components_;
+  }
+
+  inline double const NuclearComposition::getAverageMassNumber() const {
+    return avgMassNumber_;
+  }
+
+  template <class TRNG>
+  inline Code NuclearComposition::sampleTarget(std::vector<CrossSectionType> const& sigma,
+                                               TRNG&& randomStream) const {
+    if (sigma.size() != numberFractions_.size()) {
+      throw std::runtime_error("incompatible vector sigma as input");
+    }
+
+    auto zip_beg = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cbegin(), sigma.cbegin()));
+    auto zip_end = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cend(), sigma.cend()));
+    using zip_iter_type = decltype(zip_beg);
+
+    auto const mult_func = [](zip_iter_type::value_type const& zipit) -> double {
+      return zipit.get<0>() * zipit.get<1>().magnitude();
+    };
+
+    using transform_iter_type =
+        boost::transform_iterator<decltype(mult_func), zip_iter_type, double, double>;
+
+    auto trans_beg = transform_iter_type{zip_beg, mult_func};
+    auto trans_end = transform_iter_type{zip_end, mult_func};
+
+    std::discrete_distribution channelDist{trans_beg, trans_end};
+
+    auto const iChannel = channelDist(randomStream);
+    return components_[iChannel];
+  }
+
+  // Note: when this class ever modifies its internal data, the hash
+  // must be updated, too!
+  // the hash value is important to find tables, etc.
+  inline size_t NuclearComposition::getHash() const { return hash_; }
+
+  inline bool NuclearComposition::operator==(NuclearComposition const& v) const {
+    return v.hash_ == hash_;
+  }
+
+  inline void NuclearComposition::updateHash() {
+    std::vector<std::size_t> hashes;
+    for (double ifrac : this->getFractions())
+      hashes.push_back(std::hash<double>{}(ifrac));
+    for (Code icode : this->getComponents())
+      hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
+    std::size_t h = std::hash<double>{}(this->getAverageMassNumber());
+    for (std::size_t ih : hashes) h = h ^ (ih << 1);
+    hash_ = h;
+  }
+
+} // namespace corsika

corsika/detail/media/ShowerAxis.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <string>
+
+namespace corsika {
+
+  template <typename TEnvModel>
+  inline ShowerAxis::ShowerAxis(Point const& pStart, Vector<length_d> const& length,
+                                Environment<TEnvModel> const& env, bool const doThrow,
+                                int const steps)
+      : pointStart_(pStart)
+      , length_(length)
+      , throw_(doThrow)
+      , max_length_(length_.getNorm())
+      , steplength_(max_length_ / steps)
+      , axis_normalized_(length / max_length_)
+      , X_(steps + 1) {
+    auto const* const universe = env.getUniverse().get();
+
+    auto rho = [pStart, length, universe, doThrow](double x) {
+      auto const p = pStart + length * x;
+      auto const* node = universe->getContainingNode(p);
+      if (!node->hasModelProperties()) {
+        CORSIKA_LOG_CRITICAL(
+            "Unable to construct ShowerAxis. ShowerAxis includes volume "
+            "with no model properties at point {}.",
+            p);
+        if (doThrow) throw std::runtime_error("Unable to construct ShowerAxis.");
+      }
+      return node->getModelProperties().getMassDensity(p).magnitude();
+    };
+
+    double error;
+    int k = 0;
+    X_[0] = GrammageType::zero();
+    auto sum = 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 = 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 =
+          MassDensityType(phys::units::detail::magnitude_tag, r) * max_length_;
+
+      sum += result;
+      X_[i] = sum;
+
+      for (; sum > 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()));
+  }
+
+  inline GrammageType ShowerAxis::getX(LengthType l) const {
+    double const fractionalBin = l / steplength_;
+    int const lower = fractionalBin; // indices of nearest X support points
+    double const fraction = fractionalBin - lower;
+    unsigned int const upper = lower + 1;
+
+    if (fractionalBin < 0) {
+      CORSIKA_LOG_TRACE("cannot extrapolate to points behind point of injection l={} m",
+                        l / 1_m);
+      if (throw_) {
+        throw std::runtime_error(
+            "cannot extrapolate to points behind point of injection");
+      }
+      return getMinimumX();
+    }
+
+    if (upper >= X_.size()) {
+      CORSIKA_LOG_TRACE(
+          "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+          l / max_length_);
+      if (throw_) {
+        const std::string err = fmt::format(
+            "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
+            l / max_length_);
+        throw std::runtime_error(err.c_str());
+      }
+      return getMaximumX();
+    }
+    CORSIKA_LOG_TRACE("showerAxis::X frac={}, fractionalBin={}, lower={}, upper={}",
+                      fraction, fractionalBin, lower, upper);
+
+    assert(0 <= fraction && fraction <= 1.);
+
+    CORSIKA_LOG_TRACE("ShowerAxis::getX l={} m, lower={}, fraction={}, upper={}", l / 1_m,
+                      lower, fraction, upper);
+
+    // linear interpolation between getX[lower] and X[upper]
+    return X_[upper] * fraction + X_[lower] * (1 - fraction);
+  }
+
+  inline LengthType ShowerAxis::getSteplength() const { return steplength_; }
+
+  inline GrammageType ShowerAxis::getMaximumX() const { return *X_.rbegin(); }
+
+  inline GrammageType ShowerAxis::getMinimumX() const { return GrammageType::zero(); }
+
+  inline GrammageType ShowerAxis::getProjectedX(Point const& p) const {
+    auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
+    return getX(projectedLength);
+  }
+
+  inline DirectionVector const& ShowerAxis::getDirection() const {
+    return axis_normalized_;
+  }
+
+  inline Point const& ShowerAxis::getStart() const { return pointStart_; }
+
+} // namespace corsika