corsika/detail/modules/radio/ZHS.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/ZHS.hpp>
+
+namespace corsika {
+
+  template <typename TRadioDetector, typename TPropagator>
+  template <typename Particle>
+  inline ProcessReturn ZHS<TRadioDetector, TPropagator>::simulate(
+      Step<Particle> const& step) {
+    auto const startTime{step.getTimePre()};
+    auto const endTime{step.getTimePost()};
+
+    // LCOV_EXCL_START
+    if (startTime == endTime) {
+      return ProcessReturn::Ok;
+      // LCOV_EXCL_STOP
+    } else {
+
+      auto const startPoint{step.getPositionPre()};
+      auto const endPoint{step.getPositionPost()};
+      LengthType const trackLength{(startPoint - endPoint).getNorm()};
+
+      auto const betaModule{(endPoint - startPoint).getNorm() /
+                            (constants::c * (endTime - startTime))};
+      auto const beta{(endPoint - startPoint).normalized() * betaModule};
+
+      auto const charge{get_charge(step.getParticlePre().getPID())};
+
+      // // get "mid" position of the track geometrically
+
+      auto const halfVector{(startPoint - endPoint) / 2};
+      auto const midPoint{endPoint + halfVector};
+
+      // get thinning weight
+      auto const thinningWeight{step.getParticlePre().getWeight()};
+
+      auto const constants{charge * emConstant_ * thinningWeight};
+
+      // we loop over each observer in the collection
+      for (auto& observer : observers_.getObservers()) {
+        auto midPaths{this->propagator_.propagate(step.getParticlePre(), midPoint,
+                                                  observer.getLocation())};
+        // Loop over midPaths, first check Fraunhoffer limit
+        for (size_t i{0}; i < midPaths.size(); i++) {
+          double const uTimesK{beta.dot(midPaths[i].emit_) / betaModule};
+          double const sinTheta2{1. - uTimesK * uTimesK};
+          LengthType const lambda{constants::c / observer.getSampleRate()};
+          double const fraunhLimit{sinTheta2 * trackLength * trackLength /
+                                   midPaths[i].R_distance_ / lambda * 2 * M_PI};
+          // Checks if we are in fraunhoffer domain
+          if (fraunhLimit > 1.0) {
+            /// code for dividing track and calculating field.
+            double const nSubTracks{sqrt(fraunhLimit) + 1};
+            auto const step_{(endPoint - startPoint) / nSubTracks};
+            TimeType const timeStep{(endTime - startTime) / nSubTracks};
+            // energy should be divided up when it is possible to get the energy at end of
+            // track!!!!
+            auto point1{startPoint};
+            TimeType time1{startTime};
+            for (int j{0}; j < nSubTracks; j++) {
+              auto const point2{point1 + step_};
+              TimeType const time2{time1 + timeStep};
+              auto const newHalfVector{(point1 - point2) / 2.};
+              auto const newMidPoint{point2 + newHalfVector};
+              auto const newMidPaths{this->propagator_.propagate(
+                  step.getParticlePre(), newMidPoint, observer.getLocation())};
+              // A function for calculating the field should be made since it is repeated
+              // later
+              for (size_t k{0}; k < newMidPaths.size(); k++) {
+
+                double const n_source{newMidPaths[k].refractive_index_source_};
+
+                double const betaTimesK{beta.dot(newMidPaths[k].emit_)};
+                TimeType const midTime{(time1 + time2) / 2.};
+                TimeType detectionTime1{time1 + newMidPaths[k].propagation_time_ -
+                                        n_source * betaTimesK * (time1 - midTime)};
+                TimeType detectionTime2{time2 + newMidPaths[k].propagation_time_ -
+                                        n_source * betaTimesK * (time2 - midTime)};
+
+                // make detectionTime1_ be the smallest time =>
+                // changes step function order so sign is changed to account for it
+                double sign = 1.;
+                if (detectionTime1 > detectionTime2) {
+                  detectionTime1 = time2 + newMidPaths[k].propagation_time_ -
+                                   n_source * betaTimesK * (time2 - midTime);
+                  detectionTime2 = time1 + newMidPaths[k].propagation_time_ -
+                                   n_source * betaTimesK * (time1 - midTime);
+                  sign = -1.;
+                } // end if statement for time structure
+
+                double const startBin{
+                    std::floor((detectionTime1 - observer.getStartTime()) *
+                                   observer.getSampleRate() +
+                               0.5)};
+                double const endBin{
+                    std::floor((detectionTime2 - observer.getStartTime()) *
+                                   observer.getSampleRate() +
+                               0.5)};
+
+                auto const betaPerp{
+                    newMidPaths[k].emit_.cross(beta.cross(newMidPaths[k].emit_))};
+                double const denominator{1. - n_source * betaTimesK};
+
+                if (startBin == endBin) {
+                  // track contained in bin
+                  // if not in Cerenkov angle then
+                  if (std::fabs(denominator) > 1.e-15) {
+                    double const f{
+                        std::fabs((detectionTime2 * observer.getSampleRate() -
+                                   detectionTime1 * observer.getSampleRate()))};
+                    VectorPotential const Vp = betaPerp * sign * constants * f /
+                                               denominator / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime2, betaPerp, Vp);
+                  } else { // If emission in Cerenkov angle => approximation
+                    double const f{time2 * observer.getSampleRate() -
+                                   time1 * observer.getSampleRate()};
+                    VectorPotential const Vp =
+                        betaPerp * sign * constants * f / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime2, betaPerp, Vp);
+                  } // end if Cerenkov angle approx
+                } else {
+                  /*Track is contained in more than one bin*/
+                  int const numberOfBins{static_cast<int>(endBin - startBin)};
+                  // first contribution/ plus 1 bin minus 0.5 from new observer ruonding
+                  double f{std::fabs(startBin + 0.5 -
+                                     (detectionTime1 - observer.getStartTime()) *
+                                         observer.getSampleRate())};
+                  VectorPotential Vp = betaPerp * sign * constants * f / denominator /
+                                       newMidPaths[k].R_distance_;
+                  observer.receive(detectionTime1, betaPerp, Vp);
+                  // intermidiate contributions
+                  for (int it{1}; it < numberOfBins; ++it) {
+                    Vp = betaPerp * constants / denominator / newMidPaths[k].R_distance_;
+                    observer.receive(detectionTime1 + static_cast<double>(it) /
+                                                          observer.getSampleRate(),
+                                     betaPerp, Vp);
+                  } // end loop over bins in which potential vector is not zero
+                  // final contribution// f +0.5 from new observer rounding
+                  f = std::fabs((detectionTime2 - observer.getStartTime()) *
+                                    observer.getSampleRate() +
+                                0.5 - endBin);
+                  Vp = betaPerp * sign * constants * f / denominator /
+                       newMidPaths[k].R_distance_;
+                  observer.receive(detectionTime2, betaPerp, Vp);
+                } // end if statement for track in multiple bins
+
+              } // end of loop over newMidPaths
+              // update points for next sub track
+              point1 = point1 + step_;
+              time1 = time1 + timeStep;
+            }
+
+          } else // Calculate vector potential of whole track
+          {
+            double const n_source{midPaths[i].refractive_index_source_};
+
+            double const betaTimesK{beta.dot(midPaths[i].emit_)};
+            TimeType const midTime{(startTime + endTime) / 2};
+            TimeType detectionTime1{startTime + midPaths[i].propagation_time_ -
+                                    n_source * betaTimesK * (startTime - midTime)};
+            TimeType detectionTime2{endTime + midPaths[i].propagation_time_ -
+                                    n_source * betaTimesK * (endTime - midTime)};
+
+            // make detectionTime1_ be the smallest time =>
+            // changes step function order so sign is changed to account for it
+            double sign = 1.;
+            if (detectionTime1 > detectionTime2) {
+              detectionTime1 = endTime + midPaths[i].propagation_time_ -
+                               n_source * betaTimesK * (endTime - midTime);
+              detectionTime2 = startTime + midPaths[i].propagation_time_ -
+                               n_source * betaTimesK * (startTime - midTime);
+              sign = -1.;
+            } // end if statement for time structure
+
+            double const startBin{std::floor((detectionTime1 - observer.getStartTime()) *
+                                                 observer.getSampleRate() +
+                                             0.5)};
+            double const endBin{std::floor((detectionTime2 - observer.getStartTime()) *
+                                               observer.getSampleRate() +
+                                           0.5)};
+
+            auto const betaPerp{midPaths[i].emit_.cross(beta.cross(midPaths[i].emit_))};
+            double const denominator{1. -
+                                     midPaths[i].refractive_index_source_ * betaTimesK};
+
+            if (startBin == endBin) {
+              // track contained in bin
+              // if not in Cerenkov angle then
+              if (std::fabs(denominator) > 1.e-15) {
+                double const f{std::fabs((detectionTime2 * observer.getSampleRate() -
+                                          detectionTime1 * observer.getSampleRate()))};
+
+                VectorPotential const Vp = betaPerp * sign * constants * f / denominator /
+                                           midPaths[i].R_distance_;
+                observer.receive(detectionTime2, betaPerp, Vp);
+              } else { // If emission in Cerenkov angle => approximation
+                double const f{endTime * observer.getSampleRate() -
+                               startTime * observer.getSampleRate()};
+                VectorPotential const Vp =
+                    betaPerp * sign * constants * f / midPaths[i].R_distance_;
+                observer.receive(detectionTime2, betaPerp, Vp);
+              } // end if Cerenkov angle approx
+            } else {
+              /*Track is contained in more than one bin*/
+              int const numberOfBins{static_cast<int>(endBin - startBin)};
+              // TODO: should we check for Cerenkov angle?
+              // first contribution
+              double f{std::fabs(startBin + 0.5 -
+                                 (detectionTime1 - observer.getStartTime()) *
+                                     observer.getSampleRate())};
+              VectorPotential Vp =
+                  betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+              observer.receive(detectionTime1, betaPerp, Vp);
+              // intermediate contributions
+              for (int it{1}; it < numberOfBins; ++it) {
+                Vp = betaPerp * sign * constants / denominator / midPaths[i].R_distance_;
+                observer.receive(
+                    detectionTime1 + static_cast<double>(it) / observer.getSampleRate(),
+                    betaPerp, Vp);
+              } // end loop over bins in which potential vector is not zero
+              // final contribution
+              f = std::fabs((detectionTime2 - observer.getStartTime()) *
+                                observer.getSampleRate() +
+                            0.5 - endBin);
+              Vp =
+                  betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
+              observer.receive(detectionTime2, betaPerp, Vp);
+            } // end if statement for track in multiple bins
+
+          } // finish if statement of track in fraunhoffer or not
+
+        } // end loop over mid paths
+
+      } // END: loop over observers
+      return ProcessReturn::Ok;
+    }
+  } // end simulate
+
+} // namespace corsika

corsika/detail/modules/radio/detectors/ObserverCollection.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/detectors/ObserverCollection.hpp>
+
+namespace corsika {
+
+  template <typename TObserverImpl>
+  inline void ObserverCollection<TObserverImpl>::addObserver(
+      TObserverImpl const& observer) {
+    observers_.push_back(observer);
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl& ObserverCollection<TObserverImpl>::at(std::size_t const i) {
+    return observers_.at(i);
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl const& ObserverCollection<TObserverImpl>::at(
+      std::size_t const i) const {
+    return observers_.at(i);
+  }
+
+  template <typename TObserverImpl>
+  inline int ObserverCollection<TObserverImpl>::size() const {
+    return observers_.size();
+  }
+
+  template <typename TObserverImpl>
+  inline std::vector<TObserverImpl>& ObserverCollection<TObserverImpl>::getObservers() {
+    return observers_;
+  }
+
+  template <typename TObserverImpl>
+  inline void ObserverCollection<TObserverImpl>::reset() {
+    std::for_each(observers_.begin(), observers_.end(),
+                  std::mem_fn(&TObserverImpl::reset));
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/observers/Observer.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/observers/Observer.hpp>
+
+namespace corsika {
+
+  template <typename TObserverImpl>
+  inline Observer<TObserverImpl>::Observer(std::string const& name, Point const& location,
+                                           CoordinateSystemPtr const& coordinateSystem)
+      : name_(name)
+      , location_(location)
+      , coordinateSystem_(coordinateSystem) {}
+
+  template <typename TObserverImpl>
+  inline Point const& Observer<TObserverImpl>::getLocation() const {
+    return location_;
+  }
+
+  template <typename TObserverImpl>
+  inline std::string const& Observer<TObserverImpl>::getName() const {
+    return name_;
+  }
+
+  template <typename TObserverImpl>
+  inline TObserverImpl& Observer<TObserverImpl>::implementation() {
+    return static_cast<TObserverImpl&>(*this);
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/observers/TimeDomainObserver.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/observers/TimeDomainObserver.hpp>
+
+namespace corsika {
+
+  inline TimeDomainObserver::TimeDomainObserver(std::string const& name,
+                                                Point const& location,
+                                                CoordinateSystemPtr coordinateSystem,
+                                                TimeType const& start_time,
+                                                TimeType const& duration,
+                                                InverseTimeType const& sample_rate,
+                                                TimeType const ground_hit_time)
+      : Observer(name, location, coordinateSystem)
+      , start_time_(start_time)
+      , duration_(std::abs(duration / 1_s) * 1_s)
+      , sample_rate_(std::abs(sample_rate / 1_Hz) * 1_Hz)
+      , ground_hit_time_(ground_hit_time)
+      , num_bins_(static_cast<std::size_t>(duration_ * sample_rate_ + 1.5l))
+      , waveformEX_(num_bins_, 0)
+      , waveformEY_(num_bins_, 0)
+      , waveformEZ_(num_bins_, 0)
+      , time_axis_(createTimeAxis()) {
+    if (0_s == duration_) {
+      CORSIKA_LOG_WARN(
+          "Observer: \"{}\" has a duration of zero. Nothing will be injected into it",
+          name);
+    } else if (duration_ != duration) {
+      CORSIKA_LOG_WARN(
+          "Observer: \"{}\" was given a negative duration. Set to absolute value.", name);
+    }
+
+    if (sample_rate_ != sample_rate) {
+      CORSIKA_LOG_WARN(
+          "Observer: \"{}\" was given a negative sampling rate. Set to absolute value.",
+          name);
+    }
+  }
+
+  inline void TimeDomainObserver::receive(
+      const TimeType time, [[maybe_unused]] const Vector<dimensionless_d>& receive_vector,
+      const ElectricFieldVector& efield) {
+
+    if (time < start_time_ || time > (start_time_ + duration_)) {
+      return;
+    } else {
+      // figure out the correct timebin to store the E-field value.
+      // NOTE: static cast is implicitly flooring
+      auto timebin{static_cast<std::size_t>(
+          std::floor((time - start_time_) * sample_rate_ + 0.5l))};
+
+      // store the x,y,z electric field components.
+      auto const& Electric_field_components{efield.getComponents(coordinateSystem_)};
+      waveformEX_.at(timebin) += (Electric_field_components.getX() * (1_m / 1_V));
+      waveformEY_.at(timebin) += (Electric_field_components.getY() * (1_m / 1_V));
+      waveformEZ_.at(timebin) += (Electric_field_components.getZ() * (1_m / 1_V));
+      // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+      // a 3D object
+    }
+  }
+
+  inline void TimeDomainObserver::receive(
+      const TimeType time, [[maybe_unused]] const Vector<dimensionless_d>& receive_vector,
+      const VectorPotential& vectorP) {
+
+    if (time < start_time_ || time > (start_time_ + duration_)) {
+      return;
+    } else {
+      // figure out the correct timebin to store the E-field value.
+      // NOTE: static cast is implicitly flooring
+      auto timebin{static_cast<std::size_t>(
+          std::floor((time - start_time_) * sample_rate_ + 0.5l))};
+
+      // store the x,y,z electric field components.
+      auto const& Vector_potential_components{vectorP.getComponents(coordinateSystem_)};
+      waveformEX_.at(timebin) +=
+          (Vector_potential_components.getX() * (1_m / (1_V * 1_s)));
+      waveformEY_.at(timebin) +=
+          (Vector_potential_components.getY() * (1_m / (1_V * 1_s)));
+      waveformEZ_.at(timebin) +=
+          (Vector_potential_components.getZ() * (1_m / (1_V * 1_s)));
+      // TODO: Check how they are stored in memory, row-wise or column-wise? Probably use
+      // a 3D object
+    }
+  }
+
+  inline auto const& TimeDomainObserver::getWaveformX() const { return waveformEX_; }
+
+  inline auto const& TimeDomainObserver::getWaveformY() const { return waveformEY_; }
+
+  inline auto const& TimeDomainObserver::getWaveformZ() const { return waveformEZ_; }
+
+  inline std::string const TimeDomainObserver::getDomainLabel() { return "Time"; }
+
+  inline std::vector<long double> TimeDomainObserver::createTimeAxis() const {
+
+    // create a 1-D xtensor to store time values so we can print them later.
+    std::vector<long double> times(num_bins_, 0);
+
+    // calculate the sample_period
+    auto sample_period{1 / sample_rate_};
+
+    // fill in every time-value
+    for (uint64_t i = 0; i < num_bins_; i++) {
+      // create the current time in nanoseconds
+      times.at(i) = static_cast<long double>(
+          ((start_time_ - ground_hit_time_) + i * sample_period) / 1_ns);
+    }
+
+    return times;
+  }
+
+  inline auto const TimeDomainObserver::getAxis() const { return time_axis_; }
+
+  inline InverseTimeType const& TimeDomainObserver::getSampleRate() const {
+    return sample_rate_;
+  }
+
+  inline TimeType const& TimeDomainObserver::getStartTime() const { return start_time_; }
+
+  inline void TimeDomainObserver::reset() {
+    std::fill(waveformEX_.begin(), waveformEX_.end(), 0);
+    std::fill(waveformEY_.begin(), waveformEY_.end(), 0);
+    std::fill(waveformEZ_.begin(), waveformEZ_.end(), 0);
+  }
+
+  inline YAML::Node TimeDomainObserver::getConfig() const {
+
+    // top-level config
+    YAML::Node config;
+
+    config["type"] = "TimeDomainObserver";
+    config["start time"] = start_time_ / 1_ns;
+    config["duration"] = duration_ / 1_ns;
+    config["number of bins"] = duration_ * sample_rate_;
+    config["sampling frequency"] = sample_rate_ / 1_GHz;
+
+    return config;
+  }
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/DummyTestPropagator.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/propagators/DummyTestPropagator.hpp>
+
+namespace corsika {
+
+  template <typename TEnvironment>
+  inline DummyTestPropagator<TEnvironment>::DummyTestPropagator(TEnvironment const& env)
+      : RadioPropagator<DummyTestPropagator, TEnvironment>(env) {
+    rindex.reserve(2);
+  }
+
+  template <typename TEnvironment>
+  template <typename Particle>
+  inline typename DummyTestPropagator<TEnvironment>::SignalPathCollection
+  DummyTestPropagator<TEnvironment>::propagate(Particle const& particle,
+                                               Point const& source,
+                                               Point const& destination) {
+
+    /**
+     * This is the simplest case of straight propagator
+     * where no integration takes place.
+     * This can be used for fast tests and checks of the radio module.
+     *
+     */
+
+    // these are used for the direction of emission and reception of signal at the
+    // observer
+    auto const emit_{(destination - source).normalized()};
+    auto const receive_{-emit_};
+
+    // the geometrical distance from the point of emission to an observer
+    auto const distance_{(destination - source).getNorm()};
+
+    // get the universe for this environment
+    auto const* const universe{Base::env_.getUniverse().get()};
+
+    // clear the refractive index vector and points deque for this signal propagation.
+    rindex.clear();
+    points.clear();
+
+    // get and store the refractive index of the first point 'source'.
+    auto const* const nodeSource{particle.getNode()};
+    auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
+    rindex.push_back(ri_source);
+    points.push_back(source);
+
+    // add the refractive index of last point 'destination' and store it.
+    auto const* const node{universe->getContainingNode(destination)};
+    auto const ri_destination{node->getModelProperties().getRefractiveIndex(destination)};
+    rindex.push_back(ri_destination);
+    points.push_back(destination);
+
+    // compute the average refractive index.
+    auto const averageRefractiveIndex_ = (ri_source + ri_destination) * 0.5;
+
+    // compute the total time delay.
+    TimeType const time = averageRefractiveIndex_ * (distance_ / constants::c);
+
+    return std::vector<SignalPath>(
+        1, SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
+                      receive_, distance_, points));
+
+  } // END: propagate()
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/NumericalIntegratingPropagator.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/propagators/NumericalIntegratingPropagator.hpp>
+
+namespace corsika {
+
+  template <typename TEnvironment>
+  // TODO: maybe the constructor doesn't take any arguments for the environment (?)
+  inline NumericalIntegratingPropagator<TEnvironment>::NumericalIntegratingPropagator(
+      TEnvironment const& env, LengthType const stepsize)
+      : RadioPropagator<NumericalIntegratingPropagator, TEnvironment>(env)
+      , stepsize_(stepsize) {}
+
+  template <typename TEnvironment>
+  template <typename Particle>
+  inline typename NumericalIntegratingPropagator<TEnvironment>::SignalPathCollection
+  NumericalIntegratingPropagator<TEnvironment>::propagate(
+      Particle const& particle, Point const& source, Point const& destination) const {
+
+    /*
+     * get the normalized (unit) vector from `source` to `destination'.
+     * this is also the `emit` and `receive` vectors in the SignalPath class.
+     * in this case emit and receive unit vectors should be the same
+     * so they are both called direction
+     */
+
+    // these are used for the direction of emission and reception of signal at the
+    // observer
+    auto const emit{(destination - source).normalized()};
+    auto const receive{-emit};
+
+    // the distance from the point of emission to an observer
+    auto const distance{(destination - source).getNorm()};
+
+    try {
+      if (stepsize_ <= 0.5 * distance) {
+
+        // "step" is the direction vector with length `stepsize`
+        auto const step{emit * stepsize_};
+
+        // calculate the number of points (roughly) for the numerical integration
+        auto const n_points{(destination - source).getNorm() / stepsize_};
+
+        // get the universe for this environment
+        auto const* const universe{Base::env_.getUniverse().get()};
+
+        // the points that we build along the way for the numerical integration
+        std::deque<Point> points;
+
+        // store value of the refractive index at points
+        std::vector<double> rindex;
+        rindex.reserve(n_points);
+
+        // get and store the refractive index of the first point 'source'
+        auto const* const nodeSource{particle.getNode()};
+        auto const ri_source{nodeSource->getModelProperties().getRefractiveIndex(source)};
+        rindex.push_back(ri_source);
+        points.push_back(source);
+
+        // loop from `source` to `destination` to store values before Simpson's rule.
+        // this loop skips the last point 'destination' and "misses" the extra point
+        for (auto point = source + step;
+             (point - destination).getNorm() > 0.6 * stepsize_; point = point + step) {
+
+          // get the environment node at this specific 'point'
+          auto const* const node{universe->getContainingNode(point)};
+
+          // get the associated refractivity at 'point'
+          auto const refractive_index{
+              node->getModelProperties().getRefractiveIndex(point)};
+          rindex.push_back(refractive_index);
+
+          // add this 'point' to our deque collection
+          points.push_back(point);
+        }
+
+        // Get the extra points that the for loop misses until the destination
+        auto const extrapoint_{points.back() + step};
+
+        // add the refractive index of last point 'destination' and store it
+        auto const* const node{universe->getContainingNode(destination)};
+        auto const ri_destination{
+            node->getModelProperties().getRefractiveIndex(destination)};
+        rindex.push_back(ri_destination);
+        points.push_back(destination);
+
+        auto N = rindex.size();
+        std::size_t index = 0;
+        double sum = rindex.at(index);
+        auto refra = rindex.at(index);
+        TimeType time{0_s};
+
+        if ((N - 1) % 2 == 0) {
+          // Apply the standard Simpson's rule
+          auto const h = ((destination - source).getNorm()) / (N - 1);
+
+          for (std::size_t index = 1; index < (N - 1); index += 2) {
+            sum += 4 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          for (std::size_t index = 2; index < (N - 1); index += 2) {
+            sum += 2 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          index = N - 1;
+          sum = sum + rindex.at(index);
+          refra += rindex.at(index);
+
+          // compute the total time delay.
+          time = sum * (h / (3 * constants::c));
+        } else {
+          // Apply Simpson's rule for one "extra" point and then subtract the difference
+          points.pop_back();
+          rindex.pop_back();
+          auto const* const node{universe->getContainingNode(extrapoint_)};
+          auto const ri_extrapoint{
+              node->getModelProperties().getRefractiveIndex(extrapoint_)};
+          rindex.push_back(ri_extrapoint);
+          points.push_back(extrapoint_);
+          auto const extrapoint2_{extrapoint_ + step};
+          auto const* const node2{universe->getContainingNode(extrapoint2_)};
+          auto const ri_extrapoint2{
+              node2->getModelProperties().getRefractiveIndex(extrapoint2_)};
+          rindex.push_back(ri_extrapoint2);
+          points.push_back(extrapoint2_);
+          N = rindex.size();
+          auto const h = ((extrapoint2_ - source).getNorm()) / (N - 1);
+          for (std::size_t index = 1; index < (N - 1); index += 2) {
+            sum += 4 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          for (std::size_t index = 2; index < (N - 1); index += 2) {
+            sum += 2 * rindex.at(index);
+            refra += rindex.at(index);
+          }
+          index = N - 1;
+          sum = sum + rindex.at(index);
+          refra += rindex.at(index);
+
+          // compute the total time delay including the correction
+          time =
+              sum * (h / (3 * constants::c)) -
+              (ri_extrapoint2 * ((extrapoint2_ - destination).getNorm()) / constants::c);
+        }
+
+        // compute the average refractive index.
+        auto const averageRefractiveIndex = refra / N;
+
+        return std::vector<SignalPath>(
+            1, SignalPath(time, averageRefractiveIndex, ri_source, ri_destination, emit,
+                          receive, distance, points));
+      } else {
+        throw stepsize_;
+      }
+    } catch (const LengthType& s) {
+      CORSIKA_LOG_ERROR("Please choose a smaller stepsize for the numerical integration");
+    }
+
+    std::deque<Point> const defaultPoints;
+    return std::vector<SignalPath>(
+        1, SignalPath(0_s, 0, 0, 0, emit, receive, distance,
+                      defaultPoints)); // Dummy return that is never called (combat
+                                       // compile warnings)
+  }                                    // END: propagate()
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/RadioPropagator.inl

+/*
+ * (c) Copyright 2022 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/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+  template <typename TImpl, typename TEnvironment>
+  inline RadioPropagator<TImpl, TEnvironment>::RadioPropagator(TEnvironment const& env)
+      : env_(env) {}
+
+} // namespace corsika

corsika/detail/modules/radio/propagators/SignalPath.inl

+/*
+ * (c) Copyright 2022 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 {
+
+  inline SignalPath::SignalPath(
+      TimeType const propagation_time, double const average_refractive_index,
+      double const refractive_index_source, double const refractive_index_destination,
+      Vector<dimensionless_d> const& emit, Vector<dimensionless_d> const& receive,
+      LengthType const R_distance, std::deque<Point> const& points)
+      : Path(points)
+      , propagation_time_(propagation_time)
+      , average_refractive_index_(average_refractive_index)
+      , refractive_index_source_(refractive_index_source)
+      , refractive_index_destination_(refractive_index_destination)
+      , emit_(emit)
+      , receive_(receive)
+      , R_distance_(R_distance) {}
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.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/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp>
+
+namespace corsika {
+
+  template <typename TEnvironment>
+  inline TabulatedFlatAtmospherePropagator<
+      TEnvironment>::TabulatedFlatAtmospherePropagator(TEnvironment const& env,
+                                                       Point const& upperLimit,
+                                                       Point const& lowerLimit,
+                                                       LengthType const step)
+      : RadioPropagator<TabulatedFlatAtmospherePropagator, TEnvironment>(env)
+      , upperLimit_(upperLimit)
+      , lowerLimit_(lowerLimit)
+      , step_(step)
+      , inverseStep_(1 / step)
+      , maxHeight_(upperLimit_.getCoordinates().getZ())
+      , minHeight_(lowerLimit_.getCoordinates().getZ() - 1_km) {
+    auto const minX_ = lowerLimit_.getCoordinates().getX();
+    auto const minY_ = lowerLimit_.getCoordinates().getY();
+    std::size_t const nBins_ = (maxHeight_ - minHeight_) * inverseStep_ + 1;
+    refractivityTable_.reserve(nBins_);
+    heightTable_.reserve(nBins_);
+    integratedRefractivityTable_.reserve(nBins_);
+
+    // get the root coordinate system of this environment
+    CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+    // get the universe for this environment
+    auto const* const universe{Base::env_.getUniverse().get()};
+
+    for (std::size_t i = 0; i < nBins_; i++) {
+      Point point_{rootCS, minX_, minY_, minHeight_ + i * step_};
+      auto const* const node{universe->getContainingNode(point_)};
+      auto const ri_ = node->getModelProperties().getRefractiveIndex(point_);
+      refractivityTable_.push_back(ri_ - 1);
+      auto const height_ = minHeight_ + i * step_;
+      heightTable_.push_back(height_);
+    }
+
+    double const stepOverMeter_{inverseStep_ * 1_m};
+    auto const intRerfZero_ = refractivityTable_.at(0) * stepOverMeter_;
+    integratedRefractivityTable_.push_back(intRerfZero_);
+    for (std::size_t i = 1; i < nBins_; i++) {
+      auto const intRefrI_ =
+          integratedRefractivityTable_[i - 1] + refractivityTable_[i] * stepOverMeter_;
+      integratedRefractivityTable_.push_back(intRefrI_);
+    }
+
+    // this is used to interpolate refractivity and integrated refractivity for particles
+    // below the lower limit
+    slopeRefrLower_ = (refractivityTable_.at(10) - refractivityTable_.front()) / 10.;
+    slopeIntRefrLower_ =
+        (integratedRefractivityTable_.at(10) - integratedRefractivityTable_.front()) /
+        10.;
+
+    // this is used to interpolate refractivity and integrated refractivity for particles
+    // above the upper limit
+    lastElement_ = integratedRefractivityTable_.size() - 1;
+    slopeRefrUpper_ =
+        (refractivityTable_.at(lastElement_) - refractivityTable_.at(lastElement_ - 10)) /
+        10.;
+    slopeIntRefrUpper_ = (integratedRefractivityTable_.at(lastElement_) -
+                          integratedRefractivityTable_.at(lastElement_ - 10)) /
+                         10.;
+  }
+
+  template <typename TEnvironment>
+  template <typename Particle>
+  inline typename TabulatedFlatAtmospherePropagator<TEnvironment>::SignalPathCollection
+  TabulatedFlatAtmospherePropagator<TEnvironment>::propagate(
+      [[maybe_unused]] Particle const& particle, Point const& source,
+      Point const& destination) {
+
+    /**
+     * This is a simple case of straight propagator where
+     * tabulated values of refractive index are called assuming
+     * a flat atmosphere.
+     *
+     */
+
+    // these are used for the direction of emission and reception of signal at the
+    // observer
+    auto const emit_{(destination - source).normalized()};
+    auto const receive_{-emit_};
+
+    // the geometrical distance from the point of emission to an observer
+    auto const distance_{(destination - source).getNorm()};
+
+    // clear the refractive index vector and points deque for this signal propagation.
+    rindex.clear();
+    points.clear();
+    auto const sourceZ_{source.getCoordinates().getZ()};
+    auto const checkMaxZ_{(maxHeight_ - sourceZ_) / 1_m};
+    auto ri_source{1.};
+    auto intRerfr_source{1.};
+    auto ri_destination{1.};
+    auto intRerfr_destination{1.};
+    auto height_{1.};
+
+    double const sourceHeight_{(sourceZ_ - heightTable_.front()) * inverseStep_};
+    double const destinationHeight_{
+        (destination.getCoordinates().getZ() - heightTable_.front()) * inverseStep_};
+
+    if ((sourceHeight_ + 0.5) >=
+        lastElement_) { // source particle is above maximum height
+      ri_source =
+          refractivityTable_.back() +
+          slopeRefrUpper_ * std::abs((sourceZ_ - heightTable_.back()) * inverseStep_) + 1;
+      intRerfr_source =
+          integratedRefractivityTable_.back() + slopeIntRefrUpper_ * std::abs(checkMaxZ_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      ri_destination = refractivityTable_.at(indexDestination_) + 1;
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = sourceHeight_ - destinationHeight_;
+
+    } else if ((sourceHeight_ + 0.5 < lastElement_) &&
+               sourceHeight_ > 0) { // source particle in the table
+      // get and store the refractive index of the first point 'source'.
+      std::size_t const indexSource_{static_cast<std::size_t>(sourceHeight_ + 0.5)};
+      ri_source = refractivityTable_.at(indexSource_) + 1;
+      intRerfr_source = integratedRefractivityTable_.at(indexSource_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      ri_destination = refractivityTable_.at(indexDestination_) + 1;
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ =
+          (heightTable_.at(indexSource_) - heightTable_.at(indexDestination_)) / 1_m;
+      if (height_ == 0) { height_ = 1.; }
+    } else if (sourceHeight_ ==
+               0) { // source particle is exactly at the lower edge of the table
+      // get and store the refractive index of the first point 'source'.
+      std::size_t const indexSource_{static_cast<std::size_t>(sourceHeight_ + 0.5)};
+      auto const ri_source{refractivityTable_.at(indexSource_) + 1};
+      intRerfr_source = integratedRefractivityTable_.at(indexSource_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = destinationHeight_ - sourceHeight_;
+    } else if (sourceHeight_ < 0) { // source particle is in the ground.
+      // get and store the refractive index of the first point 'source'.
+      ri_source =
+          refractivityTable_.front() + slopeRefrLower_ * std::abs(sourceHeight_) + 1;
+      intRerfr_source = integratedRefractivityTable_.front() +
+                        slopeIntRefrLower_ * std::abs(sourceHeight_);
+      rindex.push_back(ri_source);
+      points.push_back(source);
+
+      // add the refractive index of last point 'destination' and store it.
+      std::size_t const indexDestination_{
+          static_cast<std::size_t>(destinationHeight_ + 0.5)};
+      auto const ri_destination{refractivityTable_.at(indexDestination_) + 1};
+      intRerfr_destination = integratedRefractivityTable_.at(indexDestination_);
+      rindex.push_back(ri_destination);
+      points.push_back(destination);
+
+      height_ = destinationHeight_ - std::abs(sourceHeight_);
+    }
+
+    // compute the average refractive index.
+    auto const averageRefractiveIndex_ = (ri_source + ri_destination) * 0.5;
+
+    // compute the total time delay.
+    TimeType const time = (1 + (intRerfr_source - intRerfr_destination) / height_) *
+                          (distance_ / constants::c);
+
+    return {SignalPath(time, averageRefractiveIndex_, ri_source, ri_destination, emit_,
+                       receive_, distance_, points)};
+
+  } // END: propagate()
+
+} // namespace corsika

corsika/detail/modules/sibyll/Decay.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/modules/sibyll/Decay.hpp>
+#include <corsika/modules/sibyll/ParticleConversion.hpp>
+#include <corsika/modules/sibyll/SibStack.hpp>
+#include <corsika/modules/Random.hpp>
+
+#include <iostream>
+#include <vector>
+
+namespace corsika::sibyll {
+
+  inline Decay::Decay(const bool sibyll_printout_on)
+      : sibyll_listing_(sibyll_printout_on) {
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+    // switch off decays to avoid internal decay chains
+    setAllStable();
+    // handle all decays by default
+    handleAllDecays_ = true;
+  }
+
+  inline Decay::Decay(std::set<Code> const& vHandled)
+      : handleAllDecays_(false)
+      , handledDecays_(vHandled) {
+    setAllStable();
+  }
+
+  inline Decay::~Decay() {
+    CORSIKA_LOGGER_TRACE(logger_, "Total number of Sibyll decays n={}", count_);
+  }
+
+  inline bool Decay::canHandleDecay(const Code vParticleCode) {
+    // if known to sibyll and not proton or neutrino it can decay
+    if (is_nucleus(vParticleCode) || vParticleCode == Code::Proton ||
+        vParticleCode == Code::AntiProton || vParticleCode == Code::NuE ||
+        vParticleCode == Code::NuMu || vParticleCode == Code::NuTau ||
+        vParticleCode == Code::NuEBar || vParticleCode == Code::NuMuBar ||
+        vParticleCode == Code::NuTauBar || vParticleCode == Code::Electron ||
+        vParticleCode == Code::Positron)
+      return false;
+    else if (corsika::sibyll::convertToSibyllRaw(
+                 vParticleCode)) // non-zero for particles known to sibyll
+      return true;
+    else
+      return false;
+  }
+
+  inline void Decay::setHandleDecay(const Code vParticleCode) {
+    handleAllDecays_ = false;
+    CORSIKA_LOGGER_DEBUG(logger_, "set to handle decay of {}", vParticleCode);
+    if (Decay::canHandleDecay(vParticleCode))
+      handledDecays_.insert(vParticleCode);
+    else
+      throw std::runtime_error("this decay can not be handled by sibyll!");
+  }
+
+  inline void Decay::setHandleDecay(std::vector<Code> const& vParticleList) {
+    handleAllDecays_ = false;
+    for (auto const p : vParticleList) Decay::setHandleDecay(p);
+  }
+
+  inline bool Decay::isDecayHandled(corsika::Code const vParticleCode) {
+    if (handleAllDecays_ && Decay::canHandleDecay(vParticleCode))
+      return true;
+    else
+      return Decay::handledDecays_.find(vParticleCode) != Decay::handledDecays_.end()
+                 ? true
+                 : false;
+  }
+
+  inline void Decay::setStable(std::vector<Code> const& vParticleList) {
+    for (auto const p : vParticleList) Decay::setStable(p);
+  }
+
+  inline void Decay::setUnstable(std::vector<Code> const& vParticleList) {
+    for (auto const p : vParticleList) Decay::setUnstable(p);
+  }
+
+  inline bool Decay::isStable(Code const vCode) {
+    return abs(sibyll::convertToSibyllRaw(vCode)) <= 0 ? true : false;
+  }
+
+  inline bool Decay::isUnstable(Code const vCode) {
+    return abs(sibyll::convertToSibyllRaw(vCode)) > 0 ? true : false;
+  }
+
+  inline void Decay::setDecay(const Code vCode, const bool vMakeUnstable) {
+    vMakeUnstable ? setUnstable(vCode) : setStable(vCode);
+  }
+
+  inline void Decay::setUnstable(Code const vCode) {
+    CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"unstable\". ", vCode);
+
+    const int s_id = abs(sibyll::convertToSibyllRaw(vCode));
+    s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
+  }
+
+  inline void Decay::setStable(Code const vCode) {
+    CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"stable\". ", vCode);
+
+    const int s_id = abs(sibyll::convertToSibyllRaw(vCode));
+    s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+  }
+
+  inline void Decay::setAllStable() {
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
+  }
+
+  inline void Decay::setAllUnstable() {
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+  }
+
+  inline void Decay::printDecayConfig([[maybe_unused]] const Code vCode) {
+    [[maybe_unused]] const int sibCode = corsika::sibyll::convertToSibyllRaw(vCode);
+    [[maybe_unused]] const int absSibCode = abs(sibCode);
+    CORSIKA_LOGGER_DEBUG(logger_, "decay configuration: {} is \"{}\"", vCode,
+                         (s_csydec_.idb[absSibCode - 1] <= 0) ? "stable" : "unstable");
+  }
+  inline void Decay::printDecayConfig() {
+    CORSIKA_LOGGER_DEBUG(logger_, "decay configuration:");
+    if (handleAllDecays_) {
+      CORSIKA_LOGGER_DEBUG(
+          logger_, "     all particles known to Sibyll are handled by Sibyll::Decay!");
+
+    } else {
+      for ([[maybe_unused]] auto const pCode : handledDecays_) {
+        CORSIKA_LOGGER_DEBUG(logger_, "      Decay of {} is handled by Sibyll!", pCode);
+      }
+    }
+  }
+
+  template <typename TParticle>
+  inline TimeType Decay::getLifetime(TParticle const& projectile) {
+
+    const Code pid = projectile.getPID();
+    if (Decay::isDecayHandled(pid)) {
+      HEPEnergyType E = projectile.getEnergy();
+      HEPMassType m = projectile.getMass();
+      const double gamma = E / m;
+      const TimeType t0 = get_lifetime(projectile.getPID());
+      auto const lifetime = gamma * t0;
+      [[maybe_unused]] const auto mkin =
+          +(E * E - projectile.getMomentum()
+                        .getSquaredNorm()); // delta_mass(projectile.getMomentum(), E, m);
+      CORSIKA_LOGGER_DEBUG(logger_, "code: {} ", projectile.getPID());
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: t0: {} ", t0);
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: energy: {} GeV ", E / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "momentum: {} GeV ",
+                           projectile.getMomentum().getComponents() / 1_GeV);
+      CORSIKA_LOGGER_DEBUG(logger_, "momentum: shell mass-kin. inv. mass {} {}",
+                           mkin / 1_GeV / 1_GeV, m / 1_GeV * m / 1_GeV);
+      [[maybe_unused]] auto sib_id =
+          corsika::sibyll::convertToSibyllRaw(projectile.getPID());
+      CORSIKA_LOGGER_DEBUG(logger_, "sib mass: {}", get_sibyll_mass2(sib_id));
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: gamma:  {}", gamma);
+      CORSIKA_LOGGER_DEBUG(logger_, "MinStep: tau {} s: ", lifetime / 1_s);
+      return lifetime;
+    }
+    return std::numeric_limits<double>::infinity() * 1_s;
+  } // namespace corsika::sibyll
+
+  template <typename TSecondaryView>
+  inline void Decay::doDecay(TSecondaryView& view) {
+
+    auto projectile = view.getProjectile();
+    const Code pCode = projectile.getPID();
+    // check if sibyll is configured to handle this decay!
+    if (!isDecayHandled(pCode))
+      throw std::runtime_error("STOP! Sibyll not configured to execute this decay!");
+
+    count_++;
+
+    // switch on decay for this particle
+    setUnstable(pCode);
+    printDecayConfig(pCode);
+
+    // call sibyll decay
+    CORSIKA_LOGGER_DEBUG(logger_, "calling Sibyll decay routine..");
+
+    // particle to pass to sibyll decay
+    int inputSibPID = sibyll::convertToSibyllRaw(pCode);
+
+    // particle momentum format: px, py, pz, e, mass. units: GeV
+    double inputMomentum[5];
+    QuantityVector<hepmomentum_d> input_components =
+        projectile.getMomentum().getComponents();
+    for (int idx = 0; idx < 3; ++idx) inputMomentum[idx] = input_components[idx] / 1_GeV;
+
+    inputMomentum[3] = projectile.getEnergy() / 1_GeV;
+    inputMomentum[4] = get_mass(pCode) / 1_GeV;
+    int nFinalParticles;
+
+    double outputMomentum[5 * 10];
+    int outputSibPID[10];
+
+    // run decay routine
+    decpar_sib_(inputSibPID, inputMomentum, nFinalParticles, outputSibPID,
+                &outputMomentum[0]);
+
+    CORSIKA_LOGGER_TRACE(logger_, "number of final state particles: {}", nFinalParticles);
+
+    // reset to stable
+    setStable(pCode);
+
+    CoordinateSystemPtr const& rootCS = get_root_CoordinateSystem();
+
+    // copy particles from sibyll ministack to corsika
+    for (int i = 0; i < nFinalParticles; ++i) {
+      QuantityVector<hepmomentum_d> components = {outputMomentum[10 * 0 + i] * 1_GeV,
+                                                  outputMomentum[10 * 1 + i] * 1_GeV,
+                                                  outputMomentum[10 * 2 + i] * 1_GeV};
+
+      auto const pid = sibyll::convertFromSibyll(
+          static_cast<corsika::sibyll::SibyllCode>(outputSibPID[i]));
+
+      CORSIKA_LOGGER_TRACE(logger_, "final particle i={} id={} p={} GeV", i, pid,
+                           components / 1_GeV);
+
+      auto const p3 = MomentumVector(rootCS, components);
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(p3.getSquaredNorm() + mass * mass) - mass;
+
+      projectile.addSecondary(std::make_tuple(pid, Ekin, p3.normalized()));
+    }
+  }
+
+} // namespace corsika::sibyll

corsika/detail/modules/sibyll/HadronInteractionModel.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/geometry/Point.hpp>
+
+#include <corsika/modules/sibyll/ParticleConversion.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/sibyll/SibStack.hpp>
+#include <corsika/modules/Random.hpp>
+
+#include <sibyll2.3d.hpp>
+
+#include <tuple>
+
+namespace corsika::sibyll {
+
+  inline void HadronInteractionModel::setVerbose(bool const flag) {
+    sibyll_listing_ = flag;
+  }
+
+  inline HadronInteractionModel::HadronInteractionModel()
+      : sibyll_listing_(false)
+      , internal_decays_(false)
+      , stable_particles_{} {
+    // initialize Sibyll
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+    static bool initialized = false;
+    if (!initialized) {
+      sibyll_ini_();
+      initialized = true;
+    }
+  }
+
+  inline HadronInteractionModel::HadronInteractionModel(std::set<Code> vlist)
+      : sibyll_listing_(false)
+      , internal_decays_(true)
+      , stable_particles_(vlist) {
+    // initialize Sibyll
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+    static bool initialized = false;
+    if (!initialized) {
+      sibyll_ini_();
+      initialized = true;
+    }
+  }
+
+  inline void HadronInteractionModel::setAllParticlesUnstable() {
+    // activate all decays in SIBYLL
+    CORSIKA_LOGGER_DEBUG(logger_, "setting all particles \"unstable\".");
+    for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+  }
+
+  inline void HadronInteractionModel::setParticleListStable(std::set<Code> vList) {
+    // de-activate specific decays in SIBYLL
+    for (auto const p : vList) {
+      CORSIKA_LOGGER_DEBUG(logger_, "setting {} as \"stable\". ", p);
+      auto const sib_code = sibyll::convertToSibyll(p);
+      if (sib_code != corsika::sibyll::SibyllCode::Unknown) {
+        int const s_id = abs(static_cast<int>(sib_code));
+        s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
+      } else {
+        CORSIKA_LOGGER_WARN(logger_,
+                            "particle {} not known to SIBYLL! Cannot set stable!", p);
+      }
+    }
+  }
+
+  inline HadronInteractionModel::~HadronInteractionModel() {
+    CORSIKA_LOGGER_DEBUG(logger_, "Sibyll::Model n={}, Nnuc={}", count_, nucCount_);
+  }
+
+  inline bool constexpr HadronInteractionModel::isValid(
+      Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn) const {
+    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
+
+    if (is_nucleus(targetId)) {
+      size_t const targA = get_nucleus_A(targetId);
+      if (targA != 1 && (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) {
+        return false;
+      }
+    } else if (targetId != Code::Proton && targetId != Code::Neutron &&
+               targetId != Code::Hydrogen) {
+      return false;
+    }
+    if (is_nucleus(projectileId) || !corsika::sibyll::canInteract(projectileId)) {
+      return false;
+    }
+    return true;
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  HadronInteractionModel::getCrossSectionInelEla(Code const projectileId,
+                                                 Code const targetId,
+                                                 FourMomentum const& projectileP4,
+                                                 FourMomentum const& targetP4) const {
+
+    int targetSibCode = 1; // nucleon or particle count
+    if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
+
+    if (!isValid(projectileId, targetId, sqrtSnn)) {
+      return {CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+
+    double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
+    int const iBeam = corsika::sibyll::getSibyllXSCode(
+        projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like,
+                       // 3:kaon-like)
+    if (!iBeam) {
+      CORSIKA_LOGGER_TRACE(
+          logger_, "The cross section for projectile {} and target {} is zero, E={}",
+          projectileId, targetId, sqrtSnn);
+      return {CrossSectionType::zero(), CrossSectionType::zero()};
+    }
+
+    double const dEcm = sqrtSnn / 1_GeV;
+    // single nucleon target (p,n, hydrogen) or 4<=A<=18
+    double sigProd = 0;
+    double sigEla = 0;
+    if (targetId == Code::Proton || targetId == Code::Hydrogen ||
+        targetId == Code::Neutron) {
+      // single nucleon target
+      sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+    } else {
+      // nuclear target
+      int const iTarget = get_nucleus_A(targetId);
+      sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+    }
+    return {sigProd * 1_mb, sigEla * 1_mb};
+  }
+
+  /**
+   * In this function SIBYLL is called to produce one event. The
+   * event is copied (and boosted) into the shower lab frame.
+   */
+
+  template <typename TSecondaryView>
+  inline void HadronInteractionModel::doInteraction(TSecondaryView& secondaries,
+                                                    Code const projectileId,
+                                                    Code const targetId,
+                                                    FourMomentum const& projectileP4,
+                                                    FourMomentum const& targetP4) {
+
+    int targetSibCode = 1; // nucleon or particle count
+    if (is_nucleus(targetId)) { targetSibCode = get_nucleus_A(targetId); }
+    CORSIKA_LOGGER_DEBUG(logger_, "sibyll code: {} (nucleon/particle count)",
+                         targetSibCode);
+
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtSnn = (projectileP4 + targetP4 / targetSibCode).getNorm();
+    COMBoost const boost(projectileP4, targetP4 / targetSibCode);
+
+    if (!isValid(projectileId, targetId, sqrtSnn)) {
+      throw std::runtime_error("Invalid target/projectile/energy combination");
+    }
+
+    CORSIKA_LOGGER_DEBUG(logger_, "pId={} tId={} sqrtSnn={}GeV", projectileId, targetId,
+                         sqrtSnn);
+
+    // beam id for sibyll
+    int const projectileSibyllCode = corsika::sibyll::convertToSibyllRaw(projectileId);
+
+    count_++;
+    // Sibyll does not know about units..
+    double const sqs = sqrtSnn / 1_GeV;
+    if (internal_decays_) {
+      setAllParticlesUnstable();
+      setParticleListStable(stable_particles_);
+    }
+    // running sibyll, filling stack
+    sibyll_(projectileSibyllCode, targetSibCode, sqs);
+    if (internal_decays_) decsib_();
+
+    if (sibyll_listing_) {
+      // print final state
+      int print_unit = 6;
+      sib_list_(print_unit);
+      nucCount_ += get_nwounded() - 1;
+    }
+
+    // ------ output and particle readout -----
+    auto const& csPrime = boost.getRotatedCS();
+
+    // add particles from sibyll to stack
+
+    // link to sibyll stack
+    SibStack ss;
+
+    auto const& originalCS = boost.getOriginalCS();
+    MomentumVector Plab_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
+    for (auto const& psib : ss) {
+      // abort on particles that have decayed in Sibyll. Should not happen!
+      if (psib.hasDecayed()) { // LCOV_EXCL_START
+        if (internal_decays_) {
+          continue;
+        } else {
+          throw std::runtime_error(
+              "internal decayse switched off! found particle that decayed in SIBYLL!");
+        } // LCOV_EXCL_STOP
+      }
+      // transform 4-momentum to lab. frame
+      // note that the momentum needs to be rotated back
+      auto const tmp = psib.getMomentum().getComponents();
+      auto const pCoM = MomentumVector(csPrime, tmp);
+      HEPEnergyType const eCoM = psib.getEnergy();
+      auto const P4lab = boost.fromCoM(FourVector{eCoM, pCoM});
+      auto const p3lab = P4lab.getSpaceLikeComponents();
+
+      Code const pid = corsika::sibyll::convertFromSibyll(psib.getPID());
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
+      // add to corsika stack
+      auto pnew =
+          secondaries.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
+
+      Plab_final += pnew.getMomentum();
+      Elab_final += pnew.getEnergy();
+      Ecm_final += psib.getEnergy();
+    }
+    { // just output
+      [[maybe_unused]] HEPEnergyType const Elab_initial =
+          static_pow<2>(sqrtSnn) / (2 * constants::nucleonMass);
+      CORSIKA_LOGGER_DEBUG(logger_,
+                           "conservation (all GeV): "
+                           "sqrtSnn={}, sqrtSnn_final={}, "
+                           "Elab_initial={}, Elab_final={}, "
+                           "diff(%)={}, "
+                           "E in nucleons={}, "
+                           "Plab_final={} ",
+                           sqrtSnn / 1_GeV, Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV,
+                           Elab_initial, Elab_final / 1_GeV,
+                           (Elab_final - Elab_initial) / Elab_initial * 100,
+                           constants::nucleonMass * get_nwounded() / 1_GeV,
+                           (Plab_final / 1_GeV).getComponents());
+    }
+  }
+} // namespace corsika::sibyll

corsika/detail/modules/sibyll/InteractionModel.inl

+/*
+ * (c) Copyright 2022 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/FourVector.hpp>
+#include <corsika/modules/sibyll/HadronInteractionModel.hpp>
+#include <corsika/modules/sibyll/NuclearInteractionModel.hpp>
+
+namespace corsika::sibyll {
+
+  inline InteractionModel::InteractionModel(std::set<Code> const& nuccomp,
+                                            std::set<Code> const& list)
+      : hadronSibyll_{list}
+      , nuclearSibyll_{hadronSibyll_, nuccomp} {}
+
+  inline HadronInteractionModel& InteractionModel::getHadronInteractionModel() {
+    return hadronSibyll_;
+  }
+
+  inline HadronInteractionModel const& InteractionModel::getHadronInteractionModel()
+      const {
+    return hadronSibyll_;
+  }
+
+  inline typename InteractionModel::nuclear_model_type&
+  InteractionModel::getNuclearInteractionModel() {
+    return nuclearSibyll_;
+  }
+
+  inline typename InteractionModel::nuclear_model_type const&
+  InteractionModel::getNuclearInteractionModel() const {
+    return nuclearSibyll_;
+  }
+
+  inline CrossSectionType InteractionModel::getCrossSection(
+      Code projCode, Code targetCode, FourMomentum const& proj4mom,
+      FourMomentum const& target4mom) const {
+    if (is_nucleus(projCode))
+      return getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                          target4mom);
+    else
+      return getHadronInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                         target4mom);
+  }
+
+  inline std::tuple<CrossSectionType, CrossSectionType>
+  InteractionModel::getCrossSectionInelEla(Code projCode, Code targetCode,
+                                           FourMomentum const& proj4mom,
+                                           FourMomentum const& target4mom) const {
+    if (is_nucleus(projCode))
+      return {getNuclearInteractionModel().getCrossSection(projCode, targetCode, proj4mom,
+                                                           target4mom),
+              CrossSectionType::zero()};
+    else
+      return getHadronInteractionModel().getCrossSectionInelEla(projCode, targetCode,
+                                                                proj4mom, target4mom);
+  }
+
+  template <typename TSecondaries>
+  inline void InteractionModel::doInteraction(TSecondaries& view, Code projCode,
+                                              Code targetCode,
+                                              FourMomentum const& proj4mom,
+                                              FourMomentum const& target4mom) {
+    if (is_nucleus(projCode))
+      return getNuclearInteractionModel().doInteraction(view, projCode, targetCode,
+                                                        proj4mom, target4mom);
+    else
+      return getHadronInteractionModel().doInteraction(view, projCode, targetCode,
+                                                       proj4mom, target4mom);
+  }
+} // namespace corsika::sibyll

corsika/detail/modules/sibyll/NuclearInteractionModel.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/media/Environment.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <nuclib.hpp>
+
+namespace corsika::sibyll {
+
+  template <typename TNucleonModel>
+  inline NuclearInteractionModel<TNucleonModel>::NuclearInteractionModel(
+      TNucleonModel& hadint, std::set<Code> const& nuccomp)
+      : hadronicInteraction_(hadint) {
+
+    // initialize nuclib
+    // TODO: make sure this does not overlap with sibyll
+    corsika::connect_random_stream("sibyll", ::sibyll::set_rng_function);
+    nuc_nuc_ini_();
+
+    // initialize cross sections
+    initializeNuclearCrossSections(nuccomp);
+  }
+
+  template <typename TNucleonModel>
+  inline NuclearInteractionModel<TNucleonModel>::~NuclearInteractionModel() {
+    CORSIKA_LOGGER_DEBUG(logger_, "nuclear interactions handled by Sibyll n={}", count_);
+  }
+
+  template <typename TNucleonModel>
+  inline bool constexpr NuclearInteractionModel<TNucleonModel>::isValid(
+      Code const projectileId, Code const targetId, HEPEnergyType const sqrtSnn) const {
+
+    // also depends on underlying model, for Proton/Neutron projectile
+    if (!hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn)) { return false; }
+
+    // projectile limits:
+    if (!is_nucleus(projectileId)) { return false; }
+    unsigned int projectileA = get_nucleus_A(projectileId);
+    if (projectileA > getMaxNucleusAProjectile() || projectileA < 2) { return false; }
+    return true;
+  } // namespace corsika::sibyll
+
+  template <typename TNucleonModel>
+  inline void NuclearInteractionModel<TNucleonModel>::printCrossSectionTable(
+      Code const pCode) const {
+    if (!hadronicInteraction_.isValid(Code::Proton, pCode, 100_GeV)) { // LCOV_EXCL_START
+      CORSIKA_LOGGER_WARN(logger_, "Invalid target type {} for hadron interaction model.",
+                          pCode);
+      return;
+    } // LCOV_EXCL_STOP
+
+    int const k = targetComponentsIndex_.at(pCode);
+    Code const pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
+                            Code::Neon,   Code::Argon,    Code::Iron};
+
+    std::ostringstream table;
+    table << "Nuclear CrossSectionTable pCode=" << pCode << " :\n en/A ";
+    for (auto& j : pNuclei) table << std::setw(9) << j;
+    table << "\n";
+
+    // loop over energy bins
+    for (unsigned int i = 0; i < getNEnergyBins(); ++i) {
+      table << " " << i << "  ";
+
+      for (auto& n : pNuclei) {
+        auto const j = get_nucleus_A(n);
+        table << " " << std::setprecision(5) << std::setw(8)
+              << cnucsignuc_.sigma[j - 1][k][i];
+      }
+      table << "\n";
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, table.str());
+  }
+
+  template <typename TNucleonModel>
+  inline void NuclearInteractionModel<TNucleonModel>::initializeNuclearCrossSections(
+      std::set<Code> const& allElementsInUniverse) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "initializing nuclear cross sections...");
+
+    // loop over target components, at most 4!!
+    int k = -1;
+    for (Code const ptarg : allElementsInUniverse) {
+      ++k;
+      CORSIKA_LOGGER_DEBUG(logger_, "init target component: {} A={}", ptarg,
+                           get_nucleus_A(ptarg));
+      int const ib = get_nucleus_A(ptarg);
+      if (!hadronicInteraction_.isValid(Code::Proton, ptarg, 100_GeV)) {
+        CORSIKA_LOGGER_WARN(
+            logger_, "Invalid target type {} for hadron interaction model.", ptarg);
+        continue;
+      }
+      targetComponentsIndex_.insert(std::pair<Code, int>(ptarg, k));
+      // loop over energies, fNEnBins log. energy bins
+      for (size_t i = 0; i < getNEnergyBins(); ++i) {
+        // hard coded energy grid, has to be aligned to definition in signuc2!!, no
+        // comment..
+        HEPEnergyType const Ecm = pow(10., 1. + 1. * i) * 1_GeV;
+        // head-on pp collision:
+        HEPEnergyType const EcmHalve = Ecm / 2;
+        HEPMomentumType const pcm =
+            sqrt(EcmHalve * EcmHalve - Proton::mass * Proton::mass);
+        CoordinateSystemPtr cs = get_root_CoordinateSystem();
+        FourMomentum projectileP4(EcmHalve, {cs, pcm, 0_eV, 0_eV});
+        FourMomentum targetP4(EcmHalve, {cs, -pcm, 0_eV, 0_eV});
+        // get p-p cross sections
+        if (!hadronicInteraction_.isValid(Code::Proton, Code::Proton, Ecm)) {
+          throw std::runtime_error("invalid (projectile,target,ecm) combination");
+        }
+        auto const [siginel, sigela] = hadronicInteraction_.getCrossSectionInelEla(
+            Code::Proton, Code::Proton, projectileP4, targetP4);
+        double const dsig = siginel / 1_mb;
+        double const dsigela = sigela / 1_mb;
+        // loop over projectiles, mass numbers from 2 to fMaxNucleusAProjectile
+        CORSIKA_LOGGER_TRACE(logger_, "Ecm={} siginel={} sigela={}", Ecm / 1_GeV, dsig,
+                             dsigela);
+        for (size_t j = 1; j < gMaxNucleusAProjectile_; ++j) {
+          const int jj = j + 1;
+          double sig_out, dsig_out, sigqe_out, dsigqe_out;
+          sigma_mc_(jj, ib, dsig, dsigela, gNSample_, sig_out, dsig_out, sigqe_out,
+                    dsigqe_out);
+          // write to table
+          cnucsignuc_.sigma[j][k][i] = sig_out;
+          cnucsignuc_.sigqe[j][k][i] = sigqe_out;
+          CORSIKA_LOGGER_TRACE(logger_, "nuc A={} sig={} qe={}", j, sig_out, sigqe_out);
+        }
+      }
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, "cross sections for {} components initialized!",
+                         targetComponentsIndex_.size());
+    for (auto& ptarg : allElementsInUniverse) { printCrossSectionTable(ptarg); }
+  }
+
+  template <typename TNucleonModel>
+  inline CrossSectionType NuclearInteractionModel<TNucleonModel>::readCrossSectionTable(
+      int const ia, Code const pTarget, HEPEnergyType const elabnuc) const {
+    CORSIKA_LOGGER_DEBUG(logger_, "ia={}, target={}, ElabNuc={} GeV", ia, pTarget,
+                         elabnuc / 1_GeV);
+    int const ib = targetComponentsIndex_.at(pTarget) + 1; // table index in fortran
+    auto const ECoMNuc = sqrt(2. * constants::nucleonMass * elabnuc);
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtSnn= {} GeV", ECoMNuc / 1_GeV);
+    if (ECoMNuc < getMinEnergyPerNucleonCoM() || ECoMNuc > getMaxEnergyPerNucleonCoM()) {
+      CORSIKA_LOGGER_WARN(
+          logger_,
+          "nucleon-nucleon energy outside range! sqrtSnn={}GeV (limits: {} .. {} GeV)",
+          ECoMNuc / 1_GeV, getMinEnergyPerNucleonCoM() / 1_GeV,
+          getMaxEnergyPerNucleonCoM() / 1_GeV);
+      // throw std::runtime_error("energy outside tabulated range!");
+    }
+    double const e0 = elabnuc / 1_GeV;
+    double sig;
+    CORSIKA_LOGGER_DEBUG(logger_, "ReadCrossSectionTable: {} {} {}", ia, ib, e0);
+    signuc2_(ia, ib, e0, sig);
+    CORSIKA_LOGGER_DEBUG(logger_, "ReadCrossSectionTable: sig={}", sig);
+    return sig * 1_mb;
+  }
+
+  template <typename TNucleonModel>
+  CrossSectionType inline NuclearInteractionModel<TNucleonModel>::getCrossSection(
+      Code const projectileId, Code const targetId, FourMomentum const& projectileP4,
+      FourMomentum const& targetP4) const {
+    CORSIKA_LOGGER_DEBUG(logger_, "projectile: E={}, p3={} \n target: E={}, p3={}",
+                         projectileP4.getTimeLikeComponent() / 1_GeV,
+                         projectileP4.getSpaceLikeComponents() / 1_GeV,
+                         targetP4.getTimeLikeComponent() / 1_GeV,
+                         targetP4.getSpaceLikeComponents() / 1_GeV);
+    // check if projectile and target are nuclei!
+    if (!is_nucleus(projectileId) || !is_nucleus(targetId)) {
+      return CrossSectionType::zero();
+    }
+
+    // calculate sqrt(Snn) (only works if projectile and target are nuclei)
+    HEPEnergyType const sqrtSnn =
+        (projectileP4 / get_nucleus_A(projectileId) + targetP4 / get_nucleus_A(targetId))
+            .getNorm();
+
+    CORSIKA_LOG_DEBUG("proj={}, targ={}, sqrtSNN={}GeV", projectileId, targetId,
+                      sqrtSnn / 1_GeV);
+    if (!isValid(projectileId, targetId, sqrtSnn)) { return CrossSectionType::zero(); }
+
+    // lab-frame energy per projectile nucleon as required by signuc2()
+    HEPEnergyType const LabEnergyPerNuc = calculate_lab_energy(
+        static_pow<2>(sqrtSnn), get_mass(projectileId) / get_nucleus_A(projectileId),
+        get_mass(targetId) / get_nucleus_A(targetId));
+    auto const sigProd =
+        readCrossSectionTable(get_nucleus_A(projectileId), targetId, LabEnergyPerNuc);
+    CORSIKA_LOGGER_DEBUG(logger_, "cross section (mb): sqrtSnn={} sig={}",
+                         sqrtSnn / 1_GeV, sigProd / 1_mb);
+    return sigProd;
+  }
+
+  template <typename TNucleonModel>
+  template <typename TSecondaryView>
+  inline void NuclearInteractionModel<TNucleonModel>::doInteraction(
+      TSecondaryView& view, Code const projectileId, Code const targetId,
+      FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+
+    // model is only designed for projectile nuclei. Collisions are broken down into
+    // "nucleon-target" collisions.
+    if (!is_nucleus(projectileId)) {
+      throw std::runtime_error("Can only handle nuclear projectiles.");
+    }
+    size_t const projectileA = get_nucleus_A(projectileId);
+
+    // this is center-of-mass for projectile_nucleon - target
+    FourMomentum const nucleonP4 = projectileP4 / projectileA;
+    HEPEnergyType const sqrtSnucleon = (nucleonP4 + targetP4).getNorm();
+    if (!isValid(projectileId, targetId, sqrtSnucleon)) {
+      throw std::runtime_error("Invalid projectile/target/energy combination.");
+    }
+    // projectile is always nucleus!
+    // Elab corresponding to sqrtSnucleon -> fixed target projectile
+    COMBoost const boost(nucleonP4, targetP4);
+
+    CORSIKA_LOGGER_DEBUG(logger_, "pId={} tId={} sqrtSnucleon={}GeV Aproj={}",
+                         projectileId, targetId, sqrtSnucleon / 1_GeV, projectileA);
+    count_++;
+
+    // lab. momentum per projectile nucleon
+    HEPMomentumType const pNucleonLab = nucleonP4.getSpaceLikeComponents().getNorm();
+    // nucleon momentum in direction of CM motion (lab system)
+    MomentumVector const p3NucleonLab(boost.getRotatedCS(), {0_GeV, 0_GeV, pNucleonLab});
+
+    /*
+      FOR NOW: allow nuclei with A<18 or protons/nucleon only.
+      when medium composition becomes more complex, approximations will have to be
+      allowed air in atmosphere also contains some Argon.
+    */
+    int kATarget = -1;
+    size_t targetA = 1;
+    if (is_nucleus(targetId)) {
+      kATarget = get_nucleus_A(targetId);
+      targetA = kATarget;
+    } else if (targetId == Code::Proton || targetId == Code::Neutron ||
+               targetId == Code::Hydrogen) {
+      kATarget = 1;
+    }
+    CORSIKA_LOGGER_DEBUG(logger_, "nuclib target code: {}", kATarget);
+
+    // end of target sampling
+
+    // superposition
+    CORSIKA_LOGGER_DEBUG(logger_, "sampling nuc. multiple interaction structure.. ");
+    // get nucleon-nucleon cross section
+    // (needed to determine number of nucleon-nucleon scatterings)
+    auto const protonId = Code::Proton;
+    auto const [prodCrossSection, elaCrossSection] =
+        hadronicInteraction_.getCrossSectionInelEla(
+            protonId, protonId, nucleonP4,
+            targetP4 / targetA); // todo check, wrong RU
+    double const sigProd = prodCrossSection / 1_mb;
+    double const sigEla = elaCrossSection / 1_mb;
+    // sample number of interactions (only input variables, output in common cnucms)
+    // nuclear multiple scattering according to glauber (r.i.p.)
+    int_nuc_(kATarget, projectileA, sigProd, sigEla);
+
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "number of nucleons in target           : {}\n"
+                         "number of wounded nucleons in target   : {}\n"
+                         "number of nucleons in projectile       : {}\n"
+                         "number of wounded nucleons in project. : {}\n"
+                         "number of inel. nuc.-nuc. interactions : {}\n"
+                         "number of elastic nucleons in target   : {}\n"
+                         "number of elastic nucleons in project. : {}\n"
+                         "impact parameter: {}",
+                         kATarget, cnucms_.na, projectileA, cnucms_.nb, cnucms_.ni,
+                         cnucms_.nael, cnucms_.nbel, cnucms_.b);
+
+    // calculate fragmentation
+    CORSIKA_LOGGER_DEBUG(logger_, "calculating nuclear fragments..");
+    // number of interactions
+    // include elastic
+    int const nElasticNucleons = cnucms_.nbel;
+    int const nInelNucleons = cnucms_.nb;
+    int const nIntProj = nInelNucleons + nElasticNucleons;
+    double const impactPar = cnucms_.b; // only needed to avoid passing common var.
+    int nFragments = 0;
+    // number of fragments is limited to 60
+    int AFragments[60];
+    // call fragmentation routine
+    // input: target A, projectile A, number of int. nucleons in projectile, impact
+    // parameter (fm) output: nFragments, AFragments in addition the momenta ar stored
+    // in pf in common fragments, neglected
+    fragm_(kATarget, projectileA, nIntProj, impactPar, nFragments, AFragments);
+
+    // this should not occur but well :)  (LCOV_EXCL_START)
+    if (nFragments > (int)getMaxNFragments()) {
+      throw std::runtime_error("Number of nuclear fragments in NUCLIB exceeded!");
+    }
+    // (LCOV_EXCL_STOP)
+
+    CORSIKA_LOGGER_DEBUG(logger_, "number of fragments: {}", nFragments);
+    CORSIKA_LOGGER_DEBUG(logger_, "adding nuclear fragments to particle stack..");
+    // put nuclear fragments on corsika stack
+    for (int j = 0; j < nFragments; ++j) {
+      CORSIKA_LOGGER_DEBUG(logger_, "fragment {}: A={} px={} py={} pz={}", j,
+                           AFragments[j], fragments_.ppp[j][0], fragments_.ppp[j][1],
+                           fragments_.ppp[j][2]);
+      auto const nuclA = AFragments[j];
+      // get Z from stability line
+      auto const nuclZ = int(nuclA / 2.15 + 0.7);
+
+      // TODO: do we need to catch single nucleons??
+      Code const specCode = (nuclA == 1 ?
+                                        // TODO: sample neutron or proton
+                                 Code::Proton
+                                        : get_nucleus_code(nuclA, nuclZ));
+      HEPMassType const mass = get_mass(specCode);
+
+      CORSIKA_LOGGER_DEBUG(logger_, "adding fragment: {}", get_name(specCode));
+      CORSIKA_LOGGER_DEBUG(logger_, "A,Z: {}, {}", nuclA, nuclZ);
+      CORSIKA_LOGGER_DEBUG(logger_, "mass: {} GeV", mass / 1_GeV);
+
+      // CORSIKA 7 way
+      // spectators inherit momentum from original projectile
+      auto const p3lab = p3NucleonLab * nuclA;
+
+      HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
+      CORSIKA_LOGGER_DEBUG(logger_, "fragment momentum {}",
+                           p3lab.getComponents() / 1_GeV);
+      view.addSecondary(std::make_tuple(specCode, Ekin, p3lab.normalized()));
+    }
+
+    // add elastic nucleons to corsika stack
+    // TODO: the elastic interaction could be external like the inelastic interaction,
+    // e.g. use existing ElasticModel
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "adding elastically scattered nucleons to particle stack..");
+    for (int j = 0; j < nElasticNucleons; ++j) {
+      // TODO: sample proton or neutron
+      Code const elaNucCode = Code::Proton;
+
+      // CORSIKA 7 way
+      // elastic nucleons inherit momentum from original projectile
+      // neglecting momentum transfer in interaction
+      auto const p3lab = p3NucleonLab;
+
+      HEPEnergyType const mass = get_mass(elaNucCode);
+      HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+
+      view.addSecondary(std::make_tuple(elaNucCode, Ekin, p3lab.normalized()));
+    }
+
+    // add inelastic interactions
+    CORSIKA_LOGGER_DEBUG(logger_, "calculate inelastic nucleon-nucleon interactions..");
+    for (int j = 0; j < nInelNucleons; ++j) {
+      // TODO: sample neutron or proton
+      auto const pCode = Code::Proton;
+      HEPEnergyType const mass = get_mass(pCode);
+      HEPEnergyType const Ekin = sqrt(p3NucleonLab.getSquaredNorm() + mass * mass) - mass;
+
+      // temporarily add to stack, will be removed after interaction in DoInteraction
+      CORSIKA_LOGGER_DEBUG(logger_, "inelastic interaction no. {}", j);
+      typename TSecondaryView::inner_stack_value_type nucleonStack;
+      Point const pDummy(boost.getOriginalCS(), {0_m, 0_m, 0_m});
+      TimeType const tDummy = 0_ns;
+      auto inelasticNucleon = nucleonStack.addParticle(
+          std::make_tuple(pCode, Ekin, p3NucleonLab.normalized(), pDummy, tDummy));
+      inelasticNucleon.setNode(view.getProjectile().getNode());
+
+      // create inelastic interaction for each nucleon
+      CORSIKA_LOGGER_TRACE(logger_, "calling HadronicInteraction...");
+      // create new StackView for each of the nucleons
+      TSecondaryView nucleon_secondaries(inelasticNucleon);
+      // all inner hadronic event generator
+      hadronicInteraction_.doInteraction(nucleon_secondaries, pCode, targetId, nucleonP4,
+                                         targetP4);
+      for (const auto& pSec : nucleon_secondaries) {
+
+        auto const p3lab = pSec.getMomentum();
+        Code const pid = pSec.getPID();
+        HEPEnergyType const mass = get_mass(pid);
+        HEPEnergyType const Ekin = sqrt(p3lab.getSquaredNorm() + mass * mass) - mass;
+        view.addSecondary(std::make_tuple(pid, Ekin, p3lab.normalized()));
+      }
+    }
+  }
+} // namespace corsika::sibyll

corsika/detail/modules/sibyll/ParticleConversion.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/ParticleProperties.hpp>
+
+#include <sibyll2.3d.hpp>
+
+namespace corsika::sibyll {
+
+  inline HEPMassType getSibyllMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
+    auto sCode = convertToSibyllRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getSibyllMass: unknown particle!");
+    else
+      return sqrt(get_sibyll_mass2(sCode)) * 1_GeV;
+  }
+} // namespace corsika::sibyll
\ No newline at end of file

corsika/detail/modules/sophia/InteractionModel.inl

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/modules/sophia/ParticleConversion.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <corsika/modules/Random.hpp>
+#include <corsika/modules/sophia/SophiaStack.hpp>
+#include <corsika/framework/core/EnergyMomentumOperations.hpp>
+
+#include <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline void InteractionModel::setVerbose(bool const flag) { sophia_listing_ = flag; }
+
+  inline InteractionModel::InteractionModel()
+      : sophia_listing_(false) {
+    corsika::connect_random_stream(RNG_, ::sophia::set_rng_function);
+    // set all particles stable in SOPHIA
+    for (int i = 0; i < 49; ++i) so_csydec_.idb[i] = -abs(so_csydec_.idb[i]);
+  }
+
+  inline InteractionModel::~InteractionModel() {
+    CORSIKA_LOG_DEBUG("Sophia::Model n={}", count_);
+  }
+
+  inline bool constexpr InteractionModel::isValid(Code const projectileId,
+                                                  Code const targetId,
+                                                  HEPEnergyType const sqrtSnn) const {
+    if ((minEnergyCoM_ > sqrtSnn) || (sqrtSnn > maxEnergyCoM_)) { return false; }
+
+    if (!(targetId == Code::Proton || targetId == Code::Neutron ||
+          targetId == Code::Hydrogen))
+      return false;
+
+    if (projectileId != Code::Photon) return false;
+
+    return true;
+  }
+
+  template <typename TSecondaryView>
+  inline void InteractionModel::doInteraction(TSecondaryView& secondaries,
+                                              Code const projectileId,
+                                              Code const targetId,
+                                              FourMomentum const& projectileP4,
+                                              FourMomentum const& targetP4) {
+
+    CORSIKA_LOGGER_DEBUG(logger_, "projectile: Id={}, E={} GeV, p3={} GeV", projectileId,
+                         projectileP4.getTimeLikeComponent() / 1_GeV,
+                         projectileP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+    CORSIKA_LOGGER_DEBUG(logger_, "target: Id={}, E={} GeV, p3={} GeV", targetId,
+                         targetP4.getTimeLikeComponent() / 1_GeV,
+                         targetP4.getSpaceLikeComponents().getComponents() / 1_GeV);
+
+    // sqrtS per target nucleon
+    HEPEnergyType const sqrtS = (projectileP4 + targetP4).getNorm();
+
+    CORSIKA_LOGGER_DEBUG(logger_, "sqrtS={}GeV", sqrtS / 1_GeV);
+
+    // accepts only photon-nucleon interactions
+    if (!isValid(projectileId, targetId, sqrtS)) {
+      CORSIKA_LOGGER_ERROR(logger_,
+                           "Invalid target/projectile/energy combination: {},{},{} GeV",
+                           projectileId, targetId, sqrtS / 1_GeV);
+      throw std::runtime_error("SOPHIA: Invalid target/projectile/energy combination");
+    }
+
+    COMBoost const boost(projectileP4, targetP4);
+
+    int nucleonSophiaCode = convertToSophiaRaw(targetId); // either proton or neutron
+    // initialize resonance spectrum
+    initial_(nucleonSophiaCode);
+    // Sophia does sqrt(1 - mass_sophia / mass_c8), so we need to make sure that E0 >=
+    // m_sophia
+    double Enucleon = std::max(corsika::sophia::getSophiaMass(targetId) / 1_GeV,
+                               targetP4.getTimeLikeComponent() / 1_GeV);
+    double Ephoton = projectileP4.getTimeLikeComponent() / 1_GeV;
+    double theta = 0.0; // set nucleon at rest in collision
+    int Imode = -1;     // overwritten inside SOPHIA
+    CORSIKA_LOGGER_DEBUG(logger_,
+                         "calling SOPHIA eventgen with L0={}, E0={}, eps={},theta={}",
+                         nucleonSophiaCode, Enucleon, Ephoton, theta);
+    count_++;
+    // call sophia
+    eventgen_(nucleonSophiaCode, Enucleon, Ephoton, theta, Imode);
+
+    if (sophia_listing_) {
+      int arg = 3;
+      print_event_(arg);
+    }
+
+    auto const& originalCS = boost.getOriginalCS();
+    // SOPHIA has photon along -z  and nucleon along +z (GZK calc..)
+    COMBoost const boostInternal(targetP4, projectileP4);
+    auto const& csPrime = boost.getRotatedCS();
+    CoordinateSystemPtr csPrimePrime =
+        make_rotation(csPrime, QuantityVector<length_d>{1_m, 0_m, 0_m}, M_PI);
+
+    SophiaStack ss;
+
+    MomentumVector P_final(originalCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+    HEPEnergyType E_final = 0_GeV;
+    for (auto& psop : ss) {
+      // abort on particles that have decayed in SOPHIA. Should not happen!
+      if (psop.hasDecayed()) { // LCOV_EXCL_START
+        throw std::runtime_error("found particle that decayed in SOPHIA!");
+      } // LCOV_EXCL_STOP
+
+      auto momentumSophia = psop.getMomentum(csPrimePrime);
+      momentumSophia.rebase(csPrime);
+      auto const energySophia = psop.getEnergy();
+      auto const P4com = boostInternal.toCoM(FourVector{energySophia, momentumSophia});
+      auto const P4lab = boost.fromCoM(P4com);
+      SophiaCode const pidSophia = psop.getPID();
+      Code const pid = convertFromSophia(pidSophia);
+      auto momentum = P4lab.getSpaceLikeComponents();
+      momentum.rebase(originalCS);
+      HEPEnergyType const Ekin =
+          calculate_kinetic_energy(momentum.getNorm(), get_mass(pid));
+
+      CORSIKA_LOGGER_TRACE(logger_, "SOPHIA: pid={}, p={} GeV", pidSophia,
+                           momentumSophia.getComponents() / 1_GeV);
+
+      CORSIKA_LOGGER_TRACE(logger_, "CORSIKA: pid={}, p={} GeV", pid,
+                           momentum.getComponents() / 1_GeV);
+
+      auto pnew =
+          secondaries.addSecondary(std::make_tuple(pid, Ekin, momentum.normalized()));
+
+      P_final += pnew.getMomentum();
+      E_final += pnew.getEnergy();
+    }
+    CORSIKA_LOGGER_TRACE(logger_, "Efinal={} GeV,Pfinal={} GeV", E_final / 1_GeV,
+                         P_final.getComponents() / 1_GeV);
+  }
+
+} // namespace corsika::sophia

corsika/detail/modules/sophia/ParticleConversion.inl

+/*
+ * (c) Copyright 2022 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 <sophia.hpp>
+
+namespace corsika::sophia {
+
+  inline HEPMassType getSophiaMass(Code const pCode) {
+    if (is_nucleus(pCode)) throw std::runtime_error("Not defined for nuclei.");
+    auto sCode = convertToSophiaRaw(pCode);
+    if (sCode == 0)
+      throw std::runtime_error("getSophiaMass: unknown particle!");
+    else
+      return sqrt(get_sophia_mass2(sCode)) * 1_GeV;
+  }
+} // namespace corsika::sophia
\ No newline at end of file

corsika/detail/modules/tauola/Decay.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.
+ */
+
+#include <Tauola/Tauola.h>
+
+#include <corsika/framework/utility/COMBoost.hpp>
+
+namespace Tpp = Tauolapp;
+
+namespace corsika::tauola {
+
+  // default constructor
+  Decay::Decay(const Helicity helicity, const bool radiative, const bool spincorr)
+      : helicity_(helicity) {
+
+    // if TAUOLA has not already been initialized,
+    if (!Tpp::Tauola::getIsTauolaIni()) {
+
+      // set the default units explicitly.
+      // we use GeV from momentum and cm for length
+      Tpp::Tauola::setUnits(Tpp::Tauola::GEV, Tpp::Tauola::CM);
+
+      // we set the lifetime as zero so that the tau decays
+      // at the vertex point given by CORSIKA since
+      // CORSIKA handles proapgating the tau until the decay point.
+      Tpp::Tauola::setTauLifetime(0.0);
+
+      // don't decay any intermediaries - add these to the stack
+      // in particular, eta's, K0, Pi-0's.
+      // See Pg. 46 in the TAUOLA manual.
+      Tpp::Tauola::setEtaK0sPi(1, 1, 0);
+
+      // use the new currents (Pg. 40 in the user manual).
+      Tpp::Tauola::setNewCurrents(1);
+
+      // whether to use QCD radiative corrections in calculating
+      // leptonic tau decay. If this is enabled, TAUOLA generates
+      // an order of magnitude more photons than Pythia (for example)
+      Tpp::Tauola::setRadiation(radiative);
+
+      // we have to set the spin correlation state after
+      // we have initialized TAUOLA
+      Tpp::Tauola::spin_correlation.setAll(spincorr);
+
+      // initialize the TAUOLA library
+      Tpp::Tauola::initialize();
+    }
+  }
+
+  Decay::~Decay() {
+    CORSIKA_LOGGER_INFO(logger_, "Number of particles decayed using TAUOLA: {}", count_);
+  }
+
+  bool Decay::isDecayHandled(const Code code) {
+    // if known to tauola, it can decay
+    if (code == Code::WPlus || code == Code::WMinus || code == Code::Z0 ||
+        code == Code::TauPlus || code == Code::TauMinus || code == Code::Eta ||
+        code == Code::K0Short || code == Code::RhoPlus || code == Code::RhoMinus ||
+        code == Code::KStarPlus || code == Code::KStarMinus)
+      return true;
+    else
+      return false;
+  }
+
+  // perform the decay
+  template <typename TSecondaryView>
+  void Decay::doDecay(TSecondaryView& view) {
+    CORSIKA_LOGGER_DEBUG(logger_, "Decaying primary particle using TAUOLA...");
+    count_++;
+
+    // get the coordinate system of the decay
+    auto projectile = view.getProjectile();
+
+    const auto& cs{projectile.getMomentum().getCoordinateSystem()};
+
+    // get the TAUOLA particle from this projectile
+    auto particle{TauolaInterfaceParticle::from_projectile(projectile)};
+
+    CORSIKA_LOGGER_TRACE(
+        logger_, "[decay primary] code: {}, momentum: ({}, {}, {}) GeV, energy: {} GeV",
+        convert_from_PDG(static_cast<PDGCode>(particle->getPdgID())), particle->getPx(),
+        particle->getPy(), particle->getPz(), particle->getE());
+
+    // decay this particle using TAUOLA
+    decay(particle.get());
+
+    std::vector<TauolaInterfaceParticle*> daughters;
+
+    // we iterate over the daughter particles and check for
+    // short-lived secondaries that we have to decay again using TAUOLA.
+    // We only ever have to go one layer down in the hierarchy.
+    for (const auto& daughter : particle->getDaughters()) {
+
+      // check if this is decay can be handled by TAUOLA
+      if (abs(daughter->getPdgID()) == 20213 ||
+          isDecayHandled(convert_from_PDG(static_cast<PDGCode>(daughter->getPdgID())))) {
+        // decay this secondary
+        decay(daughter);
+        // and now add all the daughters of this particle to our list
+        for (const auto& child : daughter->getDaughters()) {
+
+          // check if this is further decayable (K0 and Eta's typically)
+          if (abs(child->getPdgID()) == 20213 ||
+              isDecayHandled(convert_from_PDG(static_cast<PDGCode>(child->getPdgID())))) {
+            // decay the cild
+            decay(child);
+
+            // and add it's daughters to the top-level daughters vector
+            for (const auto& grandchild : child->getDaughters()) {
+              daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(grandchild));
+            }
+          } else {
+            // if we are here, the original daughter is non-decayable so add it to our
+            // struct.
+            daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(child));
+          }
+        }
+
+      } else {
+        // if we get here, `daughter` is long-lived and doesn't require a decay.
+        daughters.push_back(dynamic_cast<TauolaInterfaceParticle*>(daughter));
+      }
+    }
+
+    // we now have a complete list of the daughter particles of this decay
+    // we now iterate over them and add them to CORSIKA stack.
+    for (const auto& daughter : daughters) {
+
+      // get particle's PDG ID
+      const auto pdg{daughter->getPdgID()};
+
+      // get the CORSIKA PID from our PDG ID
+      const auto pid{convert_from_PDG(static_cast<PDGCode>(pdg))};
+
+      // construct the momentum in the ROOT coordinate system
+      const MomentumVector P(cs, {daughter->getPx() * 1_GeV, daughter->getPy() * 1_GeV,
+                                  daughter->getPz() * 1_GeV});
+
+      HEPEnergyType const mass = get_mass(pid);
+      HEPEnergyType const Ekin = sqrt(P.getSquaredNorm() + mass * mass) - mass;
+
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "[decay result] code: {}, momentum: {} GeV, energy: {} GeV",
+                           pid, P.getComponents() / 1_GeV, daughter->getE());
+
+      // and add the particle to the CORSIKA stack
+      projectile.addSecondary(std::make_tuple(pid, Ekin, P.normalized()));
+    }
+  }
+
+  template <typename InterfaceParticle>
+  void Decay::decay(InterfaceParticle* particle) {
+
+    // return if we get a nullptr here
+    // this should never happen (but to be safe!)
+    // LCOV_EXCL_START
+    if (particle == nullptr) {
+      CORSIKA_LOGGER_WARN(logger_, "Got nullptr in decay routine. Skipping decay.");
+      return;
+    }
+    // LCOV_EXCL_STOP
+
+    // get the particle's PDG
+    const auto pdg{static_cast<PDGCode>(particle->getPdgID())};
+
+    // if this is a tau- or tau+, use the helicity in the decay
+    if (pdg == PDGCode::TauMinus || pdg == PDGCode::TauPlus) {
+
+      // get helicity from enum class
+      double helic;
+      switch (helicity_) {
+        case Helicity::LeftHanded:
+          helic = -1.;
+          break;
+        case Helicity::RightHanded:
+          helic = 1.;
+          break;
+        case Helicity::Unpolarized:
+          helic = 0.;
+          break;
+          // NOT testing undefined
+        default: // LCOV_EXCL_START
+          CORSIKA_LOGGER_ERROR(logger_, "Unknown helicity setting: {}", helicity_);
+          throw std::runtime_error("Unknown helicity setting!");
+          // LCOV_EXCL_STOP
+      }
+
+      double Polx{0};
+      double Poly{0};
+      double Polz{-helic};
+
+      // and call the decay routine with the polarization
+      Tpp::Tauola::decayOne(particle, false, Polx, Poly, Polz);
+
+    } else { // this is not a tau, so just do a standard decay
+      Tpp::Tauola::decayOne(particle);
+    }
+  }
+
+  template <typename TParticle>
+  TimeType Decay::getLifetime(TParticle const& particle) {
+    // Can be moved into Cascade.inl, see Issue #271
+    auto const pid = particle.getPID();
+    if ((pid == Code::TauMinus) || (pid == Code::TauPlus)) {
+      HEPEnergyType E = particle.getEnergy();
+      HEPMassType m = particle.getMass();
+      double const gamma = E / m;
+      TimeType const t0 = get_lifetime(pid);
+      auto const lifetime = gamma * t0;
+      CORSIKA_LOGGER_TRACE(logger_,
+                           "[decay time] code: {}, t0: {} ns, momentum: {} GeV, energy: "
+                           "{} GeV, gamma: {}, tau: {} ns",
+                           particle.getPID(), t0 / 1_ns,
+                           particle.getMomentum().getComponents() / 1_GeV, E / 1_GeV,
+                           gamma, lifetime / 1_ns);
+      return lifetime;
+    } else
+      return std::numeric_limits<double>::infinity() * 1_s;
+  }
+
+  void Decay::PrintSummary() const { Tpp::Tauola::summary(); }
+
+} // namespace corsika::tauola

corsika/detail/modules/thinning/EMThinning.inl

+/*
+ * (c) Copyright 2022 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/process/SecondariesProcess.hpp>
+
+namespace corsika {
+
+  EMThinning::EMThinning(HEPEnergyType threshold, double maxWeight,
+                         bool const eraseParticles)
+      : threshold_{threshold}
+      , maxWeight_{maxWeight}
+      , eraseParticles_{eraseParticles} {}
+
+  template <typename TStackView>
+  void EMThinning::doSecondaries(TStackView& view) {
+    if (view.getSize() != 2) return;
+
+    auto projectile = view.getProjectile();
+    if (!is_em(projectile.getPID())) return;
+
+    double const parentWeight = projectile.getWeight();
+    if (parentWeight >= maxWeight_) { return; }
+
+    auto particle1 = view.begin();
+    auto particle2 = ++(view.begin());
+
+    // skip non-EM interactions
+    if (!is_em(particle1.getPID()) || !is_em(particle2.getPID())) { return; }
+
+    /* skip high-energy interactions
+     * We consider only the primary energy. A more sophisticated algorithm might
+     * sort out above-threshold secondaries and apply thinning only on the subset of
+     * those below the threshold.
+     */
+    if (projectile.getEnergy() > threshold_) { return; }
+
+    corsika::HEPEnergyType const Esum = particle1.getEnergy() + particle2.getEnergy();
+
+    double const p1 = particle1.getEnergy() / Esum;
+    double const p2 = particle2.getEnergy() / Esum;
+
+    // weight factors
+    auto const f1 = 1 / p1;
+    auto const f2 = 1 / p2;
+
+    // max. allowed weight factor
+    double const fMax = maxWeight_ / parentWeight;
+
+    if (f1 <= fMax && f2 <= fMax) {
+      /* cannot possibly reach wmax in this vertex; apply
+        Hillas thinning; select one of the two */
+      if (uniform_(rng_) <= p1) { // keep 1st with probability p1
+        particle2.setWeight(0);
+        particle1.setWeight(parentWeight * f1);
+      } else { // keep 2nd
+        particle1.setWeight(0);
+        particle2.setWeight(parentWeight * f2);
+      }
+    } else { // weight limitation kicks in, do statistical thinning
+      double const f1prime = std::min(f1, fMax);
+      double const f2prime = std::min(f2, fMax);
+
+      if (uniform_(rng_) * f1prime <= 1) { // keep 1st
+        particle1.setWeight(parentWeight * f1prime);
+      } else {
+        particle1.setWeight(0);
+      }
+      if (uniform_(rng_) * f2prime <= 1) { // keep 2nd
+        particle2.setWeight(parentWeight * f2prime);
+      } else {
+        particle2.setWeight(0);
+      }
+    }
+
+    // erase discared particles in case of multithinning
+    if (eraseParticles_) {
+      for (auto& p : view) {
+        if (auto const w = p.getWeight(); w == 0) { p.erase(); }
+      }
+    } else {
+      return;
+    }
+  }
+
+} // namespace corsika

corsika/detail/modules/tracking/Intersect.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/core/Logging.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+
+#include <limits>
+
+namespace corsika {
+
+  template <typename TDerived>
+  template <typename TParticle>
+  inline auto Intersect<TDerived>::nextIntersect(TParticle const& particle,
+                                                 TimeType const step_limit) const {
+
+    typedef typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
+
+    node_type& volumeNode =
+        *particle.getNode(); // current "logical" node, from previous tracking step
+
+    CORSIKA_LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
+                      volumeNode.getVolume().contains(particle.getPosition()));
+
+    // start values:
+    TimeType minTime = step_limit;
+    node_type* minNode = &volumeNode;
+
+    // determine the first geometric collision with any other Volume boundary
+
+    // first check, where we leave the current volume
+    // this assumes our convention, that all Volume primitives must be convex
+    // thus, the last entry is always the exit point
+    Intersections const time_intersections_curr =
+        TDerived::intersect(particle, volumeNode);
+    CORSIKA_LOG_TRACE("curr node {}, parent node {}, hasIntersections={} ",
+                      fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
+                      time_intersections_curr.hasIntersections());
+    if (time_intersections_curr.hasIntersections()) {
+      CORSIKA_LOG_DEBUG(
+          "intersection times with currentLogicalVolumeNode: entry={} s and exit={} s",
+          time_intersections_curr.getEntry() / 1_s,
+          time_intersections_curr.getExit() / 1_s);
+      if (time_intersections_curr.getExit() <= minTime) {
+        minTime =
+            time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
+        minNode = volumeNode.getParent();
+      }
+    }
+
+    // where do we collide with any of the next-tree-level volumes
+    // entirely contained by currentLogicalVolumeNode
+    for (auto const& node : volumeNode.getChildNodes()) {
+
+      Intersections const time_intersections = TDerived::intersect(particle, *node);
+      CORSIKA_LOG_DEBUG("search intersection times with child volume {}:",
+                        fmt::ptr(node));
+      if (!time_intersections.hasIntersections()) { continue; }
+      auto const t_entry = time_intersections.getEntry();
+      auto const t_exit = time_intersections.getExit();
+      CORSIKA_LOG_DEBUG("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+                        t_entry <= minTime);
+      // note, theoretically t can even be smaller than 0 since we
+      // KNOW we can't yet be in this volume yet, so we HAVE TO
+      // enter it IF exit point is not also in the "past", AND if
+      // extry point is [[much much]] closer than exit point
+      // (because we might have already numerically "exited" it)!
+      if (t_exit > 0_s && t_entry <= minTime &&
+          -t_entry < t_exit) { // protection agains numerical problem, when we already
+                               // _exited_ before
+                               // enter chile volume here
+        minTime = t_entry;
+        minNode = node.get();
+      }
+    }
+
+    // these are volumes from the previous tree-level that are cut-out partly from the
+    // current volume
+    for (node_type* node : volumeNode.getExcludedNodes()) {
+
+      Intersections const time_intersections = TDerived::intersect(particle, *node);
+      if (!time_intersections.hasIntersections()) { continue; }
+      CORSIKA_LOG_DEBUG(
+          "intersection times with exclusion volume {} : enter {} s, exit {} s",
+          fmt::ptr(node), time_intersections.getEntry() / 1_s,
+          time_intersections.getExit() / 1_s);
+      auto const t_entry = time_intersections.getEntry();
+      auto const t_exit = time_intersections.getExit();
+      CORSIKA_LOG_DEBUG("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+                        t_entry <= minTime);
+      // note, theoretically t can even be smaller than 0 since we
+      // KNOW we can't yet be in this volume yet, so we HAVE TO
+      // enter it IF exit point is not also in the "past"!
+      if (t_exit > 0_s && t_entry <= minTime) { // enter volumen child here
+        minTime = t_entry;
+        minNode = node;
+      }
+    }
+    CORSIKA_LOG_DEBUG("next time-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
+    // this branch cannot be unit-testes. This is malfunction: LCOV_EXCL_START
+    if (minTime < 0_s) {
+      if (minTime < -1e-8_s) {
+        CORSIKA_LOG_DEBUG(
+            "There is a very negative time step detected: {}. This is not physical and "
+            "may "
+            "easily crash subsequent modules. Set to 0_s, but CHECK AND FIX.",
+            minTime);
+      }
+      minTime = 0_s; // LCOV_EXCL_STOP
+    }
+    return std::make_tuple(minTime, minNode);
+  }
+
+} // namespace corsika

corsika/detail/modules/tracking/TrackingLeapFrogCurved.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/modules/tracking/TrackingStraight.hpp> // for neutral particles
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/LeapFrogTrajectory.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/modules/tracking/Intersect.hpp>
+
+#include <type_traits>
+#include <utility>
+
+namespace corsika {
+
+  namespace tracking_leapfrog_curved {
+
+    template <typename TParticle>
+    inline auto Tracking::getTrack(TParticle const& particle) {
+      VelocityVector const initialVelocity = particle.getVelocity();
+
+      auto const& position = particle.getPosition();
+      CORSIKA_LOG_DEBUG(
+          "TrackingLeapfrog_Curved pid: {}"
+          " , E = {} GeV \n"
+          "\tTracking pos: {} \n"
+          "\tTracking   p: {} GeV \n"
+          "\tTracking   v: {}",
+          particle.getPID(), particle.getEnergy() / 1_GeV, position.getCoordinates(),
+          particle.getMomentum().getComponents() / 1_GeV,
+          initialVelocity.getComponents());
+
+      typedef
+          typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
+      node_type const& volumeNode = *particle.getNode();
+
+      // for the event of magnetic fields and curved trajectories, we need to limit
+      // maximum step-length since we need to follow curved
+      // trajectories segment-wise -- at least if we don't employ concepts as "Helix
+      // Trajectories" or similar
+      MagneticFieldVector const& magneticfield =
+          volumeNode.getModelProperties().getMagneticField(position);
+      MagneticFluxType const magnitudeB = magneticfield.getNorm();
+      ElectricChargeType const charge = particle.getCharge();
+      bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
+
+      if (no_deflection) {
+        CORSIKA_LOG_TRACE("no_deflection");
+        return getLinearTrajectory(particle);
+      }
+
+      HEPMomentumType const p_perp =
+          (particle.getMomentum() -
+           particle.getMomentum().getParallelProjectionOnto(magneticfield))
+              .getNorm();
+
+      CORSIKA_LOG_TRACE("p_perp={} eV", p_perp / 1_eV);
+
+      if (p_perp < 1_eV) {
+        // particle travel along, parallel to magnetic field. Rg is
+        // "0", but for purpose of step limit we return infinity here.
+        CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
+        return getLinearTrajectory(particle);
+      }
+
+      LengthType const gyroradius = convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+                                    constants::c / (abs(charge) * magnitudeB);
+
+      if (gyroradius > 1e9_m) {
+        // this cannot be really unit-tested. It is hidden. LCOV_EXCL_START
+        CORSIKA_LOG_TRACE(
+            "CurvedLeapFrog is not very stable for extremely high gyroradius steps. "
+            "Rg={} -> straight tracking.",
+            gyroradius);
+        return getLinearTrajectory(particle);
+        // LCOV_EXCL_STOP
+      }
+
+      LengthType const steplimit = 2 * cos(maxMagneticDeflectionAngle_) *
+                                   sin(maxMagneticDeflectionAngle_) * gyroradius;
+      TimeType const steplimit_time = steplimit / initialVelocity.getNorm();
+      CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
+                        steplimit_time);
+
+      // traverse the environment volume tree and find next
+      // intersection
+      auto [minTime, minNode] = nextIntersect(particle, steplimit_time);
+
+      auto const p_norm =
+          constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+                             particle.getMomentum().getNorm()); // kg *m /s
+      // k = q/|p|
+      decltype(1 / (tesla * second)) const k =
+          charge / p_norm *
+          initialVelocity.getNorm(); // since we use steps in time and not length
+      // units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1 / (T*s)
+
+      return std::make_tuple(
+          LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
+                             minTime), // --> trajectory
+          minNode);                    // --> next volume node
+    }
+
+    template <typename TParticle>
+    inline Intersections Tracking::intersect(TParticle const& particle,
+                                             Sphere const& sphere) {
+
+      LengthType const radius = sphere.getRadius();
+      if (radius == 1_km * std::numeric_limits<double>::infinity()) {
+        return Intersections();
+      }
+
+      ElectricChargeType const charge = particle.getCharge();
+      auto const& position = particle.getPosition();
+      auto const* currentLogicalVolumeNode = particle.getNode();
+      MagneticFieldVector const& magneticfield =
+          currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
+
+      VelocityVector const velocity = particle.getVelocity();
+      DirectionVector const directionBefore = velocity.normalized();
+
+      auto const projectedDirection = directionBefore.cross(magneticfield);
+      auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm();
+      bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
+
+      CORSIKA_LOG_TRACE("projectedDirectionSqrNorm={} T^2",
+                        projectedDirectionSqrNorm / square(1_T));
+      if (isParallel) {
+        // particle moves parallel to field -> no deflection
+        return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
+      }
+
+      bool const numericallyInside = sphere.contains(particle.getPosition());
+      CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
+
+      Vector<length_d> const deltaPos = position - sphere.getCenter();
+
+      { // check extreme cases we don't want to solve analytically explicit
+        HEPMomentumType const p_perp =
+            (particle.getMomentum() -
+             particle.getMomentum().getParallelProjectionOnto(magneticfield))
+                .getNorm();
+
+        LengthType const gyroradius =
+            (convert_HEP_to_SI<MassType::dimension_type>(p_perp) * constants::c /
+             (abs(charge) * magneticfield.getNorm()));
+
+        LengthType const trackDist = abs(deltaPos.getNorm() - radius);
+        if (trackDist > gyroradius) {
+          // there is never a solution
+          return Intersections();
+        }
+
+        if (gyroradius > 1000 * trackDist) {
+          // the bending is negligible, use straight intersections instead
+          return tracking_line::Tracking::intersect(particle, sphere);
+        }
+      }
+
+      SpeedType const absVelocity = velocity.getNorm();
+      auto const p_norm =
+          constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+                             particle.getMomentum().getNorm()); // km * m /s
+      // this is: k = q/|p|
+      decltype(1 / (tesla * meter)) const k = charge / p_norm;
+
+      MagneticFieldVector const direction_x_B = directionBefore.cross(magneticfield);
+      auto const denom = 4. / (direction_x_B.getSquaredNorm() * k * k);
+      double const b = (direction_x_B.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
+      double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m);
+      LengthType const deltaPosLength = deltaPos.getNorm();
+      double const d = (deltaPosLength + radius) * (deltaPosLength - radius) * denom /
+                       (1_m * 1_m * 1_m * 1_m);
+      CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
+      // solutions of deltaL are obtained from quartic equation. Note, deltaL/2 is the
+      // length of each half step, however, the second half step is slightly longer
+      // because of the non-conservation of norm/velocity.
+      // The leap-frog length L is deltaL/2 * (1+|u_{n+1}|)
+      std::vector<double> solutions = solve_quartic_real(1, 0, b, c, d);
+      if (!solutions.size()) { return Intersections(); }
+      LengthType d_enter, d_exit;
+      int first = 0, first_entry = 0, first_exit = 0;
+      for (auto const solution : solutions) {
+        LengthType const dist = solution * 1_m;
+        CORSIKA_LOG_TRACE(
+            "Solution (real) for current Volume: deltaL/2*2={} (deltaL/2*2/v={}) ", dist,
+            dist / absVelocity);
+        if (numericallyInside) {
+          // there must be an entry (negative) and exit (positive) solution
+          if (dist < 0.0001_m) { // security margin to assure
+                                 // transfer to next logical volume
+                                 // (even if dist suggest marginal
+                                 // entry already, which we
+                                 // classify as numerical artifact)
+            if (first_entry == 0) {
+              d_enter = dist;
+            } else {
+              d_enter = std::max(d_enter, dist); // closest negative to zero >1e-4 m
+            }
+            first_entry++;
+
+          } else { // thus, dist > +0.0001_m
+
+            if (first_exit == 0) {
+              d_exit = dist;
+            } else {
+              d_exit = std::min(d_exit, dist); // closest positive to zero >1e-4 m
+            }
+            first_exit++;
+          }
+          first = int(first_exit > 0) + int(first_entry > 0);
+
+        } else { // thus, numericallyInside == false
+
+          // both physical solutions (entry, exit) must be positive, and as small as
+          // possible
+          if (dist < -0.0001_m) { // need small numerical margin, to
+                                  // assure transport. We consider
+                                  // begin marginally already inside
+                                  // next volume (besides
+                                  // numericallyInside=false) as numerical glitch.
+            // into next logical volume
+            continue;
+          }
+          if (first == 0) {
+            d_enter = dist;
+          } else {
+            if (dist < d_enter) {
+              d_exit = d_enter;
+              d_enter = dist;
+            } else {
+              d_exit = dist;
+            }
+          }
+          first++;
+        }
+      } // loop over solutions
+
+      if (first == 0) { // entry and exit points found
+        CORSIKA_LOG_DEBUG(
+            "no intersections found: count={}, first_entry={}, first_exit={}", first,
+            first_entry, first_exit);
+        return Intersections();
+      }
+      // return in units of time
+
+      return Intersections(d_enter / absVelocity, d_exit / absVelocity);
+    }
+
+    template <typename TParticle>
+    inline Intersections Tracking::intersect(TParticle const& particle,
+                                             Plane const& plane) {
+
+      CORSIKA_LOG_TRACE("intersection particle with plane");
+
+      ElectricChargeType const charge = particle.getCharge();
+
+      if (charge != ElectricChargeType::zero()) {
+
+        auto const* currentLogicalVolumeNode = particle.getNode();
+        VelocityVector const velocity = particle.getVelocity();
+        auto const absVelocity = velocity.getNorm();
+        DirectionVector const direction = velocity.normalized();
+        Point const& position = particle.getPosition();
+
+        auto const magneticfield =
+            currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
+
+        // solve:     denom x^2 + p x + q =0    for     x = delta-l
+
+        auto const direction_x_B = direction.cross(magneticfield);
+        double const denom = charge *
+                             plane.getNormal().dot(direction_x_B) // unit: C*T = kg/s
+                             / 1_kg * 1_s;
+
+        CORSIKA_LOG_TRACE("denom={}", denom);
+
+        auto const p_norm =
+            constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+                               particle.getMomentum().getNorm()); // unit: kg * m/s
+
+        double const p = (2 * p_norm * direction.dot(plane.getNormal())) // unit: kg*m/s
+                         / (1_m * 1_kg) * 1_s;
+        double const q =
+            (2 * p_norm *
+             plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m
+            / (1_m * 1_m * 1_kg) * 1_s;
+
+        // deltaL from quadratic solution return half-step length deltaL/2 for leap-frog
+        // algorithmus. Note, the leap-frog length L is longer by (1+|u_{n_1}|)/2 because
+        // the direction norm of the second half step is >1.
+
+        std::vector<double> const deltaLs = solve_quadratic_real(denom, p, q);
+
+        CORSIKA_LOG_TRACE("deltaLs=[{}]", fmt::join(deltaLs, ", "));
+
+        if (deltaLs.size() == 0) {
+          return Intersections(std::numeric_limits<double>::infinity() * 1_s);
+        }
+
+        // select smallest but positive solution
+        bool first = true;
+        LengthType maxStepLength = 0_m;
+        for (auto const& deltaL : deltaLs) {
+          if (deltaL < 0) continue;
+          if (first) {
+            first = false;
+            maxStepLength = deltaL * meter;
+          } else if (maxStepLength > deltaL * meter) {
+            maxStepLength = deltaL * meter;
+          }
+        }
+
+        // check: both intersections in past, or no valid intersection
+        if (first) {
+          return Intersections(std::numeric_limits<double>::infinity() * 1_s);
+        }
+
+        CORSIKA_LOG_TRACE("maxStepLength={} s", maxStepLength / 1_s);
+
+        // with final length correction, |direction| becomes >1 during step
+
+        return Intersections(maxStepLength / absVelocity); // unit: s
+
+      } // no charge
+
+      CORSIKA_LOG_TRACE("(plane) straight tracking with  charge={}, B={}", charge,
+                        particle.getNode()->getModelProperties().getMagneticField(
+                            particle.getPosition()));
+
+      return tracking_line::Tracking::intersect(particle, plane);
+    }
+
+    template <typename TParticle>
+    inline Intersections Tracking::intersect(TParticle const& particle,
+                                             SeparationPlane const& sepPlane) {
+      return intersect(particle, sepPlane.getPlane());
+    }
+
+    template <typename TParticle, typename TBaseNodeType>
+    inline Intersections Tracking::intersect(TParticle const& particle,
+                                             TBaseNodeType const& volumeNode) {
+      Sphere const* sphere = dynamic_cast<Sphere const*>(&volumeNode.getVolume());
+      if (sphere) { return intersect(particle, *sphere); }
+      SeparationPlane const* sepPlane =
+          dynamic_cast<SeparationPlane const*>(&volumeNode.getVolume());
+      if (sepPlane) { return intersect(particle, *sepPlane); }
+      throw std::runtime_error(
+          "The Volume type provided is not supported in "
+          "TrackingLeapFrogCurved::intersect(particle, node)");
+    }
+
+    template <typename TParticle>
+    inline auto Tracking::getLinearTrajectory(TParticle& particle) {
+
+      // perform simple linear tracking
+      auto [straightTrajectory, minNode] = straightTracking_.getTrack(particle);
+
+      // return as leap-frog trajectory
+      return std::make_tuple(
+          LeapFrogTrajectory(
+              straightTrajectory.getLine().getStartPoint(),
+              straightTrajectory.getLine().getVelocity(),
+              MagneticFieldVector(particle.getPosition().getCoordinateSystem(), 0_T, 0_T,
+                                  0_T),
+              0 * square(meter) / (square(second) * volt),
+              straightTrajectory.getDuration()), // trajectory
+          minNode);                              // next volume node
+    }
+
+  } // namespace tracking_leapfrog_curved
+
+} // namespace corsika