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 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/utility/CubicSolver.hpp>
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+
 #include <complex>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
 
 /**
- * \todo convert to class
+ * @file QuarticSolver.hpp
  */
 
-namespace corsika::quartic_solver {
-
-  const double PI = 3.141592653589793238463L;
-  const double M_2PI = 2 * PI;
-  const double eps = 1e-12;
-
-  typedef std::complex<double> DComplex;
-
-  //---------------------------------------------------------------------------
-  // useful for testing
-  inline DComplex polinom_2(DComplex x, double a, double b) {
-    // Horner's scheme for x*x + a*x + b
-    return x * (x + a) + b;
-  }
-
-  //---------------------------------------------------------------------------
-  // useful for testing
-  inline DComplex polinom_3(DComplex x, double a, double b, double c) {
-    // Horner's scheme for x*x*x + a*x*x + b*x + c;
-    return x * (x * (x + a) + b) + c;
-  }
-
-  //---------------------------------------------------------------------------
-  // useful for testing
-  inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
-    // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
-    return x * (x * (x * (x + a) + b) + c) + d;
-  }
-
-  //---------------------------------------------------------------------------
-  // 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
-  unsigned int solveP3(double* x, double a, double b, double c);
-
-  //---------------------------------------------------------------------------
-  // 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.
-  DComplex* solve_quartic(double a, double b, double c, double d);
-
-} // namespace corsika::quartic_solver
+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 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,16 +12,22 @@
 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>
-  inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
-                        std::string const& filename, bool overwrite = true);
+  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 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
 
-//#define OFFLINE_USE_GAMMA_SINGLETON
-
 namespace corsika {
   /**
    * \class Singleton Singleton.h utl/Singleton.h

corsika/framework/utility/detail/COMBoost.inl

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/geometry/CoordinateSystem.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/geometry/FourVector.hpp>
-#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/utility/sgn.hpp>
-#include <Eigen/Dense>
-#include "../../core/PhysicalUnits.hpp"
-
-namespace corsika::utl {
-
-  auto const& COMBoost::GetRotationMatrix() const { return fRotation; }
-
-  //! transforms a 4-momentum from lab frame to the center-of-mass frame
-  template <typename FourVector>
-  FourVector COMBoost::toCoM(const FourVector& p) const {
-    using namespace corsika::units::si;
-    auto pComponents = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
-    Eigen::Vector3d eVecRotated = pComponents.eVector;
-    Eigen::Vector2d lab;
-
-    lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
-        (eVecRotated(2) * (1 / 1_GeV).magnitude());
-
-    auto const boostedZ = boost_ * lab;
-    auto const E_CoM = boostedZ(0) * 1_GeV;
-
-    eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
-
-    return FourVector(E_CoM,
-                      corsika::geometry::Vector<hepmomentum_d>(rotatedCS_, eVecRotated));
-  }
-
-  //! transforms a 4-momentum from the center-of-mass frame back to lab frame
-  template <typename FourVector>
-  FourVector COMBoost::fromCoM(const FourVector& p) const {
-    using namespace corsika::units::si;
-    auto pCM = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
-    auto const Ecm = p.GetTimeLikeComponent();
-
-    Eigen::Vector2d com;
-    com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
-
-    C8LOG_TRACE(
-        "COMBoost::fromCoM Ecm={} GeV"
-        " pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
-        Ecm / 1_GeV, pCM / 1_GeV, pCM.norm() / 1_GeV, p.GetNorm() / 1_GeV);
-
-    auto const boostedZ = inverseBoost_ * com;
-    auto const E_lab = boostedZ(0) * 1_GeV;
-
-    pCM.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
-
-    geometry::Vector<typename decltype(pCM)::dimension> pLab{rotatedCS_, pCM};
-    pLab.rebase(originalCS_);
-
-    FourVector f(E_lab, pLab);
-
-    C8LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
-                " plab={} GeV (norm={} GeV) "
-                " GeV), invariant mass = {}",
-                E_lab / 1_GeV, f.GetNorm() / 1_GeV, pLab.GetComponents(),
-                pLab.norm() / 1_GeV);
-
-    return f;
-  }
-
-  COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
-                     const HEPMassType massTarget)
-      : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
-    auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
-    auto const pProjNorm = pProjectile.norm();
-    auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
-    auto const a1 = a(0), a2 = a(1);
-
-    auto const s = sgn(a(2));
-    auto const c = 1 / (1 + s * a(2));
-
-    Eigen::Matrix3d A, B;
-
-    if (s > 0) {
-      A << 1, 0, -a1,                     // comment to prevent clang-format
-          0, 1, -a2,                      // .
-          a1, a2, 1;                      // .
-      B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-          -a1 * a2 * c, -a2 * a2 * c, 0,  // .
-          0, 0, -(a1 * a1 + a2 * a2) * c; // .
-
-    } else {
-      A << 1, 0, a1,                      // comment to prevent clang-format
-          0, -1, -a2,                     // .
-          a1, a2, -1;                     // .
-      B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-          +a1 * a2 * c, +a2 * a2 * c, 0,  // .
-          0, 0, (a1 * a1 + a2 * a2) * c;  // .
-    }
-
-    fRotation = A + B;
-
-    // calculate boost
-    double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
-
-    /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
-    double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
-    //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
-    double const sinhEta = -beta * coshEta;
-
-    std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
-    std::cout << "  det = " << fRotation.determinant() - 1 << std::endl;
-
-    fBoost << coshEta, sinhEta, sinhEta, coshEta;
-
-    fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
-  }
-} // namespace corsika::utl

corsika/framework/utility/detail/CorsikaFenvFallback.inl

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/framework/utility/CorsikaFenv.hpp>
-#include <cfenv>
-
-extern "C" {
-#warning No enabling/disabling of floating point exceptions - platform needs better implementation
-
-int feenableexcept(int excepts) { return -1; }
-
-int fedisableexcept(int excepts) { return -1; }
-}

corsika/media/BaseExponential.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
@@ -12,71 +11,87 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Point.hpp>
-#include <corsika/setup/SetupTrajectory.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
+     * 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)
-     * \f], where \f$ \varrho_0 \f$ is the density at the starting point.
+     *   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;
+     *   X = \varrho_0 l
      * \f]
      */
     // clang-format on
-    GrammageType getIntegratedGrammage(setup::Trajectory const& line, LengthType vL,
+    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
+     * 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 :=  0 > 1 +
-     *     \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda}
-     *   \infty & \text{else,}
+     *   \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 with a homogeneous density:
+     * 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(setup::Trajectory const& line,
-                                        GrammageType grammage,
+    LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+                                        GrammageType const grammage,
                                         DirectionVector const& axis) const;
 
-  public:
-    BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda);
-
-    Point const& getAnchorPoint() const { return point_; }
-    MassDensityType getRho0() const { return rho0_; }
-    InverseLengthType getInvLambda() const { return invLambda_; }
-
   private:
     MassDensityType const rho0_;
     LengthType const lambda_;
     InverseLengthType const invLambda_;
     Point const point_;
+    LengthType const referenceHeight_;
 
   }; // class BaseExponential
 

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 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/media/Environment.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
@@ -17,12 +16,23 @@
 #include <corsika/media/Universe.hpp>
 
 #include <limits>
+#include <set>
 
 namespace corsika {
 
-  /** Base Evnironment class
-   *  Describes the Environment in which the shower is propagated
-   **/
+  // 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:
@@ -30,30 +40,34 @@ namespace corsika {
 
     Environment();
 
-    /** Getters for the universe stored in the Environment
+    /**
+     * Getters for the universe stored in the Environment.
      *
-     * @retval Retuns reference to a Universe object with infinite size
-     **/
+     * @retval Retuns reference to a Universe object with infinite size.
+     */
     ///@{
-    //* Get non const universe */
+    //! Get non const universe
     typename BaseNodeType::VTNUPtr& getUniverse();
-    //* Get const universe */
+    //! Get const universe
     typename BaseNodeType::VTNUPtr const& getUniverse() const;
     ///@}
 
-    /** Getter for the CoordinateSystem used in the Environment
+    /**
+     * Getter for the CoordinateSystem used in the Environment.
      *
-     * @retval Retuns a const reference to the CoordinateSystem used
-     **/
+     * @retval Retuns a const reference to the CoordinateSystem used.
+     */
     CoordinateSystemPtr const& getCoordinateSystem() const;
 
-    /** Factory method for creation of VolumeTreeNodes
+    /**
+     * 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 Retuns unique pointer to a VolumeTreeNode with the same EnvitonmentModel as
-     *this class
-     **/
+     * @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);
 

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>

corsika/media/FlatExponential.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
@@ -13,7 +12,7 @@
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/media/BaseExponential.hpp>
 #include <corsika/media/NuclearComposition.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
 
 namespace corsika {
 
@@ -21,11 +20,14 @@ namespace corsika {
   /**
    * flat exponential density distribution with
    * \f[
-   *  \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
+   *  \varrho(\vec{r}) = \varrho_0 \exp\left( \frac{1}{\lambda} (\vec{r} - \vec{p}) \cdot
    *    \vec{a} \right).
    * \f]
    * \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
-   * with the scale parameter \f$ \lambda \f$.
+   * with the scale parameter \f$ \lambda \f$, \f$ \vec{r} \f$ is the location of 
+   * the evaluation, \f$ \vec{p} \f$ is the anchor point at which \f$ \varrho_0 \f$ 
+   * is given. Thus, the unit vector \f$ \vec{a} \f$ specifies the direction of
+   * <em>decreasing</em> <b>height/altitude</b>.
    */
   // clang-format on
   template <typename T>
@@ -34,18 +36,17 @@ namespace corsika {
 
   public:
     FlatExponential(Point const& point, Vector<dimensionless_d> const& axis,
-                    MassDensityType rho, LengthType lambda,
+                    MassDensityType const rho, LengthType const lambda,
                     NuclearComposition const& nuclComp);
 
     MassDensityType getMassDensity(Point const& point) const override;
 
     NuclearComposition const& getNuclearComposition() const override;
 
-    GrammageType getIntegratedGrammage(setup::Trajectory const& line,
-                                       LengthType to) const override;
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
 
-    LengthType getArclengthFromGrammage(setup::Trajectory const& line,
-                                        GrammageType grammage) const override;
+    LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+                                        GrammageType const grammage) const override;
 
   private:
     DirectionVector const axis_;

corsika/media/GeomagneticModel.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/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <boost/filesystem.hpp>
+#include <fstream>
+#include <string>
+
+namespace corsika {
+
+  /**
+   * A magnetic field calculated with the WMM or IGRF model.
+   * WMM: https://geomag.bgs.ac.uk/documents/WMM2020_Report.pdf
+   * IGRF: https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
+   */
+
+  class GeomagneticModel {
+
+    /**
+     * Internal data structure for a single shell of the spherical harmonic
+     * parameterization.
+     */
+    struct ParameterLine {
+      int n;
+      int m;
+      double g;
+      double h;
+      double dg; //< time dependence of g in "per year"
+      double dh; //< time dependence of h in "per year"
+    };
+
+  public:
+    /**
+     * Construct a new World Magnetic Model object.
+     *
+     * @param center Center of Earth.
+     * @param data Data table to read.
+     */
+    GeomagneticModel(Point const& center, boost::filesystem::path const path =
+                                              corsika::corsika_data("GeoMag/WMM.COF"));
+
+    /**
+     * Calculates the value of the magnetic field.
+     *
+     * @param  year        Year of the evaluation, between 2020 and 2025.
+     * @param  altitude    Height of the location to evaluate the field at,
+     *                    in km between -1 and 850.
+     * @param  latitude    Latitude of the location to evaluate the field at,
+     *                   in degrees between -90 and 90 (negative for southern hemisphere).
+     * @param  longitute   Longitude of the location to evaluate the field at,
+     *                      in degrees between -180 and 180 (negative for western
+     *   hemisphere).
+     *
+     * @returns    The magnetic field vector in nT.
+     */
+    MagneticFieldVector getField(double const year, LengthType const altitude,
+                                 double const latitude, double const longitude);
+
+  private:
+    Point center_;                                         //< Center of Earth.
+    std::map<int, std::vector<ParameterLine>> parameters_; //< epoch to Parameters map
+  };
+} // namespace corsika
+
+#include <corsika/detail/media/GeomagneticModel.inl>
\ No newline at end of file

corsika/media/GladstoneDaleRefractiveIndex.hpp

+/*
+ * (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 {
+
+  /**
+   * A tabulated refractive index.
+   *
+   * This class returns the value of refractive index
+   * for a specific height bin.
+   */
+  template <typename T>
+  class GladstoneDaleRefractiveIndex : public T {
+
+    double const
+        referenceRefractivity_; ///< The constant reference refractive index minus one.
+    InverseMassDensityType const
+        referenceInvDensity_; ///< The constant inverse mass density at reference point.
+
+  public:
+    /**
+     * Construct a GladstoneDaleRefractiveIndex.
+     *
+     * This is initialized with a fixed refractive index
+     * at sea level and uses the Gladstone-Dale law to
+     * calculate the refractive index at a given point.
+     *
+     * @param seaLevelRefractiveIndex  The refractive index at sea level.
+     * @param point AA point at earth's surface.
+     */
+    template <typename... Args>
+    GladstoneDaleRefractiveIndex(double const seaLevelRefractiveIndex, Point const point,
+                                 Args&&... args);
+
+    /**
+     * Evaluate the refractive index at a given location.
+     *
+     * @param  point    The location to evaluate at.
+     * @returns    The refractive index at this point.
+     */
+    double getRefractiveIndex(Point const& point) const override;
+
+  }; // END: class GladstoneDaleRefractiveIndex
+
+} // namespace corsika
+
+#include <corsika/detail/media/GladstoneDaleRefractiveIndex.inl>

corsika/media/HomogeneousMedium.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
@@ -12,15 +11,14 @@
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/media/NuclearComposition.hpp>
-
-#include <corsika/setup/SetupTrajectory.hpp>
-
-/**
- * a homogeneous medium
- */
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
 
 namespace corsika {
 
+  /**
+   * a homogeneous medium
+   */
+
   template <typename T>
   class HomogeneousMedium : public T {
 
@@ -31,10 +29,9 @@ namespace corsika {
 
     NuclearComposition const& getNuclearComposition() const override;
 
-    GrammageType getIntegratedGrammage(setup::Trajectory const&,
-                                       LengthType to) const override;
+    GrammageType getIntegratedGrammage(BaseTrajectory const&) const override;
 
-    LengthType getArclengthFromGrammage(setup::Trajectory const&,
+    LengthType getArclengthFromGrammage(BaseTrajectory const&,
                                         GrammageType grammage) const override;
 
   private:

corsika/media/IEmpty.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
 
 #include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
 
 namespace corsika {
 
@@ -25,18 +24,23 @@ namespace corsika {
 
   class IEmpty {
   public:
-    virtual LengthType getArclengthFromGrammage(setup::Trajectory const&,
+    virtual LengthType getArclengthFromGrammage(BaseTrajectory const&,
                                                 GrammageType) const = 0;
 
+    virtual NuclearComposition const& getNuclearComposition() const = 0;
+
     virtual ~IEmpty() {}
   };
 
   template <typename TModel = IEmpty>
   class Empty : public TModel {
   public:
-    LengthType getArclengthFromGrammage(setup::Trajectory const&, GrammageType) const {
+    LengthType getArclengthFromGrammage(BaseTrajectory const&, GrammageType) const {
       return 0. * meter;
     }
+    NuclearComposition const& getNuclearComposition() const {
+      return NuclearComposition(std::vector<Code>{}, {});
+    };
   };
 
 } // namespace corsika

corsika/media/IMagneticFieldModel.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/media/IMediumModel.hpp

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -11,7 +10,7 @@
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/setup/SetupTrajectory.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
 
 namespace corsika {
 
@@ -21,12 +20,26 @@ namespace corsika {
 
     virtual MassDensityType getMassDensity(Point const&) const = 0;
 
-    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
-    // approach; for now, only lines are supported
-    virtual GrammageType getIntegratedGrammage(setup::Trajectory const&,
-                                               LengthType) const = 0;
-
-    virtual LengthType getArclengthFromGrammage(setup::Trajectory const&,
+    /**
+     * Integrate the matter density along trajectory.
+     *
+     * @return GrammageType as integrated matter density along the BaseTrajectory
+     *
+     * @todo think about the mixin inheritance of the trajectory vs the BaseTrajectory
+     *       approach; for now, only lines are supported (?).
+     */
+    virtual GrammageType getIntegratedGrammage(BaseTrajectory const&) const = 0;
+
+    /**
+     * Calculates the length along the trajectory.
+     *
+     * The length along the trajectory is determined at which the integrated matter
+     * density is reached. If the specified matter density cannot be reached (is too
+     * large) the result becomes meaningless and could be "infinity" (discuss this).
+     *
+     * @return LengthType the length corresponding to grammage.
+     */
+    virtual LengthType getArclengthFromGrammage(BaseTrajectory const&,
                                                 GrammageType) const = 0;
 
     virtual NuclearComposition const& getNuclearComposition() const = 0;

corsika/media/IMediumPropertyModel.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
@@ -16,10 +15,9 @@
 namespace corsika {
 
   /**
-   * An interface for type of media, needed e.g. to determine energy losses
+   * An interface for type of media, needed e.g. to determine energy losses.
    *
    * This is the base interface for media types.
-   *
    */
   template <typename TModel>
   class IMediumPropertyModel : public TModel {
@@ -31,7 +29,7 @@ namespace corsika {
      * @param  point    The location to evaluate at.
      * @returns    The media type
      */
-    virtual Medium getMedium(Point const&) const = 0;
+    virtual Medium getMedium() const = 0;
 
     /**
      * A virtual default destructor.