third_party/phys/units/quantity_io_weber.hpp → corsika/framework/units/quantity_io_weber.hpp

@@ -16,7 +16,7 @@
 #ifndef PHYS_UNITS_QUANTITY_IO_WEBER_HPP_INCLUDED
 #define PHYS_UNITS_QUANTITY_IO_WEBER_HPP_INCLUDED
 
-#include "phys/units/quantity_io.hpp"
+#include "corsika/framework/units/quantity_io.hpp"
 
 namespace phys { namespace units {
 

corsika/framework/utility/COMBoost.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <Eigen/Dense>
+
+namespace corsika {
+
+  /**
+   * @defgroup Utilities
+   *
+   * Collection of classes and methods to perform recurring tasks.
+   */
+
+  /**
+   * @class COMBoost
+   * @ingroup Utilities
+   *
+   * This utility class handles Lorentz boost (in one spatial direction)
+   * between different referenence frames, using FourVector.
+   *
+   * The class is initialized with projectile and optionally target
+   * energy/momentum data. During initialization, a rotation matrix is
+   * calculated to represent the projectile movement (and thus the
+   * boost) along the z-axis. Also the inverse of this rotation is
+   * calculated. The Lorentz boost matrix and its inverse are
+   * determined as 2x2 matrices considering the energy and
+   * pz-momentum.
+   *
+   * Different constructors are offered with different specialization
+   * for the cases of collisions (projectile-target) or just decays
+   * (projectile only).
+   */
+
+  class COMBoost {
+
+  public:
+    /**
+     * Construct a COMBoost given four-vector of projectile and mass of target (target at
+     * rest).
+     *
+     * The FourMomentum and mass define the lab system.
+     */
+    COMBoost(FourMomentum const& P4projectile, HEPEnergyType const massTarget);
+
+    /**
+     * Construct a COMBoost to boost into the rest frame of a particle given its
+     * 3-momentum and mass.
+     */
+    COMBoost(MomentumVector const& momentum, HEPEnergyType const mass);
+
+    /**
+     * Construct a COMBoost given two four-vectors of projectile target.
+     *
+     * The two FourMomentum can define an arbitrary system.
+     */
+    COMBoost(FourMomentum const& P4projectile, FourMomentum const& P4target);
+
+    //! transforms a 4-momentum from lab frame to the center-of-mass frame
+    template <typename FourVector>
+    FourVector toCoM(FourVector const& p4) const;
+
+    //! transforms a 4-momentum from the center-of-mass frame back to lab frame
+    template <typename FourVector>
+    FourVector fromCoM(FourVector const& p4) const;
+
+    //! returns the rotated coordinate system: +z is projectile direction
+    CoordinateSystemPtr const& getRotatedCS() const;
+
+    //! returns the original coordinate system of the projectile (lab)
+    CoordinateSystemPtr const& getOriginalCS() const;
+
+  protected:
+    //! internal method
+    void setBoost(double const coshEta, double const sinhEta);
+
+  public:
+    Eigen::Matrix2d boost_;
+    Eigen::Matrix2d inverseBoost_;
+    CoordinateSystemPtr const originalCS_;
+    CoordinateSystemPtr rotatedCS_;
+  };
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/COMBoost.inl>

corsika/framework/utility/CorsikaData.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/filesystem/path.hpp>
+
+namespace corsika {
+  /**
+   * @file CorsikaData.hpp
+   * @ingroup Utilities
+   * @{
+   * returns the full path of the file \p filename within the CORSIKA_DATA directory.
+   */
+  boost::filesystem::path corsika_data(boost::filesystem::path const& filename);
+
+  //! @}
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/CorsikaData.inl>

corsika/framework/utility/CorsikaFenv.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <cfenv>
+
+/*
+ * Same declaration of function as provided in GLIBC
+ * Repetition allowed in the case where cfenv defines the functions already, no clash.
+ */
+extern "C" {
+
+int feenableexcept(int excepts) noexcept;
+int fedisableexcept(int excepts) noexcept;
+}
+
+#ifdef CORSIKA_HAS_FEENABLEEXCEPT
+// Nothing to do, OS privides the functions
+#else
+#ifdef CORSIKA_OS_MAC
+#include <corsika/detail/framework/utility/CorsikaFenvOSX.inl>
+#else
+#include <corsika/detail/framework/utility/CorsikaFenvFallback.inl>
+#endif
+#endif

corsika/framework/utility/CrossSectionTable.hpp

+/*
+ * (c) Copyright 2024 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 <algorithm>
+#include <cmath>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include <boost/functional/identity.hpp>
+#include <boost/iterator/transform_iterator.hpp>
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+  namespace InterpolationTransforms {
+    using Identity = boost::identity;
+
+    struct Log {
+      template <typename T>
+      auto operator()(T val) const {
+        if constexpr (is_quantity_v<T>) {
+          return std::log(val.magnitude());
+        } else {
+          return std::log(val);
+        }
+      }
+    };
+  } // namespace InterpolationTransforms
+
+  template <typename Transform>
+  class CrossSectionTable {
+  public:
+    template <typename InputIt1, typename InputIt2>
+    CrossSectionTable(InputIt1 energiesFirst, InputIt1 energiesLast,
+                      InputIt2 crosssectionsFirst) {
+      static_assert(std::is_same_v<typename InputIt1::value_type, HEPEnergyType>);
+      static_assert(std::is_same_v<typename InputIt2::value_type, CrossSectionType>);
+
+      table_.reserve(std::distance(energiesFirst, energiesLast));
+
+      auto itE = boost::make_transform_iterator(std::move(energiesFirst), transform_);
+      decltype(itE) const itEEnd =
+          boost::make_transform_iterator(std::move(energiesLast), transform_);
+
+      for (auto itXS = std::move(crosssectionsFirst); itE != itEEnd; ++itE, ++itXS) {
+        table_.emplace_back(*itE, *itXS);
+      }
+
+      std::sort(table_.begin(), table_.end(), less_datapoint);
+    }
+
+    CrossSectionType interpolate(HEPEnergyType energy) const {
+      auto const transformed_val = transform_(energy);
+      auto const lb_it =
+          std::lower_bound(table_.cbegin(), table_.cend(), transformed_val, less_x);
+      if (lb_it == table_.cbegin()) {
+        throw std::runtime_error{"CrossSectionTable: value out of bounds (lower limit)"};
+      }
+      if (lb_it == table_.cend()) {
+        throw std::runtime_error{"CrossSectionTable: value out of bounds (upper limit)"};
+      }
+
+      auto const prev_it = std::prev(lb_it);
+      auto const lambda =
+          (transformed_val - prev_it->first) / (lb_it->first - prev_it->first);
+      return lambda * lb_it->second + (1 - lambda) * prev_it->second;
+    }
+
+  private:
+    using key_type = std::invoke_result_t<Transform, HEPEnergyType>;
+    using datapoint_type = std::pair<key_type, CrossSectionType>;
+    std::vector<datapoint_type> table_;
+    Transform transform_{};
+
+    static bool less_datapoint(datapoint_type const& a, datapoint_type const& b) {
+      return a.first < b.first;
+    }
+    static bool less_x(datapoint_type const& a, typename datapoint_type::first_type b) {
+      return a.first < b;
+    }
+  };
+} // namespace corsika

corsika/framework/utility/CubicSolver.hpp

+/*
+ * (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 <algorithm>
+#include <numeric>
+#include <complex>
+
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+/**
+ *   @file CubicSolver.hpp
+ */
+
+namespace corsika {
+
+  /**
+   *  @ingroup Utility
+   * @{
+   */
+
+  namespace andre {
+    /**
+       Solve a x^3 + b x^2 + c x + d = 0
+    */
+    std::vector<double> solveP3(long double a, long double b, long double c,
+                                long double d, double const epsilon = 1e-12);
+  } // namespace andre
+
+  /**
+       Solve depressed cubic: x^3 + p x + q = 0
+   */
+  inline std::vector<double> solve_cubic_depressed_disciminant_real(
+      long double p, long double q, long double const disc, double const epsilon = 1e-12);
+
+  std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+                                                 double const epsilon = 1e-12);
+
+  /**
+     Solve a x^3 + b x^2 + c x + d = 0
+
+     Analytical approach. Not very stable in all conditions.
+   */
+  std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+                                                long double c, long double d,
+                                                double const epsilon = 1e-12);
+
+  /**
+   * Cubic function.
+   *
+   * T must be a floating point type.
+   *
+   *  @return a x^3 + b x^2 + c x + d
+   */
+  template <typename T>
+  T cubic_function(T x, T a, T b, T c, T d);
+
+  /**
+     @return 3 a x^2 + 2 b x + c
+  */
+  template <typename T>
+  T cubic_function_dfdx(T x, T a, T b, T c);
+
+  /**
+     @return 6 a x + 2 b
+  */
+  template <typename T>
+  T cubic_function_d2fd2x(T x, T a, T b);
+
+  /**
+     Iterative approach to solve: a x^3 + b x^2 + c x + d = 0
+
+     This is often fastest and most precise.
+   */
+  std::vector<double> solve_cubic_real(long double a, long double b, long double c,
+                                       long double d, double const epsilon = 1e-12);
+
+  //! @}
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/CubicSolver.inl>

corsika/framework/utility/FindXmax.hpp

+/*
+ * (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 <corsika/framework/core/Logging.hpp>
+
+#include <tuple>
+
+namespace corsika {
+
+  /**
+   * Interpolates profiles around maximum.
+   *
+   * The maximum of profiles can be robustly estimated from the three maximal points by
+   * quadratic interpolation. This code is copied from CONEX and is awaiting full adaption
+   * into CORSIKA8. This is a temporary solution (just copied).
+   */
+
+  class FindXmax {
+
+    static bool invert3by3(double A[3][3]) {
+      const double kSmall = 1e-80;
+
+      double determinant = A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1]) -
+                           A[0][1] * (A[1][0] * A[2][2] - A[1][2] * A[2][0]) +
+                           A[0][2] * (A[1][0] * A[2][1] - A[1][1] * A[2][0]);
+
+      double absDet = fabs(determinant);
+
+      if (absDet < kSmall) {
+        CORSIKA_LOG_WARN("invert3by3: Error-matrix singular (absDet={})", absDet);
+        return false;
+      }
+
+      double B[3][3];
+
+      B[0][0] = A[1][1] * A[2][2] - A[1][2] * A[2][1];
+      B[1][0] = A[1][2] * A[2][0] - A[2][2] * A[1][0];
+      B[2][0] = A[1][0] * A[2][1] - A[1][1] * A[2][0];
+
+      B[0][1] = A[0][2] * A[2][1] - A[2][2] * A[0][1];
+      B[1][1] = A[0][0] * A[2][2] - A[2][0] * A[0][2];
+      B[2][1] = A[0][1] * A[2][0] - A[0][0] * A[2][1];
+
+      B[0][2] = A[0][1] * A[1][2] - A[1][1] * A[0][2];
+      B[1][2] = A[0][2] * A[1][0] - A[1][2] * A[0][0];
+      B[2][2] = A[0][0] * A[1][1] - A[0][1] * A[1][0];
+
+      for (int i = 0; i < 3; i++)
+        for (int j = 0; j < 3; j++) A[i][j] = B[i][j] / determinant;
+
+      return true;
+    }
+
+    /****************************************************
+     *
+     * solves linear system
+     *
+     *   / y[0] \   / A[0][0] A[0][1] A[0][2] \   / x[0] \
+     *   | y[1] | = | A[1][0] A[1][1] A[1][2] | * | x[1] |
+     *   \ y[2] /   \ A[2][0] A[2][1] A[2][2] /   \ x[2] /
+     *
+     *
+     * Input:  y[3] and A[3][3]
+     * Output: returns true when succeded (i.e. A is not singular)
+     *         A is overwritten with its inverse
+     *
+     * M.Unger 12/1/05
+     *
+     ****************************************************/
+    static bool solve3by3(double y[3], double A[3][3], double x[3]) {
+      if (invert3by3(A)) {
+
+        for (int i = 0; i < 3; i++)
+          x[i] = A[i][0] * y[0] + A[i][1] * y[1] + A[i][2] * y[2];
+
+        return true;
+
+      } else
+        return false;
+    }
+
+  public:
+    static std::tuple<double, double> interpolateProfile(double x1, double x2, double x3,
+                                                         double y1, double y2,
+                                                         double y3) {
+
+      // quadratic "fit" around maximum to get dEdXmax and Xmax
+      double x[3] = {x1, x2, x3};
+      double y[3] = {y1, y2, y3};
+
+      double A[3][3];
+      A[0][0] = x[0] * x[0];
+      A[0][1] = x[0];
+      A[0][2] = 1.;
+      A[1][0] = x[1] * x[1];
+      A[1][1] = x[1];
+      A[1][2] = 1.;
+      A[2][0] = x[2] * x[2];
+      A[2][1] = x[2];
+      A[2][2] = 1.;
+
+      double a[3];
+
+      solve3by3(y, A, a);
+
+      if (a[0] < 0.) {
+        double const Xmax = -a[1] / (2. * a[0]);
+        return std::make_tuple(Xmax, a[0] * Xmax * Xmax + a[1] * Xmax + a[2]);
+      }
+      return std::make_tuple(0, 0);
+    }
+  };
+} // namespace corsika
\ No newline at end of file

corsika/framework/utility/HasMethodSignature.hpp

+#pragma once
+
+/*
+ * (c) Copyright 2018 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.
+ */
+
+/**
+ * @file HasMethodSignature.hpp
+ */
+
+#include <type_traits>
+
+namespace corsika {
+
+  namespace detail {
+
+    /**
+       Helper traits class (partial) for static compile time checking.
+
+       Note, this is a poor replacement for C++20 concepts... they are
+       eagerly awaited!
+
+       It defines the default body of a generic test function returning
+       std::false_type.
+
+       In addition it defines the pattern for class-method matching with a
+       return type TReturn and function arguments TArgs... . Right now
+       both method signatures, "const" and "not const", are matched.
+    */
+    template <typename TReturn, typename... TArgs>
+    struct has_method_signature {
+
+      // the non-const version
+      template <class T>
+      static std::true_type testSignature(TReturn (T::*)(TArgs...));
+
+      // the const version
+      template <class T>
+      static std::true_type testSignature(TReturn (T::*)(TArgs...) const);
+    };
+
+  } // namespace detail
+} // namespace corsika

corsika/framework/utility/ImplementsMixin.hpp

+/*
+ * (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
+
+/**
+ * @file ImplementsMixin.hpp
+ */
+
+#include <type_traits>
+
+namespace corsika {
+
+  namespace detail {
+
+    /**
+       Helper traits class (partial) for static compile time checking.
+
+       This method checks whether a given class implements a particular
+       mixin - this is particularly useful in the media heirarchy that utilizes
+       mixins for the interface definition.
+    */
+    template <template <typename> typename Mixin, typename T>
+    struct implements_mixin {
+
+      template <typename U>
+      static std::true_type test(Mixin<U>&);
+      static std::false_type test(...);
+
+    public:
+      using type = decltype(test(std::declval<T&>()));
+      static constexpr bool value = type::value;
+    };
+
+    template <template <typename> typename Mixin, typename T>
+    bool constexpr implements_mixin_v = implements_mixin<Mixin, T>::value;
+
+  } // namespace detail
+} // namespace corsika

corsika/framework/utility/LinearSolver.hpp

+/*
+ * (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 {
+
+  std::vector<double> solve_linear_real(double a, double b);
+
+  std::vector<std::complex<double>> solve_linear(double a, double b);
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/LinearSolver.inl>

corsika/framework/utility/QuadraticSolver.hpp

+/*
+ * (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>
+
+#include <corsika/framework/utility/LinearSolver.hpp>
+
+namespace corsika {
+
+  std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+                                                    long double c,
+                                                    double const epsilon = 1e-12);
+
+  std::vector<double> solve_quadratic_real(long double a, long double b, long double c,
+                                           double const epsilon = 1e-12);
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/QuadraticSolver.inl>

corsika/framework/utility/QuarticSolver.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+
+#include <complex>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+
+/**
+ * @file QuarticSolver.hpp
+ */
+
+namespace corsika {
+
+  /**
+   * @ingroup Utilities
+   * @{
+   */
+
+  namespace andre {
+
+    /**
+     solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e
+     Attention - this function returns dynamically allocated array. It has to be
+     released afterwards.
+    */
+    std::vector<double> solve_quartic_real(long double a, long double b, long double c,
+                                           long double d, long double e,
+                                           double const epsilon = 1e-12);
+  } // namespace andre
+
+  /**
+     solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e
+  */
+  std::vector<double> solve_quartic_real(long double a, long double b, long double c,
+                                         long double d, long double e,
+                                         double const epsilon = 1e-12);
+
+  //! @}
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/QuarticSolver.inl>

corsika/framework/utility/SaveBoostHistogram.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <boost/histogram.hpp>
+
+namespace corsika {
+
+  /**
+   * @file SaveBoostHistogram.hpp
+   * @ingroup Utilities
+   *
+   * This functions saves a boost::histogram into a numpy file. Only rather basic axis
+   * types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
+   * counts (i.e. a double or int) are supported, nothing fancy like profiles.
+   *
+   * Note that this function makes a temporary, dense copy of the histogram, which could
+   * be an issue for huge sizes (e.g. for high dimensions)
+   *
+   * @param overwrite silently overwrite existing files if true, otherwise throw
+   * runtime_error
+   */
+  template <class Axes, class Storage>
+  void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
+                 std::string const& filename, bool overwrite = true);
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/SaveBoostHistogram.inl>

corsika/framework/utility/Singleton.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+  /**
+   * \class Singleton Singleton.h utl/Singleton.h
+   *
+   * \brief Curiously Recurring Template Pattern (CRTP) for Meyers singleton
+   *
+   * The singleton class is implemented as follows
+   * \code
+   * #include <utl/Singleton.hpp>
+   *
+   * class SomeClass : public utl::Singleton<SomeClass> {
+   *   ...
+   * private:
+   *   // prevent creation, destruction
+   *   SomeClass() { }
+   *   ~SomeClass() { }
+   *
+   *   friend class utl::Singleton<SomeClass>;
+   * };
+   * \endcode
+   * Singleton automatically prevents copying of the derived class.
+   *
+   * \author Darko Veberic
+   * \date 9 Aug 2006
+   * \version $Id: Singleton.h 25091 2014-01-30 09:49:57Z darko $
+   * \ingroup stl
+   */
+
+  template <typename T>
+  class Singleton {
+  public:
+    static T& getInstance() {
+      static T instance;
+      return instance;
+    }
+
+    Singleton(const Singleton&) =
+        delete; // Singleton Classes should not be copied. Removes move constructor and
+                // move assignment as well
+    Singleton& operator=(const Singleton&) =
+        delete; // Singleton Classes should not be copied.
+
+  protected:
+    // derived class can call ctor and dtor
+    Singleton() {}
+    ~Singleton() {}
+  };
+
+} // namespace corsika

corsika/media/BaseExponential.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
+#include <limits>
+
+/**
+ * @file corsika/media/BaseExponential.hpp
+ */
+
+namespace corsika {
+
+  /**
+   * This class provides the grammage/length conversion functionality for
+   * (locally) flat exponential atmospheres.
+   *
+   * The density is described according to \f[ \varrho() \f]
+   */
+  template <typename TDerived>
+  class BaseExponential {
+
+  public:
+    BaseExponential(Point const& point, LengthType const referenceHeight,
+                    MassDensityType const rho0, LengthType const lambda);
+
+    Point const& getAnchorPoint() const { return point_; }
+
+  protected:
+    auto const& getImplementation() const;
+
+    /**
+     * Returns the mass density at altitude "height".
+     *
+     * @param height
+     * @return MassDensityType
+     */
+    MassDensityType getMassDensity(LengthType const height) const;
+
+    // clang-format off
+    /**
+     * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
+     * direction \f$ \vec{u} \f$ is given by:
+     * \f[
+     *   X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right) \quad \text{,}
+     * \f]
+     * where \f$ \varrho_0 \f$ is the density at the starting point.
+     *
+     * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
+     * \f[
+     *   X = \varrho_0 l
+     * \f]
+     */
+    // clang-format on
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line,
+                                       DirectionVector const& axis) const;
+
+    // clang-format off
+    /**
+     * For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized)
+     * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by:
+     * \f[
+     *   l = \begin{cases}
+     *   \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} & Y := 1 +
+     *     \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda} > 0 \\
+     *   \infty & \text{else} & \text{,}
+     *   \end{cases}
+     * \f] where \f$ \varrho_0 \f$ is the density at the starting point.
+     *
+     * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like for a homogeneous density:
+     * \f[
+     *   l =  \frac{X}{\varrho_0}
+     * \f]
+     */
+    // clang-format on
+    LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+                                        GrammageType const grammage,
+                                        DirectionVector const& axis) const;
+
+  private:
+    MassDensityType const rho0_;
+    LengthType const lambda_;
+    InverseLengthType const invLambda_;
+    Point const point_;
+    LengthType const referenceHeight_;
+
+  }; // class BaseExponential
+
+} // namespace corsika
+
+#include <corsika/detail/media/BaseExponential.inl>

corsika/media/BaseTabular.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
+
+#include <limits>
+#include <vector>
+
+namespace corsika {
+
+  /**
+   * This class provides the grammage/length conversion functionality for
+   * (locally) flat tabulated atmospheres.
+   */
+  template <typename TDerived>
+  class BaseTabular {
+
+  public:
+    BaseTabular(Point const& point, LengthType const referenceHeight,
+                std::function<MassDensityType(LengthType)> const& rho,
+                unsigned int const nBins, LengthType const deltaHeight);
+
+    Point const& getAnchorPoint() const { return point_; }
+
+  protected:
+    auto const& getImplementation() const;
+
+    MassDensityType getMassDensity(LengthType const height) const;
+
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const;
+
+    LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+                                        GrammageType const grammage) const;
+
+  private:
+    unsigned int const nBins_;
+    LengthType deltaHeight_;
+    std::vector<MassDensityType> density_;
+    Point const point_;
+    LengthType const referenceHeight_;
+
+  }; // class BaseTabular
+
+} // namespace corsika
+
+#include <corsika/detail/media/BaseTabular.inl>

corsika/media/CORSIKA7Atmospheres.hpp

+/*
+ * (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 <corsika/media/IRefractiveIndexModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/framework/utility/ImplementsMixin.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+// for detail namespace, NoExtraModelInner, NoExtraModel and traits
+#include <corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp>
+
+namespace corsika {
+
+  /**
+   * Atmosphere Ids following the CORSIKA 7 5-layered atmosphere models.
+   *
+   * Each model corresponds to a standard 5-layered atmosphere model. 4 Layers are
+   * exponential, the outer layer is with constant density.
+   *
+   * All atmospheres are valid for heights (above Earth sea level) up to 112.8 km.
+   */
+  enum class AtmosphereId : uint8_t {
+    LinsleyUSStd = 0,
+    MiddleEuropeJan,
+    MiddleEuropeFeb,
+    MiddleEuropeMay,
+    MiddleEuropeJun,
+    MiddleEuropeAug,
+    MiddleEuropeOct,
+    MiddleEuropeDec,
+    SouthPoleMar,
+    SouthPoleJul,
+    SouthPoleOct,
+    SouthPoleDec,
+    SouthPoleJan,
+    SouthPoleAug,
+    MalargueWinterI,
+    MalargueWinterII,
+    MalargueSpring,
+    MalargueSummer,
+    MalargueAutumn,
+    USStdBK,
+    LastAtmosphere
+  };
+
+  namespace {
+    /**
+     * Struct to hold parameters of one layer of a CORSIKA7 atmosphere.
+     *
+     * The definition of each layer is according to a BaseExponential:
+     * @f[
+     *   \varrho = offset/scaleHeight \cdot
+     * \exp\left(-(height-altitude)/scaleHeight\right)
+     * @f],
+     * where @f$ altitude @f$ is the height where the atmosphere layer starts, @f$
+     * offset/scaleHeight
+     * @f$ is the density at this height.
+     */
+    struct AtmosphereLayerParameters {
+      LengthType altitude;
+      GrammageType offset;
+      LengthType scaleHeight;
+    };
+
+    /**
+     * All the 5 layers of a CORSIKA7 atmosphere.
+     */
+    typedef std::array<AtmosphereLayerParameters, 5> AtmosphereParameters;
+
+    /**
+     * Local, internal helper function to provide "Grammage" type.
+     *
+     * @param v
+     * @return value v with g/cm2 units
+     */
+    auto constexpr grammage(double const v) { return v * 1_g / (1_cm * 1_cm); }
+
+    std::array<AtmosphereParameters,
+               static_cast<uint8_t>(
+                   AtmosphereId::LastAtmosphere)> constexpr atmosphereParameterList{
+        {{{{4_km, grammage(1222.6562), 994186.38_cm},
+           {10_km, grammage(1144.9069), 878153.55_cm},
+           {40_km, grammage(1305.5948), 636143.04_cm},
+           {100_km, grammage(540.1778), 772170.16_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1173.9861), 919546_cm},
+           {10_km, grammage(1205.7625), 963267.92_cm},
+           {40_km, grammage(1386.7807), 614315_cm},
+           {100_km, grammage(555.8935), 739059.6_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1240.48), 933697_cm},
+           {10_km, grammage(1117.85), 765229_cm},
+           {40_km, grammage(1210.9), 636790_cm},
+           {100_km, grammage(608.2128), 733793.8_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1285.2782), 1088310_cm},
+           {10_km, grammage(1173.1616), 935485_cm},
+           {40_km, grammage(1320.4561), 635137_cm},
+           {100_km, grammage(680.6803), 727312.6_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1251.474), 1032310_cm},
+           {10_km, grammage(1173.321), 925528_cm},
+           {40_km, grammage(1307.826), 645330_cm},
+           {100_km, grammage(763.1139), 720851.4_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(113.3362), 923077_cm},
+           {10_km, grammage(1226.5761), 1109960_cm},
+           {40_km, grammage(1382.6933), 630217_cm},
+           {100_km, grammage(685.6073), 726901.3_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1262.7013), 1059360_cm},
+           {10_km, grammage(1139.0249), 888814_cm},
+           {10_km, grammage(1270.2886), 639902_cm},
+           {40_km, grammage(681.4061), 727251.8_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1210.4), 970276_cm},
+           {10_km, grammage(1103.8629), 820946_cm},
+           {40_km, grammage(1215.3545), 639074_cm},
+           {100_km, grammage(629.7611), 731776.5_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{4_km, grammage(1130.74), 867358_cm},
+           {10_km, grammage(1052.05), 741208_cm},
+           {40_km, grammage(1137.21), 633846_cm},
+           {100_km, grammage(442.512), 759850_cm},
+           {112.8_km, grammage(1), 5.4303203e9_cm}}},
+         {{{4_km, grammage(1183.70), 875221_cm},
+           {10_km, grammage(1108.06), 753212_cm},
+           {40_km, grammage(1424.02), 545846_cm},
+           {100_km, grammage(207.595), 793043_cm},
+           {112.8_km, grammage(1), 5.9787908e9_cm}}},
+         {{{4_km, grammage(1177.19), 861745_cm},
+           {10_km, grammage(1125.11), 765925_cm},
+           {40_km, grammage(1304.77), 581351_cm},
+           {100_km, grammage(433.823), 775155_cm},
+           {112.8_km, grammage(1), 7.4095699e9_cm}}},
+         {{{4_km, grammage(1139.99), 861913_cm},
+           {10_km, grammage(1073.82), 744955_cm},
+           {40_km, grammage(1052.96), 675928_cm},
+           {100_km, grammage(492.503), 829627_cm},
+           {112.8_km, grammage(1), 5.8587010e9_cm}}},
+         {{{2.67_km, grammage(1133.10), 861730_cm},
+           {5.33_km, grammage(1101.20), 826340_cm},
+           {8_km, grammage(1085.00), 790950_cm},
+           {100_km, grammage(1098.00), 682800_cm},
+           {112.8_km, grammage(1), 26798156e9_cm}}},
+         {{{6.67_km, grammage(1079.00), 764170_cm},
+           {13.33_km, grammage(1071.90), 699910_cm},
+           {20_km, grammage(1182.00), 635650_cm},
+           {100_km, grammage(1647.10), 551010_cm},
+           {112.8_km, grammage(1), 59.329575e9_cm}}},
+         {{{8_km, grammage(1198.5972), 945766.30_cm},
+           {18.1_km, grammage(1198.8796), 681780.12_cm},
+           {34.5_km, grammage(1419.4152), 620224.52_cm},
+           {100_km, grammage(730.6380), 728157.92_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{8.3_km, grammage(1179.5010), 939228.66_cm},
+           {12.9_km, grammage(1172.4883), 787969.34_cm},
+           {34_km, grammage(1437.4911), 620008.53_cm},
+           {100_km, grammage(761.3281), 724585.33_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{5.9_km, grammage(1202.8804), 977139.52_cm},
+           {12.0_km, grammage(1148.6275), 858087.01_cm},
+           {34.5_km, grammage(1432.0312), 614451.60_cm},
+           {100_km, grammage(696.42788), 730875.73_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{9_km, grammage(1175.3347), 986169.72_cm},
+           {14.6_km, grammage(1180.3694), 793171.45_cm},
+           {33_km, grammage(1614.5404), 600120.95_cm},
+           {100_km, grammage(755.56438), 725247.87_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{8_km, grammage(1196.9290), 985241.1_cm},
+           {13_km, grammage(1173.2537), 819245.00_cm},
+           {33.5_km, grammage(1502.1837), 611220.86_cm},
+           {100_km, grammage(750.89704), 725797.06_cm},
+           {112.8_km, grammage(1), 1e9_cm}}},
+         {{{7_km, grammage(1183.6071), 954248.34_cm},
+           {11.4_km, grammage(1143.0425), 800005.34_cm},
+           {37_km, grammage(1322.9748), 629568.93_cm},
+           {100_km, grammage(655.67307), 737521.77_cm},
+           {112.8_km, grammage(1), 1e9_cm}}}}};
+  } // namespace
+
+  /**
+   * Function to create a CORSIKA 7 5-layer atmosphere.
+   *
+   * @tparam TEnvironmentInterface
+   * @tparam TExtraEnv
+   * @tparam TEnvironment
+   * @tparam TArgs
+   * @param env
+   * @param atmId
+   * @param center
+   * @param args
+   */
+  template <typename TEnvironmentInterface,
+            template <typename> typename TExtraEnv = detail::NoExtraModel,
+            typename TEnvironment, typename... TArgs>
+  void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
+                                Point const& center, TArgs... args);
+
+  //! The standard/default air composition with fraction values based on CORSIKA 7
+  //! Composition (N2,O2,Ar) = (78.084, 20.946, 0.934)
+  //! Pfraction(Ar) = Ar/(2*N2 + 2*O2 + Ar) = 0.00469
+  static inline NuclearComposition const standardAirComposition{
+      {Code::Nitrogen, Code::Oxygen, Code::Argon}, {0.78479, .21052, 0.00469}};
+
+} // namespace corsika
+
+#include <corsika/detail/media/CORSIKA7Atmospheres.inl>

corsika/media/DensityFunction.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/LinearApproximationIntegrator.hpp>
+
+namespace corsika {
+
+  template <class TDerivableRho,
+            template <typename> class TIntegrator = LinearApproximationIntegrator>
+  class DensityFunction
+      : public TIntegrator<DensityFunction<TDerivableRho, TIntegrator>> {
+    friend class TIntegrator<DensityFunction<TDerivableRho, TIntegrator>>;
+
+    TDerivableRho rho_; //!< functor for density
+
+  public:
+    DensityFunction(TDerivableRho rho)
+        : rho_(rho) {}
+
+    MassDensityType evaluateAt(Point const& p) const { return rho_(p); }
+  };
+
+} // namespace corsika

corsika/media/Environment.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/VolumeTreeNode.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+
+#include <corsika/media/Universe.hpp>
+
+#include <limits>
+#include <set>
+
+namespace corsika {
+
+  // fwd decl
+  template <typename IEnvironmentModel>
+  class Environment;
+
+  template <typename IEnvironmentModel>
+  std::set<Code> const get_all_elements_in_universe(
+      Environment<IEnvironmentModel> const& env);
+
+  /**
+   * Base Environment class.
+   *
+   * Describes the Environment in which the shower is propagated.
+   */
+  template <typename IEnvironmentModel>
+  class Environment {
+  public:
+    using BaseNodeType = VolumeTreeNode<IEnvironmentModel>;
+
+    Environment();
+
+    /**
+     * Getters for the universe stored in the Environment.
+     *
+     * @retval Retuns reference to a Universe object with infinite size.
+     */
+    ///@{
+    //! Get non const universe
+    typename BaseNodeType::VTNUPtr& getUniverse();
+    //! Get const universe
+    typename BaseNodeType::VTNUPtr const& getUniverse() const;
+    ///@}
+
+    /**
+     * Getter for the CoordinateSystem used in the Environment.
+     *
+     * @retval Retuns a const reference to the CoordinateSystem used.
+     */
+    CoordinateSystemPtr const& getCoordinateSystem() const;
+
+    /**
+     * Factory method for creation of VolumeTreeNodes.
+     *
+     * @tparam TVolumeType Type of volume to be created
+     * @tparam TVolumeArgs Types to forward to the constructor
+     * @param args Parameter forwarded to the constructor of TVolumeType
+     * @retval Returns unique pointer to a VolumeTreeNode with the same EnvitonmentModel
+     * as this class.
+     */
+    template <class TVolumeType, typename... TVolumeArgs>
+    static std::unique_ptr<BaseNodeType> createNode(TVolumeArgs&&... args);
+
+  private:
+    CoordinateSystemPtr const coordinateSystem_;
+    typename BaseNodeType::VTNUPtr universe_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/media/Environment.inl>

corsika/media/ExponentialRefractiveIndex.hpp

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  /**
+   * An exponential refractive index.
+   *
+   * This class returns the value of an exponential refractive index
+   * for all evaluated locations.
+   *
+   */
+  template <typename T>
+  class ExponentialRefractiveIndex : public T {
+
+    double n0_;                ///< n0 constant.
+    InverseLengthType lambda_; ///< lambda parameter.
+    Point center_;             ///< center of the planet.
+    LengthType radius_;        ///< the planet radius.
+
+  public:
+    /**
+     * Construct an ExponentialRefractiveIndex.
+     *
+     * This is initialized with two parameters n_0 and lambda
+     * and returns the value of the exponential refractive index
+     * at a given point location.
+     *
+     * @param field    The refractive index to return to a given point.
+     */
+    template <typename... Args>
+    ExponentialRefractiveIndex(double const n0, InverseLengthType const lambda,
+                               Point const center, LengthType const radius,
+                               Args&&... args);
+
+    /**
+     * Evaluate the refractive index at a given location using its z-coordinate.
+     *
+     * @param  point    The location to evaluate at.
+     * @returns    The refractive index at this point.
+     */
+    double getRefractiveIndex(Point const& point) const final override;
+
+  }; // END: class ExponentialRefractiveIndex
+
+} // namespace corsika
+
+#include <corsika/detail/media/ExponentialRefractiveIndex.inl>