corsika/detail/media/GeomagneticModel.inl

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#include <corsika/framework/core/Logging.hpp>
+#include <boost/math/tr1.hpp>
+#include <boost/filesystem.hpp>
+#include <boost/filesystem/fstream.hpp>
+
+#include <stdexcept>
+#include <string>
+#include <cmath>
+
+namespace corsika {
+
+  inline GeomagneticModel::GeomagneticModel(Point const& center,
+                                            boost::filesystem::path const path)
+      : center_(center) {
+
+    // Read in coefficients
+    boost::filesystem::ifstream file(path, std::ios::in);
+
+    // Exit if file opening failed
+    if (!file.is_open()) {
+      CORSIKA_LOG_ERROR("Failed opening data file {}", path);
+      throw std::runtime_error("Cannot load GeomagneticModel data.");
+    }
+
+    // GeomagneticModel supports two types of input data: WMM.COF and IGRF.COF
+    // They have only slightly different format and content and can be easily
+    // differentiated here.
+
+    std::string line;
+    while (getline(file >> std::ws, line)) {
+
+      double epoch;
+      std::string model_name;
+      std::string release_date; // just for WMM
+      int nMax = 12; // the spherical max n (l) shell (for IGRF), for WMM this is 12
+                     // Note that n=l=0 is the monopole and is not included in the model.
+      int dummyInt;
+      double dummyDouble;
+
+      std::istringstream firstLine(line);
+
+      // check comments and ignore:
+      if (firstLine.peek() == '#' ||                     // normal comment
+          line.size() == 0 ||                            // empty line
+          line.find("9999999999999999999999999") == 0) { // crazy WMM comment
+        continue;
+      }
+
+      // check IGRF format:
+      if (firstLine >> model_name >> epoch >> nMax >> dummyInt >> dummyInt >>
+          dummyDouble >> dummyDouble >> dummyDouble >> dummyDouble >> model_name >>
+          dummyInt) {
+        static bool info = false;
+        if (!info) {
+          CORSIKA_LOG_INFO("Reading IGRF input data format.");
+          info = true;
+        }
+      } else {
+        // check WMM format:
+        firstLine.clear();
+        firstLine.seekg(0, std::ios::beg);
+        if (firstLine >> epoch >> model_name >> release_date) {
+          CORSIKA_LOG_INFO("Reading WMM input data format.");
+        } else {
+          CORSIKA_LOG_ERROR("line: {}", line);
+          throw std::runtime_error("Incompatible input data for GeomagneticModel");
+        }
+      }
+
+      int nPar = 0;
+      for (int i = 0; i < nMax; ++i) { nPar += i + 2; }
+      int iEpoch = int(epoch);
+
+      if (parameters_.count(iEpoch) != 0) {
+        throw std::runtime_error("GeomagneticModel input file has duplicate Epoch. Fix.");
+      }
+      parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
+
+      for (int i = 0; i < nPar; i++) {
+        file >> parameters_[iEpoch][i].n >> parameters_[iEpoch][i].m >>
+            parameters_[iEpoch][i].g >> parameters_[iEpoch][i].h >>
+            parameters_[iEpoch][i].dg >> parameters_[iEpoch][i].dh;
+        file.ignore(9999999, '\n');
+      }
+    }
+    file.close();
+
+    if (parameters_.size() == 0) {
+      CORSIKA_LOG_ERROR("No input data read!");
+      throw std::runtime_error("No input data read");
+    }
+  }
+
+  inline MagneticFieldVector GeomagneticModel::getField(double const year,
+                                                        LengthType const altitude,
+                                                        double const latitude,
+                                                        double const longitude) {
+
+    int iYear = int(year);
+    auto iEpoch = parameters_.rbegin();
+    for (; iEpoch != parameters_.rend(); ++iEpoch) {
+      if (iEpoch->first <= iYear) { break; }
+    }
+    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch->first, year);
+    if (iEpoch == parameters_.rend()) {
+      CORSIKA_LOG_WARN("Year {} is before first EPOCH. Results unclear.", year);
+      iEpoch--; // move one epoch back
+    }
+    if (altitude < -1_km || altitude > 600_km) {
+      CORSIKA_LOG_WARN("Altitude should be between -1_km and 600_km.");
+    }
+    if (latitude <= -90 || latitude >= 90) {
+      CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
+      throw std::runtime_error("Latitude has to be between -90 and 90 degree.");
+    } else if (latitude < -89.992 || latitude > 89.992) {
+      CORSIKA_LOG_WARN("Latitude is close to the poles.");
+    }
+    if (longitude < -180 || longitude > 180) {
+      CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
+    }
+
+    double const epoch = double(iEpoch->first);
+    auto iNextEpoch = iEpoch; // next epoch for interpolation
+    --iNextEpoch;
+    bool const lastEpoch = (iEpoch == parameters_.rbegin());
+    auto const delta_t = year - epoch;
+    CORSIKA_LOG_DEBUG(
+        "identified: t_epoch={}, delta_t={}, lastEpoch={} (false->interpolate)", epoch,
+        delta_t, lastEpoch);
+
+    double const lat_geo = latitude * constants::pi / 180;
+    double const lon = longitude * constants::pi / 180;
+
+    // Transform into spherical coordinates
+    double constexpr f = 1 / 298.257223563;
+    double constexpr e_squared = f * (2 - f);
+    LengthType R_c =
+        constants::EarthRadius::Equatorial / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
+    LengthType p = (R_c + altitude) * cos(lat_geo);
+    LengthType z = sin(lat_geo) * (altitude + R_c * (1 - e_squared));
+    LengthType r = sqrt(p * p + z * z);
+    double lat_sph = asin(z / r);
+
+    double legendre, next_legendre, derivate_legendre;
+    double magneticfield[3] = {0, 0, 0};
+
+    for (size_t j = 0; j < iEpoch->second.size(); j++) {
+
+      ParameterLine p = iEpoch->second[j];
+
+      // Time interpolation
+      if (iEpoch == parameters_.rbegin()) {
+        // this is the latest epoch in time, or time-dependence (dg/dh) was specified
+        // we use the extrapolation factors dg/dh:
+        p.g = p.g + delta_t * p.dg;
+        p.h = p.h + delta_t * p.dh;
+      } else {
+        // we linearly interpolate between two epochs
+        ParameterLine const next_p = iNextEpoch->second[j];
+        double const length = iNextEpoch->first - epoch;
+        double p_g = p.g + (next_p.g - p.g) * delta_t / length;
+        double p_h = p.h + (next_p.h - p.h) * delta_t / length;
+        CORSIKA_LOG_TRACE(
+            "interpolation: delta-g={}, delta-h={}, delta-t={}, length={} g1={} g2={} "
+            "g={} h={} ",
+            next_p.g - p.g, next_p.h - p.h, year - epoch, length, next_p.g, p.g, p_g,
+            p_h);
+        p.g = p_g;
+        p.h = p_h;
+      }
+
+      legendre = boost::math::tr1::assoc_legendre(p.n, p.m, sin(lat_sph));
+      next_legendre = boost::math::tr1::assoc_legendre(p.n + 1, p.m, sin(lat_sph));
+
+      // Schmidt semi-normalization
+      if (p.m > 0) {
+        // Note: n! = tgamma(n+1)
+        legendre *= sqrt(2 * std::tgamma(p.n - p.m + 1) / std::tgamma(p.n + p.m + 1));
+        next_legendre *=
+            sqrt(2 * std::tgamma(p.n + 1 - p.m + 1) / std::tgamma(p.n + 1 + p.m + 1));
+      }
+      derivate_legendre =
+          (p.n + 1) * tan(lat_sph) * legendre -
+          sqrt(pow(p.n + 1, 2) - pow(p.m, 2)) / cos(lat_sph) * next_legendre;
+
+      magneticfield[0] +=
+          pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+          (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * derivate_legendre;
+      magneticfield[1] +=
+          pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) * p.m *
+          (p.g * sin(p.m * lon) - p.h * cos(p.m * lon)) * legendre;
+      magneticfield[2] +=
+          (p.n + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+          (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * legendre;
+    }
+    magneticfield[0] *= -1;
+    magneticfield[1] /= cos(lat_sph);
+    magneticfield[2] *= -1;
+
+    // Transform back into geodetic coordinates
+    double magneticfield_geo[3];
+    magneticfield_geo[0] = magneticfield[0] * cos(lat_sph - lat_geo) -
+                           magneticfield[2] * sin(lat_sph - lat_geo);
+    magneticfield_geo[1] = magneticfield[1];
+    magneticfield_geo[2] = magneticfield[0] * sin(lat_sph - lat_geo) +
+                           magneticfield[2] * cos(lat_sph - lat_geo);
+
+    return MagneticFieldVector{center_.getCoordinateSystem(), magneticfield_geo[0] * 1_nT,
+                               magneticfield_geo[1] * -1_nT,
+                               magneticfield_geo[2] * -1_nT};
+  }
+
+} // namespace corsika

corsika/detail/media/GladstoneDaleRefractiveIndex.inl

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  inline GladstoneDaleRefractiveIndex<T>::GladstoneDaleRefractiveIndex(
+      double const referenceRefractiveIndex, Point const point, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , referenceRefractivity_(referenceRefractiveIndex - 1)
+      , referenceInvDensity_(1 / this->getMassDensity(point)) {}
+
+  template <typename T>
+  inline double GladstoneDaleRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return referenceRefractivity_ * (this->getMassDensity(point) * referenceInvDensity_) +
+           1.;
+  }
+
+} // namespace corsika

corsika/detail/media/HomogeneousMedium.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -32,9 +31,9 @@ namespace corsika {
   }
 
   template <typename T>
-  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(BaseTrajectory const&,
-                                                                  LengthType to) const {
-    return to * density_;
+  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+      BaseTrajectory const& track) const {
+    return track.getLength() * density_;
   }
 
   template <typename T>

corsika/detail/media/InhomogeneousMedium.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -36,8 +35,8 @@ namespace corsika {
 
   template <typename T, typename TDensityFunction>
   inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
-      BaseTrajectory const& line, LengthType to) const {
-    return densityFunction_.getIntegrateGrammage(line, to);
+      BaseTrajectory const& line) const {
+    return densityFunction_.getIntegrateGrammage(line);
   }
 
   template <typename T, typename TDensityFunction>

corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -13,13 +12,14 @@
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
 #include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
 
 namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  inline void LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
-                                                TModelArgs...>::checkRadius(LengthType r)
+  inline void LayeredSphericalAtmosphereBuilder<
+      TMediumInterface, TMediumModelExtra, TModelArgs...>::checkRadius(LengthType const r)
       const {
     if (r <= previousRadius_) {
       throw std::runtime_error("radius must be greater than previous");
@@ -36,26 +36,28 @@ namespace corsika {
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
-  inline void LayeredSphericalAtmosphereBuilder<
-      TMediumInterface, TMediumModelExtra,
-      TModelArgs...>::addExponentialLayer(GrammageType b, LengthType c,
-                                          LengthType upperBoundary) {
-
-    auto const radius = earthRadius_ + upperBoundary;
+  inline typename LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
+                                                    TModelArgs...>::volume_tree_node*
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addExponentialLayer(GrammageType const b, LengthType const scaleHeight,
+                          LengthType const upperBoundary) {
+
+    // outer radius
+    auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
 
     auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
         std::make_unique<Sphere>(center_, radius));
 
-    auto const rho0 = b / c;
+    auto const rho0 = b / scaleHeight;
 
     if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
       // helper lambda in which the last 5 arguments to make_shared<...> are bound
       auto lastBound = [&](auto... argPack) {
         return std::make_shared<
             TMediumModelExtra<SlidingPlanarExponential<TMediumInterface>>>(
-            argPack..., center_, rho0, -c, *composition_, earthRadius_);
+            argPack..., center_, rho0, -scaleHeight, *composition_, planetRadius_);
       };
 
       // now unpack the additional arguments
@@ -63,26 +65,28 @@ namespace corsika {
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
-          center_, rho0, -c, *composition_, earthRadius_);
+          center_, rho0, -scaleHeight, *composition_, planetRadius_);
     }
 
     layers_.push(std::move(node));
+    return layers_.top().get();
   }
 
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
   inline void LayeredSphericalAtmosphereBuilder<
       TMediumInterface, TMediumModelExtra,
-      TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
-    auto const radius = earthRadius_ + upperBoundary;
+      TModelArgs...>::addLinearLayer(GrammageType const b, LengthType const scaleHeight,
+                                     LengthType const upperBoundary) {
+    // outer radius
+    auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
 
     auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
         std::make_unique<Sphere>(center_, radius));
 
-    units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
-    auto const rho0 = b / c;
+    auto const rho0 = b / scaleHeight;
 
     if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
       // helper lambda in which the last 2 arguments to make_shared<...> are bound
@@ -93,7 +97,6 @@ namespace corsika {
 
       // now unpack the additional arguments
       auto model = std::apply(lastBound, additionalModelArgs_);
-
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
@@ -103,6 +106,40 @@ namespace corsika {
     layers_.push(std::move(node));
   }
 
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+                      unsigned int const nBins, LengthType const deltaHeight,
+                      LengthType const upperBoundary) {
+
+    auto const radius = planetRadius_ + upperBoundary;
+    checkRadius(radius);
+    previousRadius_ = radius;
+
+    auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+        std::make_unique<Sphere>(center_, radius));
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 5 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<
+            TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
+            argPack..., center_, funcRho, nBins, deltaHeight, *composition_,
+            planetRadius_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
+          center_, funcRho, nBins, deltaHeight, *composition_, planetRadius_);
+    }
+    layers_.push(std::move(node));
+  }
+
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
   inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
@@ -132,10 +169,11 @@ namespace corsika {
   template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
   struct make_layered_spherical_atmosphere_builder {
     template <typename... TArgs>
-    static auto create(Point const& center, LengthType earthRadius, TArgs... args) {
+    static auto create(Point const& center, LengthType const planetRadius,
+                       TArgs... args) {
       return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
                                                TArgs...>{std::forward<TArgs>(args)...,
-                                                         center, earthRadius};
+                                                         center, planetRadius};
     }
   };
 

corsika/detail/media/LinearApproximationIntegrator.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -13,21 +12,26 @@
 namespace corsika {
 
   template <typename TDerived>
-  inline auto const& LinearApproximationIntegrator<TDerived>::getImplementation() const {
+  inline TDerived const& LinearApproximationIntegrator<TDerived>::getImplementation()
+      const {
     return *static_cast<TDerived const*>(this);
   }
 
   template <typename TDerived>
-  inline auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
-      BaseTrajectory const& line, LengthType length) const {
+  inline GrammageType LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
+      BaseTrajectory const& line) const {
+    LengthType const length = line.getLength();
     auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
     auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
                                                                 line.getDirection(0));
+    CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1,
+                     line.getPosition(0), line.getDirection(0),
+                     (c0 + 0.5 * c1 * length) * length);
     return (c0 + 0.5 * c1 * length) * length;
   }
 
   template <typename TDerived>
-  inline auto LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
+  inline LengthType LinearApproximationIntegrator<TDerived>::getArclengthFromGrammage(
       BaseTrajectory const& line, GrammageType grammage) const {
     auto const c0 = getImplementation().rho_(line.getPosition(0));
     auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
@@ -37,7 +41,7 @@ namespace corsika {
   }
 
   template <typename TDerived>
-  inline auto LinearApproximationIntegrator<TDerived>::getMaximumLength(
+  inline LengthType LinearApproximationIntegrator<TDerived>::getMaximumLength(
       BaseTrajectory const& line, [[maybe_unused]] double relError) const {
     [[maybe_unused]] auto const c1 = getImplementation().rho_.getSecondDerivative(
         line.getPosition(0), line.getDirection(0));

corsika/detail/media/MediumPropertyModel.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -19,7 +18,7 @@ namespace corsika {
       , medium_(medium) {}
 
   template <typename T>
-  inline Medium MediumPropertyModel<T>::getMedium(Point const&) const {
+  inline Medium MediumPropertyModel<T>::getMedium() const {
     return medium_;
   }
 

corsika/detail/media/NuclearComposition.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -11,7 +10,8 @@
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 
-#include <corsika/media/WeightProvider.hpp>
+#include <boost/iterator/zip_iterator.hpp>
+#include <boost/iterator/transform_iterator.hpp>
 
 #include <cassert>
 #include <functional>
@@ -23,31 +23,52 @@
 namespace corsika {
 
   inline NuclearComposition::NuclearComposition(std::vector<Code> const& pComponents,
-                                                std::vector<float> const& pFractions)
+                                                std::vector<double> const& pFractions)
       : numberFractions_(pFractions)
       , components_(pComponents)
-      , avgMassNumber_(std::inner_product(
-            pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
-            std::plus<double>(), [](auto const compID, auto const fraction) -> double {
-              if (is_nucleus(compID)) {
-                return get_nucleus_A(compID) * fraction;
-              } else {
-                return get_mass(compID) / convert_SI_to_HEP(constants::u) * fraction;
-              }
-            })) {
-    assert(pComponents.size() == pFractions.size());
-    auto const sumFractions =
-        std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-    if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
+      , avgMassNumber_(getWeightedSum([](Code const compID) -> double {
+        if (is_nucleus(compID)) {
+          return get_nucleus_A(compID);
+        } else {
+          return get_mass(compID) / convert_SI_to_HEP(constants::u);
+        }
+      })) {
+    if (pComponents.size() != pFractions.size()) {
+      throw std::runtime_error(
+          "Cannot construct NuclearComposition from vectors of different sizes.");
+    }
+    auto const sumFractions = std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.);
+
+    if (!(0.999 < sumFractions && sumFractions < 1.001)) {
       throw std::runtime_error("element fractions do not add up to 1");
     }
     this->updateHash();
   }
 
   template <typename TFunction>
-  inline auto NuclearComposition::getWeightedSum(TFunction const& func) const {
-    using ResultQuantity = decltype(func(*components_.cbegin()));
+  inline auto NuclearComposition::getWeighted(TFunction func) const {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
+    auto const product = [&](auto const compID, auto const fraction) {
+      return func(compID) * fraction;
+    };
+
+    if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity::zero());
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    } else {
+      std::vector<ResultQuantity> result(components_.size(), ResultQuantity(0));
+      std::transform(components_.cbegin(), components_.cend(), numberFractions_.cbegin(),
+                     result.begin(), product);
+      return result;
+    }
+  } // namespace corsika
+
+  template <typename TFunction>
+  inline auto NuclearComposition::getWeightedSum(TFunction func) const
+      -> decltype(func(std::declval<Code>())) {
+    using ResultQuantity = decltype(func(std::declval<Code>()));
 
     auto const prod = [&](auto const compID, auto const fraction) {
       return func(compID) * fraction;
@@ -68,7 +89,7 @@ namespace corsika {
 
   inline size_t NuclearComposition::getSize() const { return numberFractions_.size(); }
 
-  inline std::vector<float> const& NuclearComposition::getFractions() const {
+  inline std::vector<double> const& NuclearComposition::getFractions() const {
     return numberFractions_;
   }
 
@@ -82,16 +103,28 @@ namespace corsika {
 
   template <class TRNG>
   inline Code NuclearComposition::sampleTarget(std::vector<CrossSectionType> const& sigma,
-                                               TRNG& randomStream) const {
+                                               TRNG&& randomStream) const {
+    if (sigma.size() != numberFractions_.size()) {
+      throw std::runtime_error("incompatible vector sigma as input");
+    }
+
+    auto zip_beg = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cbegin(), sigma.cbegin()));
+    auto zip_end = boost::make_zip_iterator(
+        boost::make_tuple(numberFractions_.cend(), sigma.cend()));
+    using zip_iter_type = decltype(zip_beg);
+
+    auto const mult_func = [](zip_iter_type::value_type const& zipit) -> double {
+      return zipit.get<0>() * zipit.get<1>().magnitude();
+    };
+
+    using transform_iter_type =
+        boost::transform_iterator<decltype(mult_func), zip_iter_type, double, double>;
 
-    assert(sigma.size() == numberFractions_.size());
+    auto trans_beg = transform_iter_type{zip_beg, mult_func};
+    auto trans_end = transform_iter_type{zip_end, mult_func};
 
-    std::discrete_distribution channelDist(
-        WeightProviderIterator<decltype(numberFractions_.begin()),
-                               decltype(sigma.begin())>(numberFractions_.begin(),
-                                                        sigma.begin()),
-        WeightProviderIterator<decltype(numberFractions_.begin()), decltype(sigma.end())>(
-            numberFractions_.end(), sigma.end()));
+    std::discrete_distribution channelDist{trans_beg, trans_end};
 
     auto const iChannel = channelDist(randomStream);
     return components_[iChannel];
@@ -99,6 +132,7 @@ namespace corsika {
 
   // Note: when this class ever modifies its internal data, the hash
   // must be updated, too!
+  // the hash value is important to find tables, etc.
   inline size_t NuclearComposition::getHash() const { return hash_; }
 
   inline bool NuclearComposition::operator==(NuclearComposition const& v) const {
@@ -107,7 +141,8 @@ namespace corsika {
 
   inline void NuclearComposition::updateHash() {
     std::vector<std::size_t> hashes;
-    for (float ifrac : this->getFractions()) hashes.push_back(std::hash<float>{}(ifrac));
+    for (double ifrac : this->getFractions())
+      hashes.push_back(std::hash<double>{}(ifrac));
     for (Code icode : this->getComponents())
       hashes.push_back(std::hash<int>{}(static_cast<int>(icode)));
     std::size_t h = std::hash<double>{}(this->getAverageMassNumber());

corsika/detail/media/ShowerAxis.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/framework/core/PhysicalUnits.hpp>
@@ -26,9 +25,16 @@ namespace corsika {
       , X_(steps + 1) {
     auto const* const universe = env.getUniverse().get();
 
-    auto rho = [pStart, length, universe](double x) {
+    auto rho = [pStart, length, universe, doThrow](double x) {
       auto const p = pStart + length * x;
       auto const* node = universe->getContainingNode(p);
+      if (!node->hasModelProperties()) {
+        CORSIKA_LOG_CRITICAL(
+            "Unable to construct ShowerAxis. ShowerAxis includes volume "
+            "with no model properties at point {}.",
+            p);
+        if (doThrow) throw std::runtime_error("Unable to construct ShowerAxis.");
+      }
       return node->getModelProperties().getMassDensity(p).magnitude();
     };
 
@@ -65,7 +71,7 @@ namespace corsika {
     unsigned int const upper = lower + 1;
 
     if (fractionalBin < 0) {
-      CORSIKA_LOG_ERROR("cannot extrapolate to points behind point of injection l={} m",
+      CORSIKA_LOG_TRACE("cannot extrapolate to points behind point of injection l={} m",
                         l / 1_m);
       if (throw_) {
         throw std::runtime_error(
@@ -75,7 +81,7 @@ namespace corsika {
     }
 
     if (upper >= X_.size()) {
-      CORSIKA_LOG_ERROR(
+      CORSIKA_LOG_TRACE(
           "shower axis too short, cannot extrapolate (l / max_length_ = {} )",
           l / max_length_);
       if (throw_) {

corsika/detail/media/SlidingPlanarExponential.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
 
-#include <corsika/media/SlidingPlanarExponential.hpp>
-
 namespace corsika {
 
-  template <typename T>
-  inline SlidingPlanarExponential<T>::SlidingPlanarExponential(
-      Point const& p0, MassDensityType rho0, LengthType lambda,
-      NuclearComposition const& nuclComp, LengthType referenceHeight)
-      : BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
-      , nuclComp_(nuclComp)
-      , referenceHeight_(referenceHeight) {}
+  template <typename TDerived>
+  inline SlidingPlanarExponential<TDerived>::SlidingPlanarExponential(
+      Point const& p0, MassDensityType const rho0, LengthType const lambda,
+      NuclearComposition const& nuclComp, LengthType const referenceHeight)
+      : BaseExponential<SlidingPlanarExponential<TDerived>>(p0, referenceHeight, rho0,
+                                                            lambda)
+      , nuclComp_(nuclComp) {}
 
-  template <typename T>
-  inline MassDensityType SlidingPlanarExponential<T>::getMassDensity(
+  template <typename TDerived>
+  inline MassDensityType SlidingPlanarExponential<TDerived>::getMassDensity(
       Point const& point) const {
-    auto const height =
-        (point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-            .getNorm() -
-        referenceHeight_;
-    return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
-           exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * height);
+    auto const heightFull =
+        (point - BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .getNorm();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getMassDensity(
+        heightFull);
   }
 
-  template <typename T>
-  inline NuclearComposition const& SlidingPlanarExponential<T>::getNuclearComposition()
-      const {
+  template <typename TDerived>
+  inline NuclearComposition const&
+  SlidingPlanarExponential<TDerived>::getNuclearComposition() const {
     return nuclComp_;
   }
 
-  template <typename T>
-  inline GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
-      BaseTrajectory const& traj, LengthType l) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, l,
-                                                                               axis);
+  template <typename TDerived>
+  inline GrammageType SlidingPlanarExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getIntegratedGrammage(
+        traj, axis);
   }
 
-  template <typename T>
-  inline LengthType SlidingPlanarExponential<T>::getArclengthFromGrammage(
+  template <typename TDerived>
+  inline LengthType SlidingPlanarExponential<TDerived>::getArclengthFromGrammage(
       BaseTrajectory const& traj, GrammageType const grammage) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getArclengthFromGrammage(
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getArclengthFromGrammage(
         traj, grammage, axis);
   }
 

corsika/detail/media/SlidingPlanarTabular.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
+
+namespace corsika {
+
+  template <typename T>
+  inline SlidingPlanarTabular<T>::SlidingPlanarTabular(
+      Point const& p0, std::function<MassDensityType(LengthType)> const& rho,
+      unsigned int const nBins, LengthType const deltaHeight,
+      NuclearComposition const& nuclComp, LengthType referenceHeight)
+      : BaseTabular<SlidingPlanarTabular<T>>(p0, referenceHeight, rho, nBins, deltaHeight)
+      , nuclComp_(nuclComp) {}
+
+  template <typename T>
+  inline MassDensityType SlidingPlanarTabular<T>::getMassDensity(
+      Point const& point) const {
+    auto const heightFull =
+        (point - BaseTabular<SlidingPlanarTabular<T>>::getAnchorPoint()).getNorm();
+    return BaseTabular<SlidingPlanarTabular<T>>::getMassDensity(heightFull);
+  }
+
+  template <typename T>
+  inline NuclearComposition const& SlidingPlanarTabular<T>::getNuclearComposition()
+      const {
+    return nuclComp_;
+  }
+
+  template <typename T>
+  inline GrammageType SlidingPlanarTabular<T>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    return BaseTabular<SlidingPlanarTabular<T>>::getIntegratedGrammage(traj);
+  }
+
+  template <typename T>
+  inline LengthType SlidingPlanarTabular<T>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+    return BaseTabular<SlidingPlanarTabular<T>>::getArclengthFromGrammage(traj, grammage);
+  }
+
+} // namespace corsika

corsika/detail/media/UniformRefractiveIndex.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -24,7 +23,7 @@ namespace corsika {
   }
 
   template <typename T>
-  inline void UniformRefractiveIndex<T>::setRefractiveIndex(double const& n) {
+  inline void UniformRefractiveIndex<T>::setRefractiveIndex(double const n) {
     n_ = n;
   }
 

corsika/detail/media/Universe.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/detail/media/VolumeTreeNode.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -52,13 +51,31 @@ namespace corsika {
     }
   }
 
+  template <typename IModelProperties>
+  inline VolumeTreeNode<IModelProperties>*
+  VolumeTreeNode<IModelProperties>::getContainingNode(Point const& p) {
+    // see Scott Meyers, Effective C++ 3rd ed., Item 3
+    return const_cast<VTN_type*>(std::as_const(*this).getContainingNode(p));
+  }
+
+  template <typename IModelProperties>
+  inline void VolumeTreeNode<IModelProperties>::addChildToContainingNode(Point const& p,
+                                                                         VTNUPtr pChild) {
+    VolumeTreeNode<IModelProperties>* node = getContainingNode(p);
+    if (!node) {
+      CORSIKA_LOG_ERROR("Adding child at {} failed!. No containing node", p);
+      throw std::runtime_error("Failed adding child node. No parent at chosen location");
+    }
+    node->addChild(std::move(pChild));
+  }
+
   template <typename IModelProperties>
   template <typename TCallable, bool preorder>
-  inline void VolumeTreeNode<IModelProperties>::walk(TCallable func) {
+  inline void VolumeTreeNode<IModelProperties>::walk(TCallable func) const {
     if constexpr (preorder) { func(*this); }
 
     std::for_each(childNodes_.begin(), childNodes_.end(),
-                  [&](auto& v) { v->walk(func); });
+                  [&](auto const& v) { v->walk(func); });
 
     if constexpr (!preorder) { func(*this); };
   }

corsika/detail/media/WeightProvider.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <vector>
-
-namespace corsika {
-
-  template <class AConstIterator, class BConstIterator>
-  inline WeightProviderIterator<AConstIterator, BConstIterator>::WeightProviderIterator(
-      AConstIterator a, BConstIterator b)
-      : aIter_(a)
-      , bIter_(b) {}
-
-  template <class AConstIterator, class BConstIterator>
-  inline typename WeightProviderIterator<AConstIterator, BConstIterator>::value_type
-      WeightProviderIterator<AConstIterator, BConstIterator>::operator*() const {
-    return ((*aIter_) * (*bIter_)).magnitude();
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  inline WeightProviderIterator<AConstIterator, BConstIterator>&
-  WeightProviderIterator<AConstIterator,
-                         BConstIterator>::operator++() { // prefix ++
-    ++aIter_;
-    ++bIter_;
-    return *this;
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  inline bool WeightProviderIterator<AConstIterator, BConstIterator>::operator==(
-      WeightProviderIterator other) {
-    return aIter_ == other.aIter_;
-  }
-
-  template <class AConstIterator, class BConstIterator>
-  inline bool WeightProviderIterator<AConstIterator, BConstIterator>::operator!=(
-      WeightProviderIterator other) {
-    return !(*this == other);
-  }
-
-} // namespace corsika

corsika/detail/modules/HadronicElasticModel.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once

corsika/detail/modules/LongitudinalProfile.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/Logging.hpp>
-
+#include <corsika/framework/core/Step.hpp>
 #include <corsika/modules/LongitudinalProfile.hpp>
 
 #include <cmath>
 #include <iomanip>
 #include <limits>
+#include <utility>
 
 namespace corsika {
 
-  inline LongitudinalProfile::LongitudinalProfile(ShowerAxis const& shower_axis,
-                                                  GrammageType dX)
-      : dX_(dX)
-      , shower_axis_{shower_axis}
-      , profiles_{static_cast<unsigned int>(shower_axis.getMaximumX() / dX_) + 1} {}
-
-  template <typename TParticle, typename TTrack>
-  inline ProcessReturn LongitudinalProfile::doContinuous(TParticle const& vP,
-                                                         TTrack const& vTrack,
-                                                         bool const) {
-    auto const pid = vP.getPID();
-
-    GrammageType const grammageStart = shower_axis_.getProjectedX(vTrack.getPosition(0));
-    GrammageType const grammageEnd = shower_axis_.getProjectedX(vTrack.getPosition(1));
-
-    CORSIKA_LOG_DEBUG("longprof: pos1={} m, pos2={}, X1={} g/cm2, X2={} g/cm2",
-                      vTrack.getPosition(0).getCoordinates() / 1_m,
-                      vTrack.getPosition(1).getCoordinates() / 1_m,
-                      grammageStart / 1_g * square(1_cm),
-                      grammageEnd / 1_g * square(1_cm));
+  template <typename TOutput>
+  template <typename... TArgs>
+  inline LongitudinalProfile<TOutput>::LongitudinalProfile(TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...) {}
 
-    // Note: particle may go also "upward", thus, grammageEnd<grammageStart
-    int const binStart = std::ceil(grammageStart / dX_);
-    int const binEnd = std::floor(grammageEnd / dX_);
-
-    for (int b = binStart; b <= binEnd; ++b) {
-      if (pid == Code::Photon) {
-        profiles_.at(b)[ProfileIndex::Photon]++;
-      } else if (pid == Code::Positron) {
-        profiles_.at(b)[ProfileIndex::Positron]++;
-      } else if (pid == Code::Electron) {
-        profiles_.at(b)[ProfileIndex::Electron]++;
-      } else if (pid == Code::MuPlus) {
-        profiles_.at(b)[ProfileIndex::MuPlus]++;
-      } else if (pid == Code::MuMinus) {
-        profiles_.at(b)[ProfileIndex::MuMinus]++;
-      } else if (is_hadron(pid)) {
-        profiles_.at(b)[ProfileIndex::Hadron]++;
-      }
-    }
+  template <typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn LongitudinalProfile<TOutput>::doContinuous(
+      Step<TParticle> const& step, bool const) {
 
+    auto const pid = step.getParticlePre().getPID();
+    this->write(step.getPositionPre(), step.getPositionPost(), pid,
+                step.getParticlePre().getWeight()); // weight hardcoded so far
     return ProcessReturn::Ok;
   }
 
-  inline void LongitudinalProfile::save(std::string const& filename, const int width,
-                                        const int precision) {
-    CORSIKA_LOG_DEBUG("Write longprof to {}", filename);
-    std::ofstream f{filename};
-    f << "# X / g·cm¯², photon, e+, e-, mu+, mu-, all hadrons" << std::endl;
-    for (size_t b = 0; b < profiles_.size(); ++b) {
-      f << std::setprecision(5) << std::setw(11) << b * (dX_ / (1_g / 1_cm / 1_cm));
-      for (auto const& N : profiles_.at(b)) {
-        f << std::setw(width) << std::setprecision(precision) << std::scientific << N;
-      }
-      f << std::endl;
-    }
+  template <typename TOutput>
+  inline YAML::Node LongitudinalProfile<TOutput>::getConfig() const {
+    YAML::Node node;
+    node["type"] = "LongitudinalProfile";
+
+    return node;
   }
+
 } // namespace corsika

corsika/detail/modules/ObservationPlane.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 namespace corsika {
 
   template <typename TTracking, typename TOutput>
+  template <typename... TArgs>
   ObservationPlane<TTracking, TOutput>::ObservationPlane(Plane const& obsPlane,
                                                          DirectionVector const& x_axis,
-                                                         bool deleteOnHit)
-      : plane_(obsPlane)
-      , deleteOnHit_(deleteOnHit)
-      , energy_ground_(0_GeV)
-      , count_ground_(0)
+                                                         bool const deleteOnHit,
+                                                         TArgs&&... args)
+      : TOutput(std::forward<TArgs>(args)...)
+      , plane_(obsPlane)
       , xAxis_(x_axis.normalized())
-      , yAxis_(obsPlane.getNormal().cross(xAxis_)) {}
+      , yAxis_(obsPlane.getNormal().cross(xAxis_))
+      , deleteOnHit_(deleteOnHit) {}
 
   template <typename TTracking, typename TOutput>
-  template <typename TParticle, typename TTrajectory>
+  template <typename TParticle>
   inline ProcessReturn ObservationPlane<TTracking, TOutput>::doContinuous(
-      TParticle& particle, TTrajectory&, bool const stepLimit) {
+      Step<TParticle>& step, bool const stepLimit) {
     /*
        The current step did not yet reach the ObservationPlane, do nothing now and wait:
      */
     if (!stepLimit) {
-#ifdef DEBUG
+      // @todo this is actually needed to fix small instabilities of the leap-frog
+      // tracking: Note, this is NOT a general solution and should be clearly revised with
+      // a more robust tracking. #ifdef DEBUG
       if (deleteOnHit_) {
+        // since this is basically a bug, it cannot be tested LCOV_EXCL_START
         LengthType const check =
-            (particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal());
+            (step.getPositionPost() - plane_.getCenter()).dot(plane_.getNormal());
         if (check < 0_m) {
-          CORSIKA_LOG_DEBUG("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
-        }
-      }
-#endif
-      return ProcessReturn::Ok;
+          CORSIKA_LOG_WARN("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
+          CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
+        } else
+          return ProcessReturn::Ok;
+        // LCOV_EXCL_STOP
+      } else
+        // #endif
+        return ProcessReturn::Ok;
     }
 
-    HEPEnergyType const energy = particle.getEnergy();
-    Point const pointOfIntersection = particle.getPosition();
+    HEPEnergyType const kineticEnergy = step.getEkinPost();
+    Point const pointOfIntersection = step.getPositionPost();
     Vector const displacement = pointOfIntersection - plane_.getCenter();
+    DirectionVector const direction = step.getDirectionPost();
 
     // add our particles to the output file stream
-    this->write(particle.getPID(), energy, displacement.dot(xAxis_),
-                displacement.dot(yAxis_));
+    double const weight = step.getParticlePre().getWeight();
+    this->write(step.getParticlePre().getPID(), kineticEnergy, displacement.dot(xAxis_),
+                displacement.dot(yAxis_), 0_m, direction.dot(xAxis_),
+                direction.dot(yAxis_), direction.dot(plane_.getNormal()),
+                step.getTimePost(), weight);
 
     CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
 
     if (deleteOnHit_) {
-      count_ground_++;
-      energy_ground_ += energy;
       return ProcessReturn::ParticleAbsorbed;
     } else {
       return ProcessReturn::Ok;
     }
-  }
+  } // namespace corsika
 
   template <typename TTracking, typename TOutput>
   template <typename TParticle, typename TTrajectory>
   inline LengthType ObservationPlane<TTracking, TOutput>::getMaxStepLength(
       TParticle const& particle, TTrajectory const& trajectory) {
 
-    CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(),
-                      particle.getPosition(), particle.getDirection(), plane_.asString());
+    CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, plane={}",
+                      particle.asString(), particle.getPosition(),
+                      particle.getDirection(), plane_.asString());
 
     auto const intersection = TTracking::intersect(particle, plane_);
 
     TimeType const timeOfIntersection = intersection.getEntry();
     CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
-    if (timeOfIntersection < TimeType::zero()) {
+    if (timeOfIntersection <= TimeType::zero()) {
       return std::numeric_limits<double>::infinity() * 1_m;
     }
     if (timeOfIntersection > trajectory.getDuration()) {
       return std::numeric_limits<double>::infinity() * 1_m;
     }
     double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
-    auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
-    auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm();
-    CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m", dist / 1_m);
-    return dist;
-  }
-
-  template <typename TTracking, typename TOutput>
-  inline void ObservationPlane<TTracking, TOutput>::showResults() const {
-    CORSIKA_LOG_INFO(
-        " ******************************\n"
-        " ObservationPlane: \n"
-        " energy an ground (GeV)     :  {}\n"
-        " no. of particles at ground :  {}\n"
-        " ******************************",
-        energy_ground_ / 1_GeV, count_ground_);
+    CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength dist={} m, pos={}",
+                      trajectory.getLength(fractionOfIntersection) / 1_m,
+                      trajectory.getPosition(fractionOfIntersection));
+    return trajectory.getLength(fractionOfIntersection);
   }
 
   template <typename TTracking, typename TOutput>
@@ -103,7 +101,9 @@ namespace corsika {
 
     // basic info
     node["type"] = "ObservationPlane";
-    node["units"] = "m"; // add default units for values
+    node["units"]["length"] = "m"; // add default units for values
+    node["units"]["energy"] = "GeV";
+    node["units"]["time"] = "s";
 
     // the center of the plane
     auto const center{plane_.getCenter()};
@@ -137,10 +137,4 @@ namespace corsika {
     return node;
   }
 
-  template <typename TTracking, typename TOutput>
-  inline void ObservationPlane<TTracking, TOutput>::reset() {
-    energy_ground_ = 0_GeV;
-    count_ground_ = 0;
-  }
-
 } // namespace corsika

corsika/detail/modules/ObservationVolume.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.
+ */
+
+namespace corsika {
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  ObservationVolume<TTracking, TVolume, TOutput>::ObservationVolume(TVolume vol)
+      : vol_(vol)
+      , energy_(0_GeV)
+      , count_(0) {}
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  template <typename TParticle>
+  inline ProcessReturn ObservationVolume<TTracking, TVolume, TOutput>::doContinuous(
+      Step<TParticle>& step, bool const stepLimit) {
+    /*
+       The current step did not yet reach the ObservationVolume, do nothing now and
+       wait:
+     */
+    if (!stepLimit) { return ProcessReturn::Ok; }
+
+    HEPEnergyType const kineticEnergy = step.getEkinPost();
+    Point const pointOfIntersection = step.getPositionPost();
+    DirectionVector const dirction = step.getDirectionPost();
+    double const weight = step.getParticlePre().getWeight();
+
+    // add particles to the output file stream
+    auto cs = vol_.getCoordinateSystem();
+    this->write(step.getParticlePre().getPID(), kineticEnergy,
+                pointOfIntersection.getX(cs), pointOfIntersection.getY(cs),
+                pointOfIntersection.getZ(cs), dirction.getX(cs), dirction.getY(cs),
+                dirction.getZ(cs), step.getTimePost(), weight);
+
+    // always absorb particles
+    count_++;
+    energy_ += kineticEnergy + get_mass(step.getParticlePre().getPID());
+    return ProcessReturn::ParticleAbsorbed;
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  template <typename TParticle, typename TTrajectory>
+  inline LengthType ObservationVolume<TTracking, TVolume, TOutput>::getMaxStepLength(
+      TParticle const& particle, TTrajectory const& trajectory) {
+
+    CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, Box={}",
+                      particle.asString(), particle.getPosition(),
+                      particle.getDirection(), vol_.asString());
+
+    auto const intersection = TTracking::intersect(particle, vol_);
+
+    TimeType const timeOfEntry = intersection.getEntry();
+    TimeType const timeOfExit = intersection.getExit();
+    CORSIKA_LOG_TRACE("timeOfEntry={}, timeOfExit={}", timeOfEntry, timeOfExit);
+    if (timeOfEntry < TimeType::zero()) {
+      if (timeOfExit < TimeType::zero()) {
+        // opposite direction
+        return std::numeric_limits<double>::infinity() * 1_m;
+      } else {
+        // inside box: not zero but 1 pm, to allow short lived particles to decay
+        return 1e-12 * 1_m;
+      }
+    }
+    if (timeOfEntry > trajectory.getDuration()) {
+      // can not reach
+      return std::numeric_limits<double>::infinity() * 1_m;
+    }
+    double const fractionOfIntersection = timeOfEntry / trajectory.getDuration();
+    CORSIKA_LOG_TRACE("ObservationVolume: getMaxStepLength dist={} m, pos={}",
+                      trajectory.getLength(fractionOfIntersection) / 1_m,
+                      trajectory.getPosition(fractionOfIntersection));
+    return trajectory.getLength(fractionOfIntersection);
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  inline void ObservationVolume<TTracking, TVolume, TOutput>::showResults() const {
+    CORSIKA_LOG_INFO(
+        " ******************************\n"
+        " ObservationVolume: \n"
+        " energy at Box (GeV)     :  {}\n"
+        " no. of particles at Box :  {}\n"
+        " ******************************",
+        energy_ / 1_GeV, count_);
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  inline YAML::Node ObservationVolume<TTracking, TVolume, TOutput>::getConfig() const {
+    using namespace units::si;
+    auto cs = vol_.getCoordinateSystem();
+    auto center = vol_.getCenter();
+
+    // construct the top-level node
+    YAML::Node node;
+
+    // basic info
+    node["type"] = "ObservationVolume";
+    node["units"]["length"] = "m";
+    node["units"]["energy"] = "GeV";
+    node["units"]["time"] = "s";
+
+    // save each component in its native coordinate system
+    auto const root_cs = get_root_CoordinateSystem();
+    node["center"].push_back(center.getX(root_cs) / 1_m);
+    node["center"].push_back(center.getY(root_cs) / 1_m);
+    node["center"].push_back(center.getZ(root_cs) / 1_m);
+
+    // the x-axis vector
+    DirectionVector const x_axis = DirectionVector{cs, {1, 0, 0}};
+    node["x-axis"].push_back(x_axis.getX(root_cs).magnitude());
+    node["x-axis"].push_back(x_axis.getY(root_cs).magnitude());
+    node["x-axis"].push_back(x_axis.getZ(root_cs).magnitude());
+
+    // the y-axis vector
+    DirectionVector const y_axis = DirectionVector{cs, {0, 1, 0}};
+    node["y-axis"].push_back(y_axis.getX(root_cs).magnitude());
+    node["y-axis"].push_back(y_axis.getY(root_cs).magnitude());
+    node["y-axis"].push_back(y_axis.getZ(root_cs).magnitude());
+
+    // the x-axis vector
+    DirectionVector const z_axis = DirectionVector{cs, {0, 0, 1}};
+    node["z-axis"].push_back(z_axis.getX(root_cs).magnitude());
+    node["z-axis"].push_back(z_axis.getY(root_cs).magnitude());
+    node["z-axis"].push_back(z_axis.getZ(root_cs).magnitude());
+
+    return node;
+  }
+
+  template <typename TTracking, typename TVolume, typename TOutput>
+  inline void ObservationVolume<TTracking, TVolume, TOutput>::reset() {
+    energy_ = 0_GeV;
+    count_ = 0;
+  }
+
+} // namespace corsika