corsika/detail/framework/stack/Stack.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/stack/StackIteratorInterface.hpp>
+
+#include <stdexcept>
+#include <string>
+#include <vector>
+#include <utility>
+#include <type_traits>
+
+namespace corsika {
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... TArgs>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::clear(
+      TArgs... args) {
+    data_.clear(args...);
+    deleted_ = std::vector<bool>(data_.getSize(), false);
+    nDeleted_ = 0;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::begin() {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[i]) break;
+    }
+    return stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::end() {
+    return stack_iterator_type(*this, getSize());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::last() {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[getSize() - 1 - i]) break;
+    }
+    return stack_iterator_type(*this, getSize() - 1 - i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::begin() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[i]) break;
+    }
+    return const_stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::end() const {
+    return const_stack_iterator_type(*this, getSize());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface,
+                 MSecondaryProducer>::const_stack_iterator_type
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::last() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[getSize() - 1 - i]) break;
+    }
+    return const_stack_iterator_type(*this, getSize() - 1 - i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::cbegin() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[i]) break;
+    }
+    return const_stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::cend() const {
+    return const_stack_iterator_type(*this, getSize());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::clast() const {
+    unsigned int i = 0;
+    for (; i < getSize(); ++i) {
+      if (!deleted_[getSize() - 1 - i]) break;
+    }
+
+    return const_stack_iterator_type(*this, getSize() - 1 - i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::at(unsigned int i) {
+    return stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<
+          StackData, MParticleInterface, MSecondaryProducer>::at(unsigned int i) const {
+    return const_stack_iterator_type(*this, i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::first() {
+    return stack_iterator_type{*this, 0};
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      const_stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                             MSecondaryProducer>::cfirst() const {
+    return const_stack_iterator_type{*this, 0};
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<StackData, MParticleInterface,
+                                       MSecondaryProducer>::getNextParticle() {
+    while (purgeLastIfDeleted()) {}
+    return last();
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... TArgs>
+  typename Stack<StackData, MParticleInterface, MSecondaryProducer>::
+      stack_iterator_type inline Stack<
+          StackData, MParticleInterface,
+          MSecondaryProducer>::addParticle(const TArgs... v) {
+    data_.incrementSize();
+    deleted_.push_back(false);
+    return stack_iterator_type(*this, getSize() - 1, v...);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::swap(
+      stack_iterator_type a, stack_iterator_type b) {
+    swap(a.getIndex(), b.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::copy(
+      stack_iterator_type a, stack_iterator_type b) {
+    copy(a.getIndex(), b.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::copy(
+      const_stack_iterator_type a, stack_iterator_type b) {
+    data_.copy(a.getIndex(), b.getIndex());
+    if (deleted_[b.getIndex()] && !deleted_[a.getIndex()]) nDeleted_--;
+    if (!deleted_[b.getIndex()] && deleted_[a.getIndex()]) nDeleted_++;
+    deleted_[b.getIndex()] = deleted_[a.getIndex()];
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::erase(
+      stack_iterator_type p) {
+    if (this->isEmpty()) {
+      throw std::runtime_error("Stack, cannot delete entry since size is zero");
+    }
+    if (deleted_[p.getIndex()]) {
+      throw std::runtime_error("Stack, cannot delete entry since already deleted");
+    }
+    this->erase(p.getIndex());
+  }
+
+  /*
+   * delete this particle
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::erase(
+      particle_interface_type p) {
+    this->erase(p.getIterator());
+  }
+
+  /*
+   * check if there are no further non-deleted particles on stack
+   */
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isEmpty() {
+    return getEntries() == 0;
+  }
+
+  /*
+   * check if this particle was already deleted
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      const stack_iterator_type& p) const {
+    return isErased(p.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      const const_stack_iterator_type& p) const {
+    return isErased(p.getIndex());
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      const particle_interface_type& p) const {
+    return isErased(p.getIterator());
+  }
+
+  /*
+   * Function to ultimatively remove the last entry from the stack,
+   * if it was marked as deleted before. If this is not the case,
+   * the function will just return false and do nothing.
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::purgeLastIfDeleted() {
+    if (!deleted_.back())
+      return false; // the last particle is not marked for deletion. Do nothing.
+
+    CORSIKA_LOG_TRACE("stack: purgeLastIfDeleted: yes");
+    data_.decrementSize();
+    nDeleted_--;
+    deleted_.pop_back();
+    return true;
+  }
+
+  /*
+   * Function to ultimatively remove all entries from the stack
+   * marked as deleted.
+   *
+   * Careful: this will re-order the entries on the stack, since
+   * "gaps" in the stack are filled with entries from the back
+   * (copied).
+   */
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::purge() {
+    unsigned int iStackFront = 0;
+    unsigned int iStackBack = getSize() - 1;
+
+    for (unsigned int iDeleted = 0; iDeleted < getErased(); ++iDeleted) {
+      // search first delete entry on stack
+      while (!deleted_[iStackFront]) { iStackFront++; }
+      // search for last non-deleted particle on stack
+      while (deleted_[iStackBack]) { iStackBack--; }
+      // copy entry from iStackBack to iStackFront
+      data_.copy(iStackBack, iStackFront);
+      data_.decrementSize();
+    }
+    deleted_.clear();
+    nDeleted_ = 0;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline unsigned int Stack<StackData, MParticleInterface, MSecondaryProducer>::getSize()
+      const {
+    return data_.getSize();
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline std::string Stack<StackData, MParticleInterface, MSecondaryProducer>::asString()
+      const {
+    std::string str(fmt::format("size {}, entries {}, deleted {} \n", getSize(),
+                                getEntries(), getErased()));
+    // we make our own begin/end since we want ALL entries
+    std::string new_line = "     ";
+    for (unsigned int iPart = 0; iPart != getSize(); ++iPart) {
+      const_stack_iterator_type itPart(*this, iPart);
+      str += fmt::format("{}{}{}", new_line, itPart.asString(),
+                         (deleted_[itPart.getIndex()] ? " [deleted]" : ""));
+      new_line = "\n     ";
+    }
+    return str;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  template <typename... TArgs>
+  inline typename Stack<StackData, MParticleInterface,
+                        MSecondaryProducer>::stack_iterator_type
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::addSecondary(
+      stack_iterator_type& parent, const TArgs... v) {
+    data_.incrementSize();
+    deleted_.push_back(false);
+    return stack_iterator_type(*this, getSize() - 1, parent, v...);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::swap(
+      unsigned int const a, unsigned int const b) {
+    data_.swap(a, b);
+    std::vector<bool>::swap(deleted_[a], deleted_[b]);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::copy(
+      unsigned int const a, unsigned int const b) {
+    data_.copy(a, b);
+    if (deleted_[b] && !deleted_[a]) nDeleted_--;
+    if (!deleted_[b] && deleted_[a]) nDeleted_++;
+    deleted_[b] = deleted_[a];
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline bool Stack<StackData, MParticleInterface, MSecondaryProducer>::isErased(
+      unsigned int const i) const {
+    if (i >= deleted_.size()) return false;
+    return deleted_.at(i);
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::erase(
+      unsigned int const i) {
+    deleted_[i] = true;
+    nDeleted_++;
+  }
+
+  /*
+   * will remove from storage the element i. This is a helper
+   * function for SecondaryView.
+   */
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::purge(
+      unsigned int i) {
+    unsigned int iStackBack = getSize() - 1;
+    // search for last non-deleted particle on stack
+    while (deleted_[iStackBack]) { iStackBack--; }
+    // copy entry from iStackBack to iStackFront
+    data_.copy(iStackBack, i);
+    if (deleted_[i]) nDeleted_--;
+    deleted_[i] = deleted_[iStackBack];
+    data_.decrementSize();
+    deleted_.pop_back();
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline unsigned int
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::getIndexFromIterator(
+      const unsigned int vI) const {
+    return vI;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline typename Stack<StackData, MParticleInterface, MSecondaryProducer>::value_type&
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::getStackData() {
+    return data_;
+  }
+
+  template <typename StackData, template <typename> typename MParticleInterface,
+            template <typename T1, template <class> class T2> class MSecondaryProducer>
+  inline const typename Stack<StackData, MParticleInterface,
+                              MSecondaryProducer>::value_type&
+  Stack<StackData, MParticleInterface, MSecondaryProducer>::getStackData() const {
+    return data_;
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/COMBoost.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <cmath>
+
+#include <Eigen/Dense>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  inline COMBoost::COMBoost(FourMomentum const& P4projectile,
+                            HEPMassType const massTarget)
+      : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
+      , rotatedCS_{make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents())} {
+    auto const& pProjectile = P4projectile.getSpaceLikeComponents();
+    auto const pProjNormSquared = pProjectile.getSquaredNorm();
+    auto const pProjNorm = sqrt(pProjNormSquared);
+
+    auto const eProjectile = P4projectile.getTimeLikeComponent();
+    auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
+    auto const s =
+        massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
+
+    auto const sqrtS = sqrt(s);
+    auto const sinhEta = -pProjNorm / sqrtS;
+    auto const coshEta = sqrt(1 + pProjNormSquared / s);
+
+    setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
+
+  inline COMBoost::COMBoost(FourMomentum const& P4projectile,
+                            FourMomentum const& P4target)
+      : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()} {
+
+    // this is the center-of-momentum CM frame
+    auto const pCM =
+        P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents();
+    auto const pCM2 = pCM.getSquaredNorm();
+    auto const pCMnorm = sqrt(pCM2);
+    if (pCMnorm == 0_eV) {
+      // CM is at reset
+      rotatedCS_ = originalCS_;
+    } else {
+      rotatedCS_ = make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents() +
+                                                     P4target.getSpaceLikeComponents());
+    }
+
+    auto const s = (P4projectile + P4target).getNormSqr();
+    auto const sqrtS = sqrt(s);
+    auto const sinhEta = -pCMnorm / sqrtS;
+    auto const coshEta = sqrt(1 + pCM2 / s);
+
+    setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
+
+  inline COMBoost::COMBoost(MomentumVector const& momentum, HEPEnergyType const mass)
+      : originalCS_{momentum.getCoordinateSystem()}
+      , rotatedCS_{make_rotationToZ(momentum.getCoordinateSystem(), momentum)} {
+    auto const squaredNorm = momentum.getSquaredNorm();
+    auto const norm = sqrt(squaredNorm);
+    auto const sinhEta = -norm / mass;
+    auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
+    setBoost(coshEta, sinhEta);
+    CORSIKA_LOG_TRACE("COMBoost (1-beta)={}, gamma={}, det={}", 1 - sinhEta / coshEta,
+                      coshEta, boost_.determinant() - 1);
+  }
+
+  template <typename FourVector>
+  inline FourVector COMBoost::toCoM(FourVector const& p4) const {
+    auto const pComponents = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+    Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
+    Eigen::Vector2d lab;
+
+    lab << (p4.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();
+
+    CORSIKA_LOG_TRACE("E0={}, p={}, E0'={}, p'={}", p4.getTimeLikeComponent() / 1_GeV,
+                      eVecRotated(2) * (1 / 1_GeV).magnitude(), E_CoM / 1_GeV, boostedZ);
+
+    return FourVector(E_CoM, MomentumVector(rotatedCS_, eVecRotated));
+  }
+
+  template <typename FourVector>
+  inline FourVector COMBoost::fromCoM(FourVector const& p4) const {
+    auto pCM = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+    auto const Ecm = p4.getTimeLikeComponent();
+
+    Eigen::Vector2d com;
+    com << (Ecm * (1 / 1_GeV)), (pCM.getEigenVector()(2) * (1 / 1_GeV).magnitude());
+
+    CORSIKA_LOG_TRACE("Ecm={} GeV, pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
+                      Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV,
+                      p4.getNorm() / 1_GeV);
+
+    auto const boostedZ = inverseBoost_ * com;
+    auto const E_lab = boostedZ(0) * 1_GeV;
+
+    pCM.getEigenVector()[2] = boostedZ(1) * (1_GeV).magnitude();
+
+    Vector<typename decltype(pCM)::dimension_type> pLab{rotatedCS_, pCM};
+    pLab.rebase(originalCS_);
+
+    CORSIKA_LOG_TRACE("COMBoost::fromCoM --> Elab={} GeV",
+                      " plab={} GeV (norm={} GeV) "
+                      " GeV), invariant mass = {}",
+                      E_lab / 1_GeV, FourVector{E_lab, pLab}.getNorm() / 1_GeV,
+                      pLab.getComponents(), pLab.getNorm() / 1_GeV);
+
+    return FourVector{E_lab, pLab};
+  }
+
+  inline void COMBoost::setBoost(double const coshEta, double const sinhEta) {
+    boost_ << coshEta, sinhEta, sinhEta, coshEta;
+    inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
+  }
+
+  inline CoordinateSystemPtr const& COMBoost::getRotatedCS() const { return rotatedCS_; }
+
+  inline CoordinateSystemPtr const& COMBoost::getOriginalCS() const {
+    return originalCS_;
+  }
+
+} // namespace corsika

corsika/detail/framework/utility/CorsikaData.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <corsika/corsika.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <boost/filesystem/path.hpp>
+
+#include <cstdlib>
+#include <stdexcept>
+#include <string>
+
+namespace corsika {
+
+  inline boost::filesystem::path corsika_data(boost::filesystem::path const& key) {
+    boost::filesystem::path fname =
+        boost::filesystem::path(corsika::CORSIKA_DATA_DIR) / key;
+    // LCOV_EXCL_START, this cannot be easily tested system-independently
+    if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
+      fname = boost::filesystem::path(p) / key;
+    }
+    // LCOV_EXCL_STOP
+    CORSIKA_LOG_INFO("opening data file={}", fname);
+    return fname;
+  }
+} // namespace corsika

corsika/detail/framework/utility/CorsikaFenvFallback.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * 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
+
+inline int feenableexcept(int /*excepts*/) noexcept { return -1; }
+
+inline int fedisableexcept(int /*excepts*/) noexcept { return -1; }
+}

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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+/**
+ * Import public domain code
+ *
+ * Provide portable or fallback versions of feenableexcept() and fedisableexcept()
+ * Exist by default in glibc since version 2.2, but not in the standard
+ * fenv.h / cfenv headers for C 99 or C++ 11
+ *
+ * \author Lukas Nellen
+ * \date 14 Jan 2019
+ *
+ */
+
+#include <cfenv>
+
+// Implementation for OS X on intel X64_86
+// code from
+// https://stackoverflow.com/questions/37819235/how-do-you-enable-floating-point-exceptions-for-clang-in-os-x
+// based on
+// http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c
+
+extern "C" {
+
+inline int feenableexcept(int excepts) noexcept {
+  static fenv_t fenv;
+  int new_excepts = excepts & FE_ALL_EXCEPT;
+  // previous masks
+  int old_excepts;
+
+  if (fegetenv(&fenv)) { return -1; }
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
+
+  // unmask
+  fenv.__control &= ~new_excepts;
+  fenv.__mxcsr &= ~(new_excepts << 7);
+
+  return fesetenv(&fenv) ? -1 : old_excepts;
+}
+
+inline int fedisableexcept(int excepts) noexcept {
+  static fenv_t fenv;
+  int new_excepts = excepts & FE_ALL_EXCEPT;
+  // all previous masks
+  int old_excepts;
+
+  if (fegetenv(&fenv)) { return -1; }
+  old_excepts = fenv.__control & FE_ALL_EXCEPT;
+
+  // mask
+  fenv.__control |= new_excepts;
+  fenv.__mxcsr |= new_excepts << 7;
+
+  return fesetenv(&fenv) ? -1 : old_excepts;
+}
+}

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 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 {
+
+    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];
+      }
+
+      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
+
+  inline std::vector<double> solve_quartic_depressed_real(long double p, long double q,
+                                                          long double r,
+                                                          double const epsilon) {
+
+    CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f},  epsilon={}", p, q, r,
+                      epsilon);
+
+    long double const p2 = static_pow<2>(p);
+    long double const q2 = static_pow<2>(q);
+
+    std::vector<double> const resolve_cubic =
+        solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
+
+    CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(),
+                      fmt::join(resolve_cubic, ", "));
+
+    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
+      }
+    }
+    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) {
+
+    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);
+    }
+
+    if ((std::abs(a - 1) < epsilon) &&
+        (std::abs(b) < epsilon)) { // this is a depressed quartic
+      return solve_quartic_depressed_real(c, d, e, epsilon);
+    }
+
+    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

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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#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>
+#include <numeric>
+#include <utility>
+#include <vector>
+#include <string>
+
+namespace corsika {
+
+  template <class Axes, class Storage>
+  inline void save_hist(boost::histogram::histogram<Axes, Storage> const& h,
+                        std::string const& filename, bool overwrite) {
+    if (boost::filesystem::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());
+      }
+    }
+
+    unsigned const rank = h.rank();
+
+    std::vector<size_t> axes_dims;
+    axes_dims.reserve(rank);
+
+    auto overflow = std::make_unique<bool[]>(rank);
+    auto underflow = std::make_unique<bool[]>(rank);
+
+    std::vector<char> axis_types;
+    axis_types.reserve(rank);
+
+    for (unsigned int i = 0; i < rank; ++i) {
+      auto const& ax = h.axis(i);
+      unsigned const has_underflow = (ax.options() & 0x01) ? 1 : 0;
+      unsigned const has_overflow = (ax.options() & 0x02) ? 1 : 0;
+
+      underflow[i] = has_underflow;
+      overflow[i] = has_overflow;
+
+      axes_dims.emplace_back(ax.size() + has_underflow + has_overflow);
+
+      if (ax.continuous()) {
+        axis_types.push_back('c');
+        std::vector<double> ax_edges;
+        ax_edges.reserve(ax.size() + 1);
+
+        for (int j = 0; j < ax.size(); ++j) { ax_edges.push_back(ax.bin(j).lower()); }
+
+        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()}, "a");
+      } else {
+        axis_types.push_back('d');
+        std::vector<int64_t> bins; // we assume that discrete axes have integer bins
+        bins.reserve(ax.size());
+
+        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()}, "a");
+      }
+    }
+
+    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<>());
+    auto temp = std::make_unique<float[]>(prod_axis_size);
+
+    assert(prod_axis_size == h.size());
+
+    // reduce multi-dim. to 1-dim, row-major (i.e., last axis index is contiguous in
+    // memory) take special care of underflow bins which have -1 as index and thus need
+    // to be shifted by +1
+
+    for (auto&& x : indexed(h, boost::histogram::coverage::all)) {
+      int p = 0;
+
+      for (unsigned axis_index = 0; axis_index < rank; ++axis_index) {
+        int const offset_underflow = (h.axis(axis_index).options() & 0x01) ? 1 : 0;
+        auto k = x.index(axis_index) + offset_underflow;
+        p = k + p * axes_dims.at(axis_index);
+      }
+
+      temp[p] = *x;
+    }
+
+    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

corsika/detail/media/BaseExponential.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>
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline BaseExponential<TDerived>::BaseExponential(Point const& point,
+                                                    LengthType const referenceHeight,
+                                                    MassDensityType const rho0,
+                                                    LengthType const lambda)
+      : rho0_(rho0)
+      , lambda_(lambda)
+      , invLambda_(1 / lambda)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {}
+
+  template <typename TDerived>
+  inline auto const& BaseExponential<TDerived>::getImplementation() const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+      LengthType const height) const {
+    return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+  }
+
+  template <typename TDerived>
+  inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj, 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 length * rhoStart;
+    } else {
+      return rhoStart * (lambda_ / uDotA) * expm1(uDotA * length * invLambda_);
+    }
+  }
+
+  template <typename TDerived>
+  inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage,
+      DirectionVector const& axis) const {
+    // 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;
+      if (logArg > -1) {
+        return lambda_ / uDotA * log1p(logArg);
+      } else {
+        return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
+               meter;
+      }
+    }
+  }
+
+} // 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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/Environment.hpp>
+
+namespace corsika {
+
+  template <typename IEnvironmentModel>
+  inline Environment<IEnvironmentModel>::Environment()
+      : coordinateSystem_{get_root_CoordinateSystem()}
+      , universe_(std::make_unique<BaseNodeType>(
+            std::make_unique<Universe>(coordinateSystem_))) {}
+
+  template <typename IEnvironmentModel>
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr&
+  Environment<IEnvironmentModel>::getUniverse() {
+    return universe_;
+  }
+
+  template <typename IEnvironmentModel>
+  inline typename Environment<IEnvironmentModel>::BaseNodeType::VTNUPtr const&
+  Environment<IEnvironmentModel>::getUniverse() const {
+    return universe_;
+  }
+
+  template <typename IEnvironmentModel>
+  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> > inline Environment<
+      IEnvironmentModel>::createNode(TVolumeArgs&&... args) {
+    static_assert(std::is_base_of_v<IVolume, TVolumeType>,
+                  "unusable type provided, needs to be derived from "
+                  "\"Volume\"");
+
+    return std::make_unique<BaseNodeType>(
+        std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
+  }
+
+  template <typename IEnvironmentModel>
+  std::set<Code> const get_all_elements_in_universe(
+      Environment<IEnvironmentModel> const& env) {
+
+    auto const& universe = *(env.getUniverse());
+    auto const allElementsInUniverse = std::invoke([&]() {
+      std::set<Code> allElementsInUniverse;
+      auto collectElements = [&](auto& vtn) {
+        if (vtn.hasModelProperties()) {
+          auto const& comp =
+              vtn.getModelProperties().getNuclearComposition().getComponents();
+          for (auto const c : comp) allElementsInUniverse.insert(c);
+        }
+      };
+      universe.walk(collectElements);
+      return allElementsInUniverse;
+    });
+    return allElementsInUniverse;
+  }
+
+} // namespace corsika

corsika/detail/media/ExponentialRefractiveIndex.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 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 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/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/BaseExponential.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  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>
+  inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
+    return BaseExponential<FlatExponential<T>>::getMassDensity(
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
+  }
+
+  template <typename T>
+  inline NuclearComposition const& FlatExponential<T>::getNuclearComposition() const {
+    return nuclComp_;
+  }
+
+  template <typename T>
+  inline GrammageType FlatExponential<T>::getIntegratedGrammage(
+      BaseTrajectory const& line) const {
+    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
+  }
+
+  template <typename T>
+  inline LengthType FlatExponential<T>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType const grammage) const {
+    return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
+                                                                         axis_);
+  }
+
+} // namespace corsika
+
+#include <corsika/detail/media/BaseExponential.inl>

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 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/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  inline HomogeneousMedium<T>::HomogeneousMedium(MassDensityType density,
+                                                 NuclearComposition const& nuclComp)
+      : density_(density)
+      , nuclComp_(nuclComp) {}
+
+  template <typename T>
+  inline MassDensityType HomogeneousMedium<T>::getMassDensity(Point const&) const {
+    return density_;
+  }
+
+  template <typename T>
+  inline NuclearComposition const& HomogeneousMedium<T>::getNuclearComposition() const {
+    return nuclComp_;
+  }
+
+  template <typename T>
+  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+      BaseTrajectory const& track) const {
+    return track.getLength() * density_;
+  }
+
+  template <typename T>
+  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 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/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+
+namespace corsika {
+
+  template <typename T, typename TDensityFunction>
+  template <typename... TArgs>
+  inline InhomogeneousMedium<T, TDensityFunction>::InhomogeneousMedium(
+      NuclearComposition const& nuclComp, TArgs&&... rhoTArgs)
+      : nuclComp_(nuclComp)
+      , densityFunction_(rhoTArgs...) {}
+
+  template <typename T, typename TDensityFunction>
+  inline MassDensityType InhomogeneousMedium<T, TDensityFunction>::getMassDensity(
+      Point const& point) const {
+    return densityFunction_.evaluateAt(point);
+  }
+
+  template <typename T, typename TDensityFunction>
+  inline NuclearComposition const&
+  InhomogeneousMedium<T, TDensityFunction>::getNuclearComposition() const {
+    return nuclComp_;
+  }
+
+  template <typename T, typename TDensityFunction>
+  inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
+      BaseTrajectory const& line) const {
+    return densityFunction_.getIntegrateGrammage(line);
+  }
+
+  template <typename T, typename TDensityFunction>
+  inline LengthType InhomogeneousMedium<T, TDensityFunction>::getArclengthFromGrammage(
+      BaseTrajectory const& line, GrammageType grammage) const {
+    return densityFunction_.getArclengthFromGrammage(line, grammage);
+  }
+
+} // namespace corsika