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 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/framework/core/PhysicalUnits.hpp>
@@ -62,7 +61,7 @@ namespace corsika {
     long double f = std::fma(b, b, -w);
     long double radicant = f + e;
 
-    CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4);
+    CORSIKA_LOG_TRACE("radicant={} {} ", radicant, b * b - a * c * 4);
 
     if (std::abs(radicant) < epsilon) { // just one real solution
       return {double(-b / (2 * a))};

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 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
@@ -35,11 +34,14 @@ namespace corsika {
       // 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 (fabs(x3[1]) > fabs(y)) y = x3[1];
-        if (fabs(x3[2]) > fabs(y)) y = x3[2];
+        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;
@@ -47,13 +49,13 @@ namespace corsika {
 
       long double Det = y * y - 4 * e;
       CORSIKA_LOG_TRACE("Det={}", Det);
-      if (fabs(Det) < epsilon) // in other words - D==0
+      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 (fabs(Det) < epsilon) { // in other words - D==0
+        if (std::abs(Det) < epsilon) { // in other words - D==0
           p1 = p2 = b * 0.5;
         } else {
           if (Det < 0) return {};
@@ -74,10 +76,10 @@ namespace corsika {
       // 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);
-      std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
+      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 val : quad2) quad1.push_back(val);
+        for (auto const val : quad2) quad1.push_back(val);
       }
       return quad1;
     }
@@ -109,17 +111,37 @@ namespace corsika {
     }
     CORSIKA_LOG_TRACE("check m={}", m);
     if (m == 0) { return {0}; }
-
-    CORSIKA_LOG_TRACE("check m={}", m);
-
+    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, epsilon);
+        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, epsilon);
+        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;
   }

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.
  */
 
 #pragma once
@@ -25,8 +24,7 @@ namespace corsika {
   template <class Axes, class Storage>
   inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
                         std::string const& filename, bool overwrite) {
-    if (boost::filesystem::status(filename).type() !=
-        boost::filesystem::file_type::file_not_found) {
+    if (boost::filesystem::exists(filename)) {
       if (overwrite) {
         boost::filesystem::remove(filename);
       } else {

corsika/detail/media/BaseExponential.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.
  */
 
 #pragma once
@@ -14,39 +13,61 @@
 
 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 <typename TDerived>
   inline auto const& BaseExponential<TDerived>::getImplementation() const {
     return *static_cast<TDerived const*>(this);
   }
 
+  template <typename TDerived>
+  inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+      LengthType const height) const {
+    return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+  }
+
   template <typename TDerived>
   inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
-      BaseTrajectory const& traj, LengthType vL, DirectionVector const& axis) const {
-    if (vL == LengthType::zero()) { return GrammageType::zero(); }
+      BaseTrajectory const& traj, DirectionVector const& axis) const {
+    LengthType const length = traj.getLength();
+    if (length == LengthType::zero()) { return GrammageType::zero(); }
 
-    auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
-    auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis);
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
 
     if (uDotA == 0) {
-      return vL * rhoStart;
+      return length * rhoStart;
     } else {
-      return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1);
+      return rhoStart * (lambda_ / uDotA) * expm1(uDotA * length * invLambda_);
     }
   }
 
   template <typename TDerived>
   inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
-      BaseTrajectory const& traj, GrammageType grammage,
+      BaseTrajectory const& traj, GrammageType const grammage,
       DirectionVector const& axis) const {
-    auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
-    auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis);
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
 
     if (uDotA == 0) {
       return grammage / rhoStart;
     } else {
-      auto const logArg = grammage * invLambda_ * uDotA / rhoStart + 1;
-      if (logArg > 0) {
-        return lambda_ / uDotA * log(logArg);
+      auto const logArg = grammage * invLambda_ * uDotA / rhoStart;
+      if (logArg > -1) {
+        return lambda_ / uDotA * log1p(logArg);
       } else {
         return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
                meter;
@@ -54,13 +75,4 @@ namespace corsika {
     }
   }
 
-  template <typename TDerived>
-  inline BaseExponential<TDerived>::BaseExponential(Point const& point,
-                                                    MassDensityType rho0,
-                                                    LengthType lambda)
-      : rho0_(rho0)
-      , lambda_(lambda)
-      , invLambda_(1 / lambda)
-      , point_(point) {}
-
 } // 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
  *
- * 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
@@ -49,4 +48,24 @@ namespace corsika {
         std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
   }
 
+  template <typename IEnvironmentModel>
+  std::set<Code> const get_all_elements_in_universe(
+      Environment<IEnvironmentModel> const& env) {
+
+    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;
+  }
+
 } // 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 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
@@ -15,14 +14,18 @@ namespace corsika {
   template <typename T>
   template <typename... Args>
   ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
-      double const n0, InverseLengthType const lambda, Args&&... args)
+      double const n0, InverseLengthType const lambda, Point const center,
+      LengthType const radius, Args&&... args)
       : T(std::forward<Args>(args)...)
-      , n_0(n0)
-      , lambda_(lambda) {}
+      , n0_(n0)
+      , lambda_(lambda)
+      , center_(center)
+      , radius_(radius) {}
 
   template <typename T>
-  double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
-    return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+  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 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.
  */
 
 #pragma once
@@ -18,19 +17,18 @@ namespace corsika {
 
   template <typename T>
   inline FlatExponential<T>::FlatExponential(Point const& point,
-                                             Vector<dimensionless_d> const& axis,
-                                             MassDensityType rho, LengthType lambda,
+                                             DirectionVector const& axis,
+                                             MassDensityType const rho,
+                                             LengthType const lambda,
                                              NuclearComposition const& nuclComp)
-      : BaseExponential<FlatExponential<T>>(point, rho, lambda)
+      : 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>>::getRho0() *
-           exp(BaseExponential<FlatExponential<T>>::getInvLambda() *
-               (point - BaseExponential<FlatExponential<T>>::getAnchorPoint())
-                   .dot(axis_));
+    return BaseExponential<FlatExponential<T>>::getMassDensity(
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
   }
 
   template <typename T>
@@ -40,13 +38,13 @@ namespace corsika {
 
   template <typename T>
   inline GrammageType FlatExponential<T>::getIntegratedGrammage(
-      BaseTrajectory const& line, LengthType to) const {
-    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
+      BaseTrajectory const& line) const {
+    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
   }
 
   template <typename T>
   inline LengthType FlatExponential<T>::getArclengthFromGrammage(
-      BaseTrajectory const& line, GrammageType grammage) const {
+      BaseTrajectory const& line, GrammageType const grammage) const {
     return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
                                                                          axis_);
   }

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 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
@@ -32,9 +31,9 @@ namespace corsika {
   }
 
   template <typename T>
-  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(BaseTrajectory const&,
-                                                                  LengthType to) const {
-    return to * density_;
+  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+      BaseTrajectory const& track) const {
+    return track.getLength() * density_;
   }
 
   template <typename T>

corsika/detail/media/InhomogeneousMedium.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.
  */
 
 #pragma once
@@ -36,8 +35,8 @@ namespace corsika {
 
   template <typename T, typename TDensityFunction>
   inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
-      BaseTrajectory const& line, LengthType to) const {
-    return densityFunction_.getIntegrateGrammage(line, to);
+      BaseTrajectory const& line) const {
+    return densityFunction_.getIntegrateGrammage(line);
   }
 
   template <typename T, typename TDensityFunction>

corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp

 /*
  * (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.
  */
 
 #pragma once

corsika/detail/media/LayeredSphericalAtmosphereBuilder.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.
  */
 
 #pragma once
@@ -13,13 +12,14 @@
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
 #include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
 
 namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  inline void LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
-                                                TModelArgs...>::checkRadius(LengthType r)
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::checkRadius(LengthType const r)
       const {
     if (r <= previousRadius_) {
       throw std::runtime_error("radius must be greater than previous");
@@ -36,26 +36,28 @@ namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  inline void LayeredSphericalAtmosphereBuilder<
-      TMediumInterface, TMediumModelExtra,
-      TModelArgs...>::addExponentialLayer(GrammageType b, LengthType c,
-                                          LengthType upperBoundary) {
-
-    auto const radius = earthRadius_ + upperBoundary;
+  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 / c;
+    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, -c, *composition_, earthRadius_);
+            argPack..., center_, rho0, -scaleHeight, *composition_, planetRadius_);
       };
 
       // now unpack the additional arguments
@@ -63,26 +65,28 @@ namespace corsika {
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
-          center_, rho0, -c, *composition_, earthRadius_);
+          center_, rho0, -scaleHeight, *composition_, planetRadius_);
     }
 
     layers_.push(std::move(node));
+    return layers_.top().get();
   }
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
   inline void LayeredSphericalAtmosphereBuilder<
       TMediumInterface, TMediumModelExtra,
-      TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
-    auto const radius = earthRadius_ + upperBoundary;
+      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));
 
-    units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
-    auto const rho0 = b / c;
+    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
@@ -93,7 +97,6 @@ namespace corsika {
 
       // now unpack the additional arguments
       auto model = std::apply(lastBound, additionalModelArgs_);
-
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
@@ -103,6 +106,40 @@ namespace corsika {
     layers_.push(std::move(node));
   }
 
+  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));
+  }
+
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
   inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
@@ -132,10 +169,11 @@ namespace corsika {
   template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
   struct make_layered_spherical_atmosphere_builder {
     template <typename... TArgs>
-    static auto create(Point const& center, LengthType earthRadius, TArgs... args) {
+    static auto create(Point const& center, LengthType const planetRadius,
+                       TArgs... args) {
       return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
                                                TArgs...>{std::forward<TArgs>(args)...,
-                                                         center, earthRadius};
+                                                         center, planetRadius};
     }
   };
 

corsika/detail/media/LinearApproximationIntegrator.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.
  */
 
 #pragma once
@@ -13,21 +12,26 @@
 namespace corsika {
 
   template <typename TDerived>
-  inline auto const& LinearApproximationIntegrator<TDerived>::getImplementation() const {
+  inline TDerived const& LinearApproximationIntegrator<TDerived>::getImplementation()
+      const {
     return *static_cast<TDerived const*>(this);
   }
 
   template <typename TDerived>
-  inline auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
-      BaseTrajectory const& line, LengthType length) const {
+  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 auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
+  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),
@@ -37,7 +41,7 @@ namespace corsika {
   }
 
   template <typename TDerived>
-  inline auto LinearApproximationIntegrator<TDerived>::getMaximumLength(
+  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));

corsika/detail/media/MediumPropertyModel.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.
  */
 
 #pragma once
@@ -19,7 +18,7 @@ namespace corsika {
       , medium_(medium) {}
 
   template <typename T>
-  inline Medium MediumPropertyModel<T>::getMedium(Point const&) const {
+  inline Medium MediumPropertyModel<T>::getMedium() const {
     return medium_;
   }
 

corsika/detail/media/NuclearComposition.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.
  */
 
 #pragma once
@@ -11,7 +10,8 @@
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
-#include <corsika/media/WeightProvider.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
 
 #include <cassert>
 #include <functional>
@@ -23,31 +23,52 @@
 namespace corsika {
 
   inline NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
-                                                std::vector<float> const& pFractions)
+                                                std::vector<double> const& pFractions)
       : numberFractions_(pFractions)
       , components_(pComponents)
-      , avgMassNumber_(std::inner_product(
-            pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
-            std::plus<double>(), [](auto const compID, auto const fraction) -> double {
-              if (is_nucleus(compID)) {
-                return get_nucleus_A(compID) * fraction;
-              } else {
-                return get_mass(compID) / convert_SI_to_HEP(constants::u) * fraction;
-              }
-            })) {
-    assert(pComponents.size() == pFractions.size());
-    auto const sumFractions =
-        std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-    if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
+      , 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::getWeightedSum(TFunction const& func) const {
-    using ResultQuantity = decltype(func(*components_.cbegin()));
+  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;
@@ -68,7 +89,7 @@ namespace corsika {
 
   inline size_t NuclearComposition::getSize() const { return numberFractions_.size(); }
 
-  inline std::vector<float> const& NuclearComposition::getFractions() const {
+  inline std::vector<double> const& NuclearComposition::getFractions() const {
     return numberFractions_;
   }
 
@@ -82,16 +103,28 @@ namespace corsika {
 
   template <class TRNG>
   inline Code NuclearComposition::sampleTarget(std::vector<CrossSectionType> const& sigma,
-                                               TRNG& randomStream) const {
+                                               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>;
 
-    assert(sigma.size() == numberFractions_.size());
+    auto trans_beg = transform_iter_type{zip_beg, mult_func};
+    auto trans_end = transform_iter_type{zip_end, mult_func};
 
-    std::discrete_distribution channelDist(
-        WeightProviderIterator<decltype(numberFractions_.begin()),
-                               decltype(sigma.begin())>(numberFractions_.begin(),
-                                                        sigma.begin()),
-        WeightProviderIterator<decltype(numberFractions_.begin()), decltype(sigma.end())>(
-            numberFractions_.end(), sigma.end()));
+    std::discrete_distribution channelDist{trans_beg, trans_end};
 
     auto const iChannel = channelDist(randomStream);
     return components_[iChannel];
@@ -99,6 +132,7 @@ namespace corsika {
 
   // 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 {
@@ -107,7 +141,8 @@ namespace corsika {
 
   inline void NuclearComposition::updateHash() {
     std::vector<std::size_t> hashes;
-    for (float ifrac : this->getFractions()) hashes.push_back(std::hash<float>{}(ifrac));
+    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());

corsika/detail/media/ShowerAxis.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/framework/core/PhysicalUnits.hpp>
@@ -26,9 +25,16 @@ namespace corsika {
       , X_(steps + 1) {
     auto const* const universe = env.getUniverse().get();
 
-    auto rho = [pStart, length, universe](double x) {
+    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();
     };
 
@@ -65,7 +71,7 @@ namespace corsika {
     unsigned int const upper = lower + 1;
 
     if (fractionalBin < 0) {
-      CORSIKA_LOG_ERROR("cannot extrapolate to points behind point of injection l={} m",
+      CORSIKA_LOG_TRACE("cannot extrapolate to points behind point of injection l={} m",
                         l / 1_m);
       if (throw_) {
         throw std::runtime_error(
@@ -75,7 +81,7 @@ namespace corsika {
     }
 
     if (upper >= X_.size()) {
-      CORSIKA_LOG_ERROR(
+      CORSIKA_LOG_TRACE(
           "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
           l / max_length_);
       if (throw_) {

corsika/detail/media/SlidingPlanarExponential.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.
  */
 
 #pragma once
 
-#include <corsika/media/SlidingPlanarExponential.hpp>
-
 namespace corsika {
 
-  template <typename T>
-  inline SlidingPlanarExponential<T>::SlidingPlanarExponential(
-      Point const& p0, MassDensityType rho0, LengthType lambda,
-      NuclearComposition const& nuclComp, LengthType referenceHeight)
-      : BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
-      , nuclComp_(nuclComp)
-      , referenceHeight_(referenceHeight) {}
+  template <typename TDerived>
+  inline SlidingPlanarExponential<TDerived>::SlidingPlanarExponential(
+      Point const& p0, MassDensityType const rho0, LengthType const lambda,
+      NuclearComposition const& nuclComp, LengthType const referenceHeight)
+      : BaseExponential<SlidingPlanarExponential<TDerived>>(p0, referenceHeight, rho0,
+                                                            lambda)
+      , nuclComp_(nuclComp) {}
 
-  template <typename T>
-  inline MassDensityType SlidingPlanarExponential<T>::getMassDensity(
+  template <typename TDerived>
+  inline MassDensityType SlidingPlanarExponential<TDerived>::getMassDensity(
       Point const& point) const {
-    auto const height =
-        (point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-            .getNorm() -
-        referenceHeight_;
-    return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
-           exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * height);
+    auto const heightFull =
+        (point - BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .getNorm();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getMassDensity(
+        heightFull);
   }
 
-  template <typename T>
-  inline NuclearComposition const& SlidingPlanarExponential<T>::getNuclearComposition()
-      const {
+  template <typename TDerived>
+  inline NuclearComposition const&
+  SlidingPlanarExponential<TDerived>::getNuclearComposition() const {
     return nuclComp_;
   }
 
-  template <typename T>
-  inline GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
-      BaseTrajectory const& traj, LengthType l) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, l,
-                                                                               axis);
+  template <typename TDerived>
+  inline GrammageType SlidingPlanarExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getIntegratedGrammage(
+        traj, axis);
   }
 
-  template <typename T>
-  inline LengthType SlidingPlanarExponential<T>::getArclengthFromGrammage(
+  template <typename TDerived>
+  inline LengthType SlidingPlanarExponential<TDerived>::getArclengthFromGrammage(
       BaseTrajectory const& traj, GrammageType const grammage) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getArclengthFromGrammage(
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getArclengthFromGrammage(
         traj, grammage, axis);
   }