corsika/detail/framework/random/random_iterator/detail/uint128.hpp

 /*
- * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ * (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.
  */
 
 /*
@@ -391,6 +390,7 @@ namespace random_iterator {
     void _init_dec(const char* s);
     void _init_oct(const char* s);
 
+#if defined(__powerpc64__) || defined(__x86_64__)
     static inline uint128_t mul128(uint128_t const x, uint128_t const y) {
       uint128_t z;
 #ifdef __powerpc64__
@@ -401,6 +401,7 @@ namespace random_iterator {
       z.UPPER += (x.UPPER * y.LOWER) + (x.LOWER * y.UPPER);
       return z;
     }
+#endif
 
 #ifdef __BIG_ENDIAN__
     uint64_t UPPER, LOWER;

corsika/detail/framework/random/random_iterator/detail/uint128.inl

 /*
- * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ * (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.
  */
 
 /*
@@ -211,11 +210,11 @@ namespace random_iterator {
   }
 
   inline bool uint128_t::operator>=(const uint128_t& rhs) const {
-    return ((*this > rhs) | (*this == rhs));
+    return ((*this > rhs) || (*this == rhs));
   }
 
   inline bool uint128_t::operator<=(const uint128_t& rhs) const {
-    return ((*this < rhs) | (*this == rhs));
+    return ((*this < rhs) || (*this == rhs));
   }
 
   inline uint128_t uint128_t::operator+(const uint128_t& rhs) const {

corsika/detail/framework/stack/CombinedStack.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

corsika/detail/framework/stack/SecondaryView.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

corsika/detail/framework/stack/Stack.inl

 /*
  * (c) Copyright 2018 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
@@ -348,7 +347,7 @@ namespace corsika {
   inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::swap(
       unsigned int const a, unsigned int const b) {
     data_.swap(a, b);
-    std::swap(deleted_[a], deleted_[b]);
+    std::vector<bool>::swap(deleted_[a], deleted_[b]);
   }
 
   template <typename StackData, template <typename> typename MParticleInterface,

corsika/detail/framework/utility/COMBoost.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,17 +18,15 @@
 
 namespace corsika {
 
-  inline COMBoost::COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
+  inline COMBoost::COMBoost(FourMomentum const& P4projectile,
                             HEPMassType const massTarget)
-      : originalCS_{Pprojectile.getSpaceLikeComponents().getCoordinateSystem()}
-      , rotatedCS_{
-            make_rotationToZ(Pprojectile.getSpaceLikeComponents().getCoordinateSystem(),
-                             Pprojectile.getSpaceLikeComponents())} {
-    auto const pProjectile = Pprojectile.getSpaceLikeComponents();
+      : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
+      , rotatedCS_{make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents())} {
+    auto const& pProjectile = P4projectile.getSpaceLikeComponents();
     auto const pProjNormSquared = pProjectile.getSquaredNorm();
     auto const pProjNorm = sqrt(pProjNormSquared);
 
-    auto const eProjectile = Pprojectile.getTimeLikeComponent();
+    auto const eProjectile = P4projectile.getTimeLikeComponent();
     auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
     auto const s =
         massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
@@ -39,12 +36,38 @@ namespace corsika {
     auto const coshEta = sqrt(1 + pProjNormSquared / s);
 
     setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
 
+  inline COMBoost::COMBoost(FourMomentum const& P4projectile,
+                            FourMomentum const& P4target)
+      : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()} {
+
+    // this is the center-of-momentum CM frame
+    auto const pCM =
+        P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents();
+    auto const pCM2 = pCM.getSquaredNorm();
+    auto const pCMnorm = sqrt(pCM2);
+    if (pCMnorm == 0_eV) {
+      // CM is at reset
+      rotatedCS_ = originalCS_;
+    } else {
+      rotatedCS_ = make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents() +
+                                                     P4target.getSpaceLikeComponents());
+    }
+
+    auto const s = (P4projectile + P4target).getNormSqr();
+    auto const sqrtS = sqrt(s);
+    auto const sinhEta = -pCMnorm / sqrtS;
+    auto const coshEta = sqrt(1 + pCM2 / s);
+
+    setBoost(coshEta, sinhEta);
     CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
                       coshEta, boost_.determinant() - 1);
   }
 
-  inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType mass)
+  inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType const mass)
       : originalCS_{momentum.getCoordinateSystem()}
       , rotatedCS_{make_rotationToZ(momentum.getCoordinateSystem(), momentum)} {
     auto const squaredNorm = momentum.getSquaredNorm();
@@ -52,15 +75,17 @@ namespace corsika {
     auto const sinhEta = -norm / mass;
     auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
     setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
   }
 
   template <typename FourVector>
-  inline FourVector COMBoost::toCoM(FourVector const& p) const {
-    auto pComponents = p.getSpaceLikeComponents().getComponents(rotatedCS_);
+  inline FourVector COMBoost::toCoM(FourVector const& p4) const {
+    auto const pComponents = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
     Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
     Eigen::Vector2d lab;
 
-    lab << (p.getTimeLikeComponent() * (1 / 1_GeV)),
+    lab << (p4.getTimeLikeComponent() * (1 / 1_GeV)),
         (eVecRotated(2) * (1 / 1_GeV).magnitude());
 
     auto const boostedZ = boost_ * lab;
@@ -68,21 +93,23 @@ namespace corsika {
 
     eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
 
+    CORSIKA_LOG_TRACE("E0={}, p={}, E0'={}, p'={}", p4.getTimeLikeComponent() / 1_GeV,
+                      eVecRotated(2) * (1 / 1_GeV).magnitude(), E_CoM / 1_GeV, boostedZ);
+
     return FourVector(E_CoM, MomentumVector(rotatedCS_, eVecRotated));
   }
 
   template <typename FourVector>
-  inline FourVector COMBoost::fromCoM(FourVector const& p) const {
-    auto pCM = p.getSpaceLikeComponents().getComponents(rotatedCS_);
-    auto const Ecm = p.getTimeLikeComponent();
+  inline FourVector COMBoost::fromCoM(FourVector const& p4) const {
+    auto pCM = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+    auto const Ecm = p4.getTimeLikeComponent();
 
     Eigen::Vector2d com;
     com << (Ecm * (1 / 1_GeV)), (pCM.getEigenVector()(2) * (1 / 1_GeV).magnitude());
 
-    CORSIKA_LOG_TRACE(
-        "COMBoost::fromCoM Ecm={} GeV"
-        " pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
-        Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV, p.getNorm() / 1_GeV);
+    CORSIKA_LOG_TRACE("Ecm={} GeV, pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
+                      Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV,
+                      p4.getNorm() / 1_GeV);
 
     auto const boostedZ = inverseBoost_ * com;
     auto const E_lab = boostedZ(0) * 1_GeV;
@@ -92,22 +119,24 @@ namespace corsika {
     Vector<typename decltype(pCM)::dimension_type> pLab{rotatedCS_, pCM};
     pLab.rebase(originalCS_);
 
-    FourVector f(E_lab, pLab);
-
     CORSIKA_LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
                       " plab={} GeV (norm={} GeV) "
                       " GeV), invariant mass = {}",
-                      E_lab / 1_GeV, f.getNorm() / 1_GeV, pLab.getComponents(),
-                      pLab.getNorm() / 1_GeV);
+                      E_lab / 1_GeV, FourVector{E_lab, pLab}.getNorm() / 1_GeV,
+                      pLab.getComponents(), pLab.getNorm() / 1_GeV);
 
-    return f;
+    return FourVector{E_lab, pLab};
   }
 
-  inline void COMBoost::setBoost(double coshEta, double sinhEta) {
+  inline void COMBoost::setBoost(double const coshEta, double const sinhEta) {
     boost_ << coshEta, sinhEta, sinhEta, coshEta;
     inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
   }
 
-  inline CoordinateSystemPtr COMBoost::getRotatedCS() const { return rotatedCS_; }
+  inline CoordinateSystemPtr const& COMBoost::getRotatedCS() const { return rotatedCS_; }
+
+  inline CoordinateSystemPtr const& COMBoost::getOriginalCS() const {
+    return originalCS_;
+  }
 
 } // namespace corsika

corsika/detail/framework/utility/CorsikaData.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/corsika.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
 #include <boost/filesystem/path.hpp>
 
 #include <cstdlib>
 #include <stdexcept>
 #include <string>
 
-inline boost::filesystem::path corsika::corsika_data(boost::filesystem::path const& key) {
-  if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
-    return boost::filesystem::path(p) / key;
-  } else { // LCOV_EXCL_START, this cannot be easily tested system-independently
-    throw std::runtime_error("CORSIKA_DATA not set");
-  } // LCOV_EXCL_STOP
-}
+namespace corsika {
+
+  inline boost::filesystem::path corsika_data(boost::filesystem::path const& key) {
+    boost::filesystem::path fname =
+        boost::filesystem::path(corsika::CORSIKA_DATA_DIR) / key;
+    // LCOV_EXCL_START, this cannot be easily tested system-independently
+    if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
+      fname = boost::filesystem::path(p) / key;
+    }
+    // LCOV_EXCL_STOP
+    CORSIKA_LOG_INFO("opening data file={}", fname);
+    return fname;
+  }
+} // namespace corsika

corsika/detail/framework/utility/CorsikaFenvFallback.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/utility/CorsikaFenv.hpp>

corsika/detail/framework/utility/CorsikaFenvOSX.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.
  */
 
 /**

corsika/detail/framework/utility/CubicSolver.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,6 +17,8 @@ namespace corsika {
 
     //---------------------------------------------------------------------------
     // solve cubic equation A x^3 + B*x^2 + C*x + D = 0
+    //                        x^3 + a*x^2 + b*x + c = 0
+    // mainly along WolframAlpha formulas
     inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
                                                          long double C, long double D,
                                                          double const epsilon) {
@@ -29,8 +30,8 @@ namespace corsika {
       long double c = D / A;
 
       long double a2 = a * a;
-      long double q = (a2 - 3 * b) / 9;
-      long double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+      long double q = (3 * b - a2) / 9;
+      long double r = (a * (9 * b - 2 * a2) - 27 * c) / 54;
       long double q3 = q * q * q;
 
       // disc = q**3 + r**2
@@ -38,8 +39,8 @@ namespace corsika {
       long double w = r * r;
       long double e = std::fma(r, r, -w);
       // s:t =  q*q exactly
-      long double s = -q * q;
-      long double t = std::fma(-q, q, -s);
+      long double s = q * q;
+      long double t = std::fma(q, q, -s);
       // s:t * q + w:e = s*q + w + t*q +e = s*q+w + u:v + e = f + u:v + e
       long double f = std::fma(s, q, w);
       long double u = t * q;
@@ -51,30 +52,34 @@ namespace corsika {
       au = f - au;
       ab = u - ab;
       // sum all terms into final result
-      long double const disc = -(((e + uf) + au) + ab) + v;
+      long double const disc = (((e + uf) + au) + ab) + v;
 
-      if (disc >= 0) {
-        long double t = r / std::sqrt(q3);
+      CORSIKA_LOG_TRACE("disc={} {}", disc, q3 + r * r);
+
+      if (std::abs(disc) < epsilon) {
+
+        a /= 3;
+        long double const cbrtR = std::cbrt(r);
+        return {double(2 * cbrtR - a), double(-cbrtR - a)}; // 2nd solution is doublet
+
+      } else if (disc > 0) {
+
+        long double const S = std::cbrt(r + std::sqrt(disc));
+        long double const T = std::cbrt(r - std::sqrt(disc));
+        a /= 3;
+        return {double((S + T) - a)}; // plus two imaginary solution
+
+      } else { // disc < 0
+
+        long double t = r / std::sqrt(-q3);
         if (t < -1) t = -1;
         if (t > 1) t = 1;
         t = std::acos(t);
         a /= 3;
-        q = -2 * std::sqrt(q);
+        q = 2 * std::sqrt(-q);
         return {double(q * std::cos(t / 3) - a),
                 double(q * std::cos((t + 2 * M_PI) / 3) - a),
-                double(q * std::cos((t - 2 * M_PI) / 3) - a)};
-      } else {
-        long double term1 = -cbrt(std::fabs(r) + std::sqrt(-disc));
-        if (r < 0) term1 = -term1;
-        long double term2 = (0 == term1 ? 0 : q / term1);
-
-        a /= 3;
-        long double test = 0.5 * std::sqrt(3.) * (term1 - term2);
-        if (std::fabs(test) < epsilon) {
-          return {double((term1 + term2) - 1), double(-0.5 * (term1 + term2) - a)};
-        }
-
-        return {double((term1 + term2) - a)};
+                double(q * std::cos((t + 4 * M_PI) / 3) - a)};
       }
     }
   } // namespace andre
@@ -228,7 +233,8 @@ namespace corsika {
   }
 
   /**
-   * Iterative approach.
+   * Iterative approach. https://en.wikipedia.org/wiki/Halley%27s_method
+   *  Halley's method
    */
 
   inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
@@ -238,72 +244,84 @@ namespace corsika {
                       a, b, c, d, epsilon, (std::abs(a - 1) < epsilon),
                       (std::abs(b) < epsilon));
 
-#ifdef DEBUG
-    {
-      auto test = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
-
-      for (long double test_v : test) {
-        CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
-                          cubic_function(test_v, a, b, c, d));
-      }
-    }
-#endif
-
     if (std::abs(a) < epsilon) { // this is just a quadratic
       return solve_quadratic_real(b, c, d, epsilon);
     }
 
-    long double const dist = std::fma(b / a, b / a, -3 * c / a);
-    long double const xinfl = -b / (a * 3);
+    auto pre_opt = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+    long double x1 = 0; // start value
 
-    long double x1 = xinfl;
-    long double f_x1 = cubic_function(xinfl, a, b, c, d);
-
-    if (std::abs(f_x1) > epsilon) {
-      if (std::abs(dist) < epsilon) {
-        x1 = xinfl - std::cbrt(f_x1);
-      } else if (dist > 0) {
-        if (f_x1 > 0)
-          x1 = xinfl - 2 / 3 * std::sqrt(dist);
-        else
-          x1 = xinfl + 2 / 3 * std::sqrt(dist);
+    if (pre_opt.size()) {
+      x1 = pre_opt[0]; //*std::max_element(pre_opt.begin(), pre_opt.end());
+#ifdef _C8_DEBUG_
+      for (long double test_v : pre_opt) {
+        CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
+                          cubic_function(test_v, a, b, c, d));
       }
-
-      int niter = 0;
-      const int maxiter = 100;
-      do {
-        long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
-        long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
-        // if (potential) saddle point... avoid
-        if (std::abs(f_prime_x1) < epsilon) {
-          x1 -= std::cbrt(f_x1);
-        } else {
-          x1 -=
-              f_x1 * f_prime_x1 / (static_pow<2>(f_prime_x1) - 0.5 * f_x1 * f_prime2_x1);
+#endif
+    } else {
+      // this is only if the former solve_cubic_real_analytic would not result
+      // in any solution. We have no test case for this. This is excluded from tests:
+      // LCOV_EXCL_START
+      long double const dist = std::fma(b / a, b / a, -3 * c / a);
+      long double const xinfl = -b / (a * 3);
+
+      x1 = xinfl;
+      long double f_test = cubic_function(xinfl, a, b, c, d);
+
+      if (std::abs(f_test) > epsilon) {
+        if (std::abs(dist) < epsilon) {
+          x1 = xinfl - std::cbrt(f_test);
+        } else if (dist > 0) {
+          if (f_test > 0)
+            x1 = xinfl - 2 / 3 * std::sqrt(dist);
+          else
+            x1 = xinfl + 2 / 3 * std::sqrt(dist);
         }
-        f_x1 = cubic_function(x1, a, b, c, d);
-        CORSIKA_LOG_TRACE("niter={} x1={} f_x1={} f_prime={} f_prime2={} eps={}", niter,
-                          x1, f_x1, f_prime_x1, f_prime2_x1, epsilon);
-      } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon));
-
-      CORSIKA_LOG_TRACE("niter={}", niter);
-      if (niter >= maxiter) {
-        CORSIKA_LOG_TRACE("failure, no solution");
-        //        return std::vector<double>{};
       }
+      // LCOV_EXCL_STOP
+    }
+
+    long double f_x1 = cubic_function(x1, a, b, c, d);
+    long double dx1 = 0;
+
+    int niter = 0;
+    const int maxiter = 100;
+    do {
+      long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
+      long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
+      // if (potential) saddle point... avoid
+      if (std::abs(f_prime_x1) < epsilon) {
+        dx1 = std::cbrt(f_x1);
+      } else {
+        dx1 = f_x1 * f_prime_x1 * 2 / (f_prime_x1 * f_prime_x1 * 2 - f_x1 * f_prime2_x1);
+      }
+      x1 -= dx1;
+      f_x1 = cubic_function(x1, a, b, c, d);
+      CORSIKA_LOG_TRACE(
+          "niter={} x1={:.20f} f_x1={:.20f} f_prime={:.20f} f_prime2={:.20f} dx1={}",
+          niter, x1, f_x1, f_prime_x1, f_prime2_x1,
+          f_x1 * f_prime_x1 / (f_prime_x1 * f_prime_x1 - f_x1 * f_prime2_x1 * 0.5));
+    } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon * 1000) &&
+             (std::abs(dx1) > epsilon));
+
+    CORSIKA_LOG_TRACE("niter={}", niter);
+    if (niter >= maxiter) {
+      CORSIKA_LOG_DEBUG("niter reached max iterations {}", niter);
+      return andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
     }
 
     CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1);
 
     double const b1 = x1 + b / a;
     double const b0 = b1 * x1 + c / a;
-    std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+    std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, 1e-3);
     CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
                       cubic_function(x1, a, b, c, d));
 
     quad_check.push_back(x1);
     CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(quad_check, ", "));
     return quad_check;
-  }
+  } // namespace corsika
 
 } // namespace corsika

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 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/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
@@ -36,13 +35,13 @@ namespace corsika {
 
       std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
       if (!x3.size()) {
-        return {}; // no solution, numeric problem
+        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;
@@ -50,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 {};
@@ -77,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;
     }
@@ -112,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

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