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.
  */
 
 /**
@@ -28,7 +27,7 @@
 
 extern "C" {
 
-int feenableexcept(int excepts) throw() {
+inline int feenableexcept(int excepts) noexcept {
   static fenv_t fenv;
   int new_excepts = excepts & FE_ALL_EXCEPT;
   // previous masks
@@ -44,7 +43,7 @@ int feenableexcept(int excepts) throw() {
   return fesetenv(&fenv) ? -1 : old_excepts;
 }
 
-int fedisableexcept(int excepts) throw() {
+inline int fedisableexcept(int excepts) noexcept {
   static fenv_t fenv;
   int new_excepts = excepts & FE_ALL_EXCEPT;
   // all previous masks

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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <cmath>
+
+namespace corsika {
+
+  namespace andre {
+
+    //---------------------------------------------------------------------------
+    // 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) {
+
+      if (std::abs(A) < epsilon) { return solve_quadratic_real(B, C, epsilon); }
+
+      long double a = B / A;
+      long double b = C / A;
+      long double c = D / A;
+
+      long double a2 = a * a;
+      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
+      // w:e = r*r exactly
+      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);
+      // 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;
+      long double v = std::fma(t, q, -u);
+      // error-free sum f+u. See Knuth, TAOCP vol. 2
+      long double uf = u + f;
+      long double au = uf - u;
+      long double ab = uf - au;
+      au = f - au;
+      ab = u - ab;
+      // sum all terms into final result
+      long double const disc = (((e + uf) + au) + ab) + v;
+
+      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);
+        return {double(q * std::cos(t / 3) - a),
+                double(q * std::cos((t + 2 * M_PI) / 3) - a),
+                double(q * std::cos((t + 4 * M_PI) / 3) - a)};
+      }
+    }
+  } // namespace andre
+
+  inline std::vector<double> solve_cubic_depressed_disciminant_real(
+      long double p, long double q, long double const disc, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("p={}, q={}, disc={}", p, q, disc);
+
+    if (std::abs(disc) < epsilon) { // disc==0   multiple roots !
+      if (std::abs(p) < epsilon) {  // tripple root
+        return {0};
+      }
+      // double root, single root
+      CORSIKA_LOG_TRACE("cubic double root");
+      return {double(-3 * q / (2 * p)), double(3 * q / p)};
+    }
+
+    if (std::abs(p) < epsilon) { // p==0  --> x^3 + q = 0
+      return {double(std::cbrt(-q))};
+    }
+
+    if (disc > 0) { // three real roots
+      CORSIKA_LOG_TRACE("cubic three roots");
+      long double const p_third = p / 3;
+      long double const sqrt_minus_p_third = std::sqrt(-p_third);
+
+      long double const arg = std::acos(q / (2 * p_third) / sqrt_minus_p_third) / 3;
+
+      long double constexpr pi = M_PI;
+      return {double(2 * sqrt_minus_p_third * std::cos(arg)),
+              double(2 * sqrt_minus_p_third * std::cos(arg - 2 * pi / 3)),
+              double(2 * sqrt_minus_p_third * std::cos(arg - 4 * pi / 3))};
+    }
+
+    // thus disc < 0  ->  one real root
+    if (p < 0) {
+      CORSIKA_LOG_TRACE("cubic cosh");
+      long double const abs_q = std::abs(q);
+      long double const p_third = p / 3;
+      long double const sqrt_minus_p_third = std::sqrt(-p_third);
+      CORSIKA_LOG_TRACE("sqrt_minus_p_third={}, arcosh({})={}", sqrt_minus_p_third,
+                        -abs_q / (2 * p_third) / sqrt_minus_p_third,
+                        std::acosh(-abs_q / (2 * p_third) / sqrt_minus_p_third));
+      CORSIKA_LOG_TRACE(
+          "complex: {}",
+          -2 * abs_q / q * sqrt_minus_p_third *
+              std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3));
+      double const z =
+          -2 * abs_q / q * sqrt_minus_p_third *
+          std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3);
+      return {z};
+    } else { // p > 0
+      CORSIKA_LOG_TRACE("cubic sinh");
+      long double const p_third = p / 3;
+      long double const sqrt_p_third = std::sqrt(p_third);
+      return {double(-2 * sqrt_p_third *
+                     std::sinh(std::asinh(q / (2 * p_third * sqrt_p_third)) / 3))};
+    }
+  }
+
+  inline std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+                                                        double const epsilon) {
+
+    // thanks!
+    // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+    // long double const disc = -(4 * p * p * p + 27 * q * q);
+    // disc = (p/3)**3 + (q/2)**2
+    long double p_third = p / 3;
+    long double q_half = q / 2;
+    // w:e = (q/2)*(q/2) exactly
+    long double w = q_half * q_half;
+    long double e = std::fma(q_half, q_half, -w);
+    // s:t =  (p/3)*(p/3) exactly
+    long double s = p_third * p_third;
+    long double t = std::fma(p_third, p_third, -s);
+    // s:t * p + w:e = s*p + w + t*p +e = s*p+w + u:v + e = f + u:v + e
+    long double f = std::fma(s, p_third, w);
+    long double u = t * p_third;
+    long double v = std::fma(t, p_third, -u);
+    // error-free sum f+u. See Knuth, TAOCP vol. 2
+    long double a = u + f;
+    long double b = a - u;
+    long double c = a - b;
+    b = f - b;
+    c = u - c;
+    // sum all terms into final result
+    long double const disc = -(((e + a) + b) + c) + v;
+    return solve_cubic_depressed_disciminant_real(p, q, disc, epsilon);
+  }
+
+  /**
+   * Analytical approach. Not very stable in all conditions.
+   */
+  inline std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+                                                       long double c, long double d,
+                                                       double const epsilon) {
+
+    CORSIKA_LOG_TRACE("cubic: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b, c,
+                      d, epsilon, (std::abs(a - 1) < epsilon), (std::abs(b) < epsilon));
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_quadratic_real(b, c, d, epsilon);
+    }
+
+    if ((std::abs(a - 1) < epsilon) &&
+        (std::abs(b) < epsilon)) { // this is a depressed cubic
+      return solve_cubic_depressed_real(c, d, epsilon);
+    }
+
+    // p = (3ac - b^2) / (3a^2) = 3 * ( c/(3*a) - b^2/(9*a^2) )
+    long double b_over_a = b / a;
+    long double const p_third = std::fma(-b_over_a, b_over_a / 9, c / (a * 3));
+
+    // p = std::fma(a * 3, c, -b2) / (3 * a2);
+    // q = (2*b^3 - 9*abc + 27*a^2*d) / (27a^3) = 2 * ( b^3/(27*a^3) - bc/(6*a^2) +
+    // d/(2*a) )
+    long double const q_half_term1 = std::fma(b_over_a, b_over_a / 27, -c / (a * 6));
+    long double const q_half = std::fma(b_over_a, q_half_term1, d / (a * 2));
+
+    std::vector<double> zs = solve_cubic_depressed_real(3 * p_third, 2 * q_half, epsilon);
+    CORSIKA_LOG_TRACE("cubic: solve_depressed={}, b/3a={}", fmt::join(zs, ", "),
+                      b / (a * 3));
+    for (auto& z : zs) {
+      z -= b / (a * 3);
+      double const b1 = z + b / a;
+      double const b0 = b1 * z + c / a;
+      std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+      CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+                        static_pow<3>(z) * a + static_pow<2>(z) * b + z * c + d);
+    }
+    CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(zs, ", "));
+    return zs;
+  }
+
+  template <typename T> // T must be floating point type
+  inline T cubic_function(T x, T a, T b, T c, T d) {
+    T x2 = x * x;
+    return x2 * x * a + x2 * b + x * c + d;
+  }
+
+  template <typename T> // T must be floating point type
+  inline T cubic_function_dfdx(T x, T a, T b, T c) {
+    T x2 = x * x;
+    return x2 * a * 3 + x * b * 2 + c;
+  }
+
+  template <typename T> // T must be floating point type
+  inline T cubic_function_d2fd2x(T x, T a, T b) {
+    return x * a * 6 + b * 2;
+  }
+
+  /**
+   * 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,
+                                              long double d, double const epsilon) {
+
+    CORSIKA_LOG_TRACE("cubic_iterative: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}",
+                      a, b, c, d, epsilon, (std::abs(a - 1) < epsilon),
+                      (std::abs(b) < epsilon));
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_quadratic_real(b, c, d, epsilon);
+    }
+
+    auto pre_opt = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+    long double x1 = 0; // start value
+
+    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));
+      }
+#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);
+        }
+      }
+      // 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, 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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+  inline std::vector<double> solve_linear_real(double a, double b) {
+
+    if (a == 0) {
+      return {}; // no (b!=0), or infinite number (b==0) of solutions....
+    }
+
+    return {-b / a};
+  }
+
+  inline std::vector<std::complex<double>> solve_linear(double a, double b) {
+
+    if (std::abs(a) == 0) {
+      return {}; // no (b!=0), or infinite number (b==0) of solutions....
+    }
+
+    return {std::complex<double>(-b / a, 0)};
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/QuadraticSolver.inl

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

corsika/detail/framework/utility/QuarticSolver.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
 #include <cmath>
 
-namespace corsika::quartic_solver {
-
-  //---------------------------------------------------------------------------
-  // solve cubic equation x^3 + a*x^2 + b*x + c = 0
-  // x - array of size 3
-  // In case 3 real roots: => x[0], x[1], x[2], return 3
-  //         2 real roots: x[0], x[1],          return 2
-  //         1 real root : x[0], x[1] ± i*x[2], return 1
-  inline unsigned int solveP3(double* x, double a, double b, double c) {
-    double a2 = a * a;
-    double q = (a2 - 3 * b) / 9;
-    double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
-    double r2 = r * r;
-    double q3 = q * q * q;
-    double A, B;
-    if (r2 < q3) {
-      double t = r / sqrt(q3);
-      if (t < -1) t = -1;
-      if (t > 1) t = 1;
-      t = acos(t);
-      a /= 3;
-      q = -2 * sqrt(q);
-      x[0] = q * cos(t / 3) - a;
-      x[1] = q * cos((t + M_2PI) / 3) - a;
-      x[2] = q * cos((t - M_2PI) / 3) - a;
-      return 3;
-    } else {
-      A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
-      if (r < 0) A = -A;
-      B = (0 == A ? 0 : q / A);
-
-      a /= 3;
-      x[0] = (A + B) - a;
-      x[1] = -0.5 * (A + B) - a;
-      x[2] = 0.5 * sqrt(3.) * (A - B);
-      if (fabs(x[2]) < eps) {
-        x[2] = x[1];
-        return 2;
+namespace corsika {
+
+  namespace andre {
+
+    inline std::vector<double> solve_quartic_real(long double a, long double b,
+                                                  long double c, long double d,
+                                                  long double e, double const epsilon) {
+
+      if (std::abs(a) < epsilon) { return solve_cubic_real(b, c, d, e, epsilon); }
+
+      b /= a;
+      c /= a;
+      d /= a;
+      e /= a;
+
+      long double a3 = -c;
+      long double b3 = b * d - 4. * e;
+      long double c3 = -b * b * e - d * d + 4. * c * e;
+
+      // cubic resolvent
+      // y^3 − c*y^2 + (bd−4e)*y − b^2*e−d^2+4*c*e = 0
+
+      std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
+      if (!x3.size()) {
+        return {}; // no solution, numeric problem (LCOV_EXCL_LINE)
+      }
+      long double y = x3[0]; // there is always at least one solution
+      // The essence - choosing Y with maximal absolute value.
+      if (x3.size() == 3) {
+        if (std::abs(x3[1]) > std::abs(y)) y = x3[1];
+        if (std::abs(x3[2]) > std::abs(y)) y = x3[2];
       }
 
-      return 1;
+      long double q1, q2, p1, p2;
+      // h1+h2 = y && h1*h2 = e  <=>  h^2 -y*h + e = 0    (h === q)
+
+      long double Det = y * y - 4 * e;
+      CORSIKA_LOG_TRACE("Det={}", Det);
+      if (std::abs(Det) < epsilon) // in other words - D==0
+      {
+        q1 = q2 = y * 0.5;
+        // g1+g2 = b && g1+g2 = c-y   <=>   g^2 - b*g + c-y = 0    (p === g)
+        Det = b * b - 4 * (c - y);
+        CORSIKA_LOG_TRACE("Det={}", Det);
+        if (std::abs(Det) < epsilon) { // in other words - D==0
+          p1 = p2 = b * 0.5;
+        } else {
+          if (Det < 0) return {};
+          long double sqDet = sqrt(Det);
+          p1 = (b + sqDet) * 0.5;
+          p2 = (b - sqDet) * 0.5;
+        }
+      } else {
+        if (Det < 0) return {};
+        long double sqDet1 = sqrt(Det);
+        q1 = (y + sqDet1) * 0.5;
+        q2 = (y - sqDet1) * 0.5;
+        // g1+g2 = b && g1*h2 + g2*h1 = c       ( && g === p )  Krammer
+        p1 = (b * q1 - d) / (q1 - q2);
+        p2 = (d - b * q2) / (q1 - q2);
+      }
+
+      // solving quadratic eqs.  x^2 + p1*x + q1 = 0
+      //                         x^2 + p2*x + q2 = 0
+
+      std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
+      std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
+      if (quad2.size() > 0) {
+        for (auto const val : quad2) quad1.push_back(val);
+      }
+      return quad1;
     }
-  }
+  } // namespace andre
 
-  //---------------------------------------------------------------------------
-  // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
-  // Attention - this function returns dynamically allocated array. It has to be released
-  // afterwards.
-  inline DComplex* solve_quartic(double a, double b, double c, double d) {
-    double a3 = -b;
-    double b3 = a * c - 4. * d;
-    double c3 = -a * a * d - c * c + 4. * b * d;
+  inline std::vector<double> solve_quartic_depressed_real(long double p, long double q,
+                                                          long double r,
+                                                          double const epsilon) {
 
-    // cubic resolvent
-    // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+    CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f},  epsilon={}", p, q, r,
+                      epsilon);
 
-    double x3[3];
-    unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+    long double const p2 = static_pow<2>(p);
+    long double const q2 = static_pow<2>(q);
 
-    double q1, q2, p1, p2, D, sqD, y;
+    std::vector<double> const resolve_cubic =
+        solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
 
-    y = x3[0];
-    // The essence - choosing Y with maximal absolute value.
-    if (iZeroes != 1) {
-      if (fabs(x3[1]) > fabs(y)) y = x3[1];
-      if (fabs(x3[2]) > fabs(y)) y = x3[2];
-    }
+    CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(),
+                      fmt::join(resolve_cubic, ", "));
 
-    // h1+h2 = y && h1*h2 = d  <=>  h^2 -y*h + d = 0    (h === q)
-
-    D = y * y - 4 * d;
-    if (fabs(D) < eps) // in other words - D==0
-    {
-      q1 = q2 = y * 0.5;
-      // g1+g2 = a && g1+g2 = b-y   <=>   g^2 - a*g + b-y = 0    (p === g)
-      D = a * a - 4 * (b - y);
-      if (fabs(D) < eps) // in other words - D==0
-        p1 = p2 = a * 0.5;
-
-      else {
-        sqD = sqrt(D);
-        p1 = (a + sqD) * 0.5;
-        p2 = (a - sqD) * 0.5;
+    if (!resolve_cubic.size()) return {};
+
+    long double m = 0;
+    for (auto const& v : resolve_cubic) {
+      CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+                                             v * (p2 / 4 - r) - q2 / 8));
+      if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
+    }
+    CORSIKA_LOG_TRACE("check m={}", m);
+    if (m == 0) { return {0}; }
+    if (m < 0) {
+      // this is a rare numerical instability
+      // first: try analytical solution, second: discard (curved->straight tracking)
+      std::vector<double> const resolve_cubic_analytic =
+          andre::solve_cubic_real_analytic(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+      CORSIKA_LOG_TRACE("andre::resolve_cubic_analytic: N={}, m=[{}]",
+                        resolve_cubic_analytic.size(),
+                        fmt::join(resolve_cubic_analytic, ", "));
+
+      if (!resolve_cubic_analytic.size()) return {};
+
+      for (auto const& v : resolve_cubic_analytic) {
+        CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
+                                               v * (p2 / 4 - r) - q2 / 8));
+        if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
+      }
+      CORSIKA_LOG_TRACE("check m={}", m);
+      if (m == 0) { return {0}; }
+      if (m < 0) {
+        return {}; // now we are out of options, cannot solve: curved->straight tracking
       }
-    } else {
-      sqD = sqrt(D);
-      q1 = (y + sqD) * 0.5;
-      q2 = (y - sqD) * 0.5;
-      // g1+g2 = a && g1*h2 + g2*h1 = c       ( && g === p )  Krammer
-      p1 = (a * q1 - c) / (q1 - q2);
-      p2 = (c - a * q2) / (q1 - q2);
     }
+    long double const quad_term1 = p / 2 + m;
+    long double const quad_term2 = std::sqrt(2 * m);
+    long double const quad_term3 = q / (2 * quad_term2);
+
+    std::vector<double> z_quad1 =
+        solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, 1e-5);
+    std::vector<double> z_quad2 =
+        solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, 1e-5);
+    for (auto const& z : z_quad2) z_quad1.push_back(z);
+    return z_quad1;
+  }
+
+  inline std::vector<double> solve_quartic_real(long double a, long double b,
+                                                long double c, long double d,
+                                                long double e, double const epsilon) {
 
-    DComplex* retval = new DComplex[4];
-
-    // solving quadratic eq. - x^2 + p1*x + q1 = 0
-    D = p1 * p1 - 4 * q1;
-    if (D < 0.0) {
-      retval[0].real(-p1 * 0.5);
-      retval[0].imag(sqrt(-D) * 0.5);
-      retval[1] = std::conj(retval[0]);
-    } else {
-      sqD = sqrt(D);
-      retval[0].real((-p1 + sqD) * 0.5);
-      retval[1].real((-p1 - sqD) * 0.5);
+    CORSIKA_LOG_TRACE("quartic: a={:f}, b={:f}, c={:f}, d={:f}, e={:f}, epsilon={}", a, b,
+                      c, d, e, epsilon);
+
+    if (std::abs(a) < epsilon) { // this is just a quadratic
+      return solve_cubic_real(b, c, d, e, epsilon);
     }
 
-    // solving quadratic eq. - x^2 + p2*x + q2 = 0
-    D = p2 * p2 - 4 * q2;
-    if (D < 0.0) {
-      retval[2].real(-p2 * 0.5);
-      retval[2].imag(sqrt(-D) * 0.5);
-      retval[3] = std::conj(retval[2]);
-    } else {
-      sqD = sqrt(D);
-      retval[2].real((-p2 + sqD) * 0.5);
-      retval[3].real((-p2 - sqD) * 0.5);
+    if ((std::abs(a - 1) < epsilon) &&
+        (std::abs(b) < epsilon)) { // this is a depressed quartic
+      return solve_quartic_depressed_real(c, d, e, epsilon);
     }
 
-    return retval;
+    long double const b2 = static_pow<2>(b);
+    long double const b3 = static_pow<3>(b);
+    long double const b4 = static_pow<4>(b);
+    long double const a2 = static_pow<2>(a);
+    long double const a3 = static_pow<3>(a);
+    long double const a4 = static_pow<4>(a);
+
+    long double const p = (c * a * 8 - b2 * 3) / (a4 * 8);
+    long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8);
+    long double const r =
+        (-b4 * 3 + e * a3 * 256 - b * d * a2 * 64 + b2 * c * a * 16) / (a4 * 256);
+
+    std::vector<double> zs = solve_quartic_depressed_real(p, q, r, epsilon);
+    CORSIKA_LOG_TRACE("quartic: solve_depressed={}, b/4a={}", fmt::join(zs, ", "),
+                      b / (4 * a));
+    for (auto& z : zs) { z -= b / (4 * a); }
+    CORSIKA_LOG_TRACE("quartic: solve_quartic_real returns={}", fmt::join(zs, ", "));
+    return zs;
   }
-
-} // namespace corsika::quartic_solver
\ No newline at end of file
+} // namespace corsika

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
@@ -11,6 +10,7 @@
 #include <cnpy.hpp>
 
 #include <boost/histogram.hpp>
+#include <boost/filesystem.hpp> // can be changed to std::filesystem if compiler supports it
 
 #include <functional>
 #include <memory>
@@ -23,11 +23,18 @@ namespace corsika {
 
   template <class Axes, class Storage>
   inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
-                        std::string const& filename, SaveMode mode) {
-    unsigned const rank = h.rank();
+                        std::string const& filename, bool overwrite) {
+    if (boost::filesystem::exists(filename)) {
+      if (overwrite) {
+        boost::filesystem::remove(filename);
+      } else {
+        using namespace std::literals;
+        throw std::runtime_error(
+            ("save_hist(): "s + filename + " already exists"s).c_str());
+      }
+    }
 
-    // append vs. overwrite
-    const std::string mode_str = (mode == SaveMode::append ? "a" : "w");
+    unsigned const rank = h.rank();
 
     std::vector<size_t> axes_dims;
     axes_dims.reserve(rank);
@@ -58,7 +65,7 @@ namespace corsika {
         ax_edges.push_back(ax.bin(ax.size() - 1).upper());
 
         cnpy::npz_save(filename, std::string{"binedges_"} + std::to_string(i),
-                       ax_edges.data(), {ax_edges.size()}, mode_str);
+                       ax_edges.data(), {ax_edges.size()}, "a");
       } else {
         axis_types.push_back('d');
         std::vector<int64_t> bins; // we assume that discrete axes have integer bins
@@ -67,14 +74,13 @@ namespace corsika {
         for (int j = 0; j < ax.size(); ++j) { bins.push_back(ax.bin(j).lower()); }
 
         cnpy::npz_save(filename, std::string{"bins_"} + std::to_string(i), bins.data(),
-                       {bins.size()}, mode_str);
+                       {bins.size()}, "a");
       }
     }
 
-    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank},
-                   mode_str);
-    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, mode_str);
-    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, mode_str);
+    cnpy::npz_save(filename, std::string{"axistypes"}, axis_types.data(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"overflow"}, overflow.get(), {rank}, "a");
+    cnpy::npz_save(filename, std::string{"underflow"}, underflow.get(), {rank}, "a");
 
     auto const prod_axis_size = std::accumulate(axes_dims.cbegin(), axes_dims.cend(),
                                                 unsigned{1}, std::multiplies<>());
@@ -98,9 +104,9 @@ namespace corsika {
       temp[p] = *x;
     }
 
-    cnpy::npz_save(filename, "data", temp.get(), axes_dims, mode_str);
+    cnpy::npz_save(filename, "data", temp.get(), axes_dims, "a");
     // In Python this array can directly be assigned to a histogram view if that
     // histogram has its axes correspondingly: hist.view(flow=True)[:] = file['data']
   } // namespace corsika
 
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika

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
@@ -15,38 +14,60 @@
 namespace corsika {
 
   template <typename TDerived>
-  auto const& BaseExponential<TDerived>::getImplementation() const {
+  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>
-  GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
-      setup::Trajectory const& traj, LengthType vL, DirectionVector const& axis) const {
-    if (vL == LengthType::zero()) { return GrammageType::zero(); }
+  inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+      LengthType const height) const {
+    return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+  }
 
-    auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
-    auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+  template <typename TDerived>
+  inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj, DirectionVector const& axis) const {
+    LengthType const length = traj.getLength();
+    if (length == LengthType::zero()) { return GrammageType::zero(); }
+
+    // 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>
-  LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
-      setup::Trajectory const& traj, GrammageType grammage,
+  inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
+      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,12 +75,4 @@ namespace corsika {
     }
   }
 
-  template <typename TDerived>
-  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
@@ -13,33 +12,34 @@
 namespace corsika {
 
   template <typename IEnvironmentModel>
-  Environment<IEnvironmentModel>::Environment()
+  inline Environment<IEnvironmentModel>::Environment()
       : coordinateSystem_{get_root_CoordinateSystem()}
       , universe_(std::make_unique<BaseNodeType>(
             std::make_unique<Universe>(coordinateSystem_))) {}
 
   template <typename IEnvironmentModel>
-  typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr&
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr&
   Environment<IEnvironmentModel>::getUniverse() {
     return universe_;
   }
 
   template <typename IEnvironmentModel>
-  typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr const&
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr const&
   Environment<IEnvironmentModel>::getUniverse() const {
     return universe_;
   }
 
   template <typename IEnvironmentModel>
-  CoordinateSystemPtr const& Environment<IEnvironmentModel>::getCoordinateSystem() const {
+  inline CoordinateSystemPtr const& Environment<IEnvironmentModel>::getCoordinateSystem()
+      const {
     return coordinateSystem_;
   }
 
   // factory method for creation of VolumeTreeNodes
   template <typename IEnvironmentModel>
   template <class TVolumeType, typename... TVolumeArgs>
-  std::unique_ptr<VolumeTreeNode<IEnvironmentModel> >
-  Environment<IEnvironmentModel>::createNode(TVolumeArgs&&... args) {
+  std::unique_ptr<VolumeTreeNode<IEnvironmentModel> > inline Environment<
+      IEnvironmentModel>::createNode(TVolumeArgs&&... args) {
     static_assert(std::is_base_of_v<IVolume, TVolumeType>,
                   "unusable type provided, needs to be derived from "
                   "\"Volume\"");
@@ -48,4 +48,24 @@ namespace corsika {
         std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
   }
 
-} // namespace corsika
\ No newline at end of file
+  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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
+      double const n0, InverseLengthType const lambda, Point const center,
+      LengthType const radius, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , n0_(n0)
+      , lambda_(lambda)
+      , center_(center)
+      , radius_(radius) {}
+
+  template <typename T>
+  inline double ExponentialRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
+  }
+
+} // namespace corsika

corsika/detail/media/FlatExponential.inl

 /*
  * (c) Copyright 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
@@ -17,36 +16,35 @@
 namespace corsika {
 
   template <typename T>
-  FlatExponential<T>::FlatExponential(Point const& point,
-                                      Vector<dimensionless_d> const& axis,
-                                      MassDensityType rho, LengthType lambda,
-                                      NuclearComposition const& nuclComp)
-      : BaseExponential<FlatExponential<T>>(point, rho, lambda)
+  inline FlatExponential<T>::FlatExponential(Point const& point,
+                                             DirectionVector const& axis,
+                                             MassDensityType const rho,
+                                             LengthType const lambda,
+                                             NuclearComposition const& nuclComp)
+      : BaseExponential<FlatExponential<T>>(point, 0_m, rho, lambda)
       , axis_(axis)
       , nuclComp_(nuclComp) {}
 
   template <typename T>
-  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_));
+  inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
+    return BaseExponential<FlatExponential<T>>::getMassDensity(
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
   }
 
   template <typename T>
-  NuclearComposition const& FlatExponential<T>::getNuclearComposition() const {
+  inline NuclearComposition const& FlatExponential<T>::getNuclearComposition() const {
     return nuclComp_;
   }
 
   template <typename T>
-  GrammageType FlatExponential<T>::getIntegratedGrammage(setup::Trajectory const& line,
-                                                         LengthType to) const {
-    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
+  inline GrammageType FlatExponential<T>::getIntegratedGrammage(
+      BaseTrajectory const& line) const {
+    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
   }
 
   template <typename T>
-  LengthType FlatExponential<T>::getArclengthFromGrammage(setup::Trajectory const& line,
-                                                          GrammageType grammage) const {
+  inline LengthType FlatExponential<T>::getArclengthFromGrammage(
+      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
@@ -16,29 +15,30 @@
 namespace corsika {
 
   template <typename T>
-  HomogeneousMedium<T>::HomogeneousMedium(MassDensityType density,
-                                          NuclearComposition const& nuclComp)
+  inline HomogeneousMedium<T>::HomogeneousMedium(MassDensityType density,
+                                                 NuclearComposition const& nuclComp)
       : density_(density)
       , nuclComp_(nuclComp) {}
 
   template <typename T>
-  MassDensityType HomogeneousMedium<T>::getMassDensity(Point const&) const {
+  inline MassDensityType HomogeneousMedium<T>::getMassDensity(Point const&) const {
     return density_;
   }
+
   template <typename T>
-  NuclearComposition const& HomogeneousMedium<T>::getNuclearComposition() const {
+  inline NuclearComposition const& HomogeneousMedium<T>::getNuclearComposition() const {
     return nuclComp_;
   }
 
   template <typename T>
-  GrammageType HomogeneousMedium<T>::getIntegratedGrammage(setup::Trajectory const&,
-                                                           LengthType to) const {
-    return to * density_;
+  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+      BaseTrajectory const& track) const {
+    return track.getLength() * density_;
   }
 
   template <typename T>
-  LengthType HomogeneousMedium<T>::getArclengthFromGrammage(setup::Trajectory const&,
-                                                            GrammageType grammage) const {
+  inline LengthType HomogeneousMedium<T>::getArclengthFromGrammage(
+      BaseTrajectory const&, GrammageType grammage) const {
     return grammage / density_;
   }
 } // namespace corsika

corsika/detail/media/InhomogeneousMedium.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
@@ -17,32 +16,32 @@ namespace corsika {
 
   template <typename T, typename TDensityFunction>
   template <typename... TArgs>
-  InhomogeneousMedium<T, TDensityFunction>::InhomogeneousMedium(
+  inline InhomogeneousMedium<T, TDensityFunction>::InhomogeneousMedium(
       NuclearComposition const& nuclComp, TArgs&&... rhoTArgs)
       : nuclComp_(nuclComp)
       , densityFunction_(rhoTArgs...) {}
 
   template <typename T, typename TDensityFunction>
-  MassDensityType InhomogeneousMedium<T, TDensityFunction>::getMassDensity(
+  inline MassDensityType InhomogeneousMedium<T, TDensityFunction>::getMassDensity(
       Point const& point) const {
     return densityFunction_.evaluateAt(point);
   }
 
   template <typename T, typename TDensityFunction>
-  NuclearComposition const&
+  inline NuclearComposition const&
   InhomogeneousMedium<T, TDensityFunction>::getNuclearComposition() const {
     return nuclComp_;
   }
 
   template <typename T, typename TDensityFunction>
-  GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
-      setup::Trajectory const& line, LengthType to) const {
-    return densityFunction_.getIntegrateGrammage(line, to);
+  inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
+      BaseTrajectory const& line) const {
+    return densityFunction_.getIntegrateGrammage(line);
   }
 
   template <typename T, typename TDensityFunction>
-  LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
-      setup::Trajectory const& line, GrammageType grammage) const {
+  inline LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType grammage) const {
     return densityFunction_.getArclengthFromGrammage(line, grammage);
   }
 

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,15 @@
 #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>
-  void LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
-                                         TModelArgs...>::checkRadius(LengthType r) const {
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::checkRadius(LengthType const r)
+      const {
     if (r <= previousRadius_) {
       throw std::runtime_error("radius must be greater than previous");
     }
@@ -27,7 +28,7 @@ namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  void LayeredSphericalAtmosphereBuilder<
+  inline void LayeredSphericalAtmosphereBuilder<
       TMediumInterface, TMediumModelExtra,
       TModelArgs...>::setNuclearComposition(NuclearComposition const& composition) {
     composition_ = std::make_unique<NuclearComposition>(composition);
@@ -35,26 +36,28 @@ namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  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
@@ -62,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>
-  void LayeredSphericalAtmosphereBuilder<
+  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
@@ -92,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>>(
@@ -104,7 +108,41 @@ namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
+  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<
       TMediumInterface, TMediumModelExtra, TModelArgs...>::assemble() {
     Environment<TMediumInterface> env;
     assemble(env);
@@ -113,7 +151,7 @@ namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  void LayeredSphericalAtmosphereBuilder<
+  inline void LayeredSphericalAtmosphereBuilder<
       TMediumInterface, TMediumModelExtra,
       TModelArgs...>::assemble(Environment<TMediumInterface>& env) {
     auto& universe = env.getUniverse();
@@ -131,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,22 +12,27 @@
 namespace corsika {
 
   template <typename TDerived>
-  auto const& LinearApproximationIntegrator<TDerived>::getImplementation() const {
+  inline TDerived const& LinearApproximationIntegrator<TDerived>::getImplementation()
+      const {
     return *static_cast<TDerived const*>(this);
   }
 
   template <typename TDerived>
-  auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
-      setup::Trajectory 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>
-  auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
-      setup::Trajectory const& line, GrammageType grammage) const {
+  inline LengthType LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType grammage) const {
     auto const c0 = getImplementation().rho_(line.getPosition(0));
     auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
                                                                 line.getDirection(0));
@@ -37,8 +41,8 @@ namespace corsika {
   }
 
   template <typename TDerived>
-  auto LinearApproximationIntegrator<TDerived>::getMaximumLength(
-      setup::Trajectory const& line, [[maybe_unused]] double relError) const {
+  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
@@ -14,17 +13,17 @@ namespace corsika {
 
   template <typename T>
   template <typename... Args>
-  MediumPropertyModel<T>::MediumPropertyModel(Medium const medium, Args&&... args)
+  inline MediumPropertyModel<T>::MediumPropertyModel(Medium const medium, Args&&... args)
       : T(std::forward<Args>(args)...)
       , medium_(medium) {}
 
   template <typename T>
-  Medium MediumPropertyModel<T>::getMedium(Point const&) const {
+  inline Medium MediumPropertyModel<T>::getMedium() const {
     return medium_;
   }
 
   template <typename T>
-  void MediumPropertyModel<T>::setMedium(Medium const medium) {
+  inline void MediumPropertyModel<T>::setMedium(Medium const medium) {
     medium_ = medium;
   }