corsika/modules/radio/CoREAS.hpp

+/*
+ * (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/RadioProcess.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  template <typename TRadioDetector, typename TPropagator>
+  class CoREAS final
+      : public RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>,
+                            TPropagator> {
+
+  public:
+    // an identifier for which algorithm was used
+    static constexpr auto algorithm = "CoREAS";
+
+    /**
+     * Construct a new CoREAS instance.
+     *
+     * This forwards the detector and other arguments to
+     * the RadioProcess parent.
+     *
+     */
+    template <typename... TArgs>
+    CoREAS(TRadioDetector& detector, TArgs&&... args)
+        : RadioProcess<TRadioDetector, CoREAS, TPropagator>(detector, args...) {}
+
+    /**
+     * Simulate the radio emission from a particle across a track.
+     *
+     * This must be provided by the TRadioImpl.
+     *
+     * @param particle    The current particle.
+     * @param track       The current track.
+     *
+     */
+    template <typename Particle>
+    ProcessReturn simulate(Step<Particle> const& step);
+    static constexpr auto emConstant_{1.0 / (4.0 * M_PI) / (constants::epsilonZero) /
+                                      constants::c};
+
+    using Base =
+        RadioProcess<TRadioDetector, CoREAS<TRadioDetector, TPropagator>, TPropagator>;
+    using Base::observers_;
+
+  }; // end of class CoREAS
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/CoREAS.inl>

corsika/modules/radio/LimitedRadioProcess.hpp

+/*
+ * (c) Copyright 2024 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/modules/radio/RadioProcess.hpp>
+
+namespace corsika {
+  /*!
+   * Wrapper for a radio process which turns the production of signal on and off depending
+   * on where the currently-emitting particle is in space. A function with signature `bool
+   * Fcn(Point const&)` must be given which turns this on and off
+   */
+
+  template <class TUnderlyingProcess>
+  class LimitedRadioProcess
+      : public ContinuousProcess<LimitedRadioProcess<TUnderlyingProcess>> {
+
+    //! define the default constructor for the function signature
+    static auto constexpr default_functor = [](Point const&) { return true; };
+
+  public:
+    //! The signature of the if-statement that decides if this RadioProcess should be
+    //! called
+    using functor_signature = bool(Point const&);
+
+    LimitedRadioProcess(TUnderlyingProcess&& process,
+                        std::function<functor_signature> runThis);
+
+    // This function wraps the one from RadioProcess and runs it depending on output of
+    // `runThis_`
+    template <typename Particle>
+    ProcessReturn doContinuous(Step<Particle> const& step, bool const);
+
+    // This process just returns the one from RadioProcess since it would be more work to
+    // turn this on and off than to just feed-through
+    template <typename Particle, typename Track>
+    LengthType getMaxStepLength(Particle const& vParticle, Track const& vTrack) const;
+
+  private:
+    TUnderlyingProcess process_;
+    std::function<functor_signature> const runThis_{default_functor};
+  };
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/LimitedRadioProcess.inl>

corsika/modules/radio/RadioProcess.hpp

+/*
+ * (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/output/BaseOutput.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+#include <corsika/framework/process/ContinuousProcess.hpp>
+#include <corsika/framework/core/Step.hpp>
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+namespace corsika {
+
+  /**
+   * The base interface for radio emission processes.
+   *
+   * TRadioImpl is the concrete implementation of the radio algorithm.
+   * TObserverCollection is the detector instance that stores observers
+   * and is responsible for managing the output writing.
+   */
+  template <typename TObserverCollection, typename TRadioImpl, typename TPropagator>
+  class RadioProcess : public ContinuousProcess<
+                           RadioProcess<TObserverCollection, TRadioImpl, TPropagator>>,
+                       public BaseOutput {
+
+    /*
+     * A collection of filter objects for deciding on valid particles and tracks.
+     */
+    // std::vector<std::function<bool(ParticleType&, TrackType const&)>> filters_;
+
+    /**
+     * Get a reference to the underlying radio implementation.
+     */
+    TRadioImpl& implementation();
+
+    /**
+     *  Get a const reference to the underlying implementation.
+     */
+    TRadioImpl const& implementation() const;
+
+  protected:
+    TObserverCollection& observers_; ///< The radio observers we store into.
+    TPropagator propagator_;         ///< The propagator implementation.
+    unsigned int showerId_{0};       ///< The current event ID.
+    ParquetStreamer output_;         //!< The parquet streamer for this process.
+
+  public:
+    using axistype = std::vector<long double>;
+    /**
+     * Construct a new RadioProcess.
+     */
+    RadioProcess(TObserverCollection& observers, TPropagator& propagator);
+
+    /**
+     * Perform the continuous process (radio emission).
+     *
+     * This handles filtering individual particle tracks
+     * before passing them to `Simulate`.`
+     *
+     * @param particle    The current particle.
+     * @param track       The current track.
+     */
+    template <typename Particle>
+    ProcessReturn doContinuous(Step<Particle> const& step, bool const);
+
+    /**
+     * Return the maximum step length for this particle and track.
+     *
+     * This must be provided by the TRadioImpl.
+     *
+     * @param particle    The current particle.
+     * @param track       The current track.
+     *
+     * @returns The maximum length of this track.
+     */
+    template <typename Particle, typename Track>
+    LengthType getMaxStepLength(Particle const& vParticle, Track const& vTrack) const;
+
+    /**
+     * Called at the start of each library.
+     */
+    void startOfLibrary(boost::filesystem::path const& directory) final override;
+
+    /**
+     * Called at the end of each shower.
+     */
+    virtual void endOfShower(unsigned int const) final override;
+
+    /**
+     * Called at the end of each library.
+     *
+     */
+    void endOfLibrary() final override {}
+
+    /**
+     * Get the configuration of this output.
+     */
+    YAML::Node getConfig() const final;
+
+  }; // END: class RadioProcess
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/RadioProcess.inl>

corsika/modules/radio/ZHS.hpp

+/*
+ * (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/RadioProcess.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  /**
+   * A concrete implementation of the ZHS algorithm.
+   */
+  template <typename TRadioDetector, typename TPropagator>
+  class ZHS final : public RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>,
+                                        TPropagator> {
+
+  public:
+    // an identifier for which algorithm was used
+    static constexpr auto algorithm = "ZHS";
+
+    /**
+     * Construct a new ZHS instance.
+     *
+     * This forwards the detector and other arguments to
+     * the RadioProcess parent.
+     *
+     */
+    template <typename... TArgs>
+    ZHS(TRadioDetector& detector, TArgs&&... args)
+        : RadioProcess<TRadioDetector, ZHS, TPropagator>(detector, args...) {}
+
+    /**
+     * Simulate the radio emission from a particle across a track.
+     *
+     * This must be provided by the TRadioImpl.
+     *
+     * @param particle    The current particle.
+     * @param track       The current track.
+     *
+     */
+    template <typename Particle>
+    ProcessReturn simulate(Step<Particle> const& step);
+    static constexpr auto emConstant_{1.0 / (4.0 * M_PI) / (constants::epsilonZero) /
+                                      constants::c};
+
+  private:
+    using Base =
+        RadioProcess<TRadioDetector, ZHS<TRadioDetector, TPropagator>, TPropagator>;
+    using Base::observers_;
+
+  }; // END: class ZHS
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/ZHS.inl>

corsika/modules/radio/detectors/ObserverCollection.hpp

+/*
+ * (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 {
+
+  /**
+   * The base interface for radio detectors.
+   * At the moment it is a collection of observers with the same implementation.
+   */
+
+  template <typename TObserverImpl>
+  class ObserverCollection {
+
+    /**
+     * The collection of observers used in this simulation.
+     */
+    std::vector<TObserverImpl> observers_;
+
+  public:
+    /**
+     * Add an observer to this radio process.
+     *
+     * @param observer    The observer to add
+     */
+    void addObserver(TObserverImpl const& observer);
+
+    /**
+     * Get the specific observer at that place in the collection
+     *
+     * @param index in the collection
+     */
+    TObserverImpl& at(std::size_t const i);
+
+    TObserverImpl const& at(std::size_t const i) const;
+
+    /**
+     * Get the number of observerss in the collection
+     */
+    int size() const;
+
+    /**
+     * Get a *non*-const reference to the collection of observers.
+     *
+     * @returns    An iterable mutable reference to the observers.
+     */
+    std::vector<TObserverImpl>& getObservers();
+
+    /**
+     * Get a const reference to the collection of observers.
+     *
+     * @returns    An iterable mutable reference to the observers.
+     */
+    std::vector<TObserverImpl> const& getObservers() const;
+
+    /**
+     * Reset all the observer waveforms.
+     */
+    void reset();
+  }; // END: class RadioDetector
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/detectors/ObserverCollection.inl>

corsika/modules/radio/observers/Observer.hpp

+/*
+ * (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 <boost/filesystem.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalGeometry.hpp>
+#include <corsika/output/ParquetStreamer.hpp>
+
+namespace corsika {
+
+  /**
+   * A common abstract interface for radio oberservers.
+   *
+   * All concrete observer implementations should be of
+   * type Observer<T> where T is a concrete observer implementation.
+   *
+   */
+  template <typename TObserverImpl>
+  class Observer {
+
+  protected:
+    std::string const name_; ///< The name/identifier of this observer.
+    Point const location_;   ///< The location of this observer.
+    CoordinateSystemPtr const
+        coordinateSystem_; ///< The coordinate system of the observer
+
+  public:
+    using axistype = std::vector<long double>;
+
+    /**
+     * \brief Construct a base observer instance.
+     *
+     * @param name    A name for this observer.
+     * @param location    The location of this observer.
+     *
+     */
+    Observer(std::string const& name, Point const& location,
+             CoordinateSystemPtr const& coordinateSystem);
+
+    /**
+     * Receive a signal at this observer.
+     *
+     * This is a general implementation call that must be specialized
+     * for the particular observer implementation and usage.
+     *
+     */
+    template <typename... TVArgs>
+    void receive(TVArgs&&... args);
+
+    /**
+     * Get the location of this observer.
+     */
+    Point const& getLocation() const;
+
+    /**
+     * Get the name of this name observer.
+     *
+     * This is used in producing the output data file.
+     */
+    std::string const& getName() const;
+
+    /**
+     * Reset the observer before starting a new simulation.
+     */
+    void reset();
+
+    /**
+     * Return a reference to the x-axis labels (i.e. time or frequency).
+     *
+     * This should be an xtensor-convertible type with
+     * a ->data() method that converts to a raw pointer.
+     */
+    axistype getAxis() const;
+
+    /**
+     * Return a reference to the underlying waveform data for X polarization.
+     *
+     * This is used when writing the observer information to disk
+     * and will be converted to a 32-bit float before writing.
+     */
+    std::vector<double> const& getWaveformX() const;
+
+    /**
+     * Return a reference to the underlying waveform data for Y polarization.
+     *
+     * This is used when writing the observer information to disk
+     * and will be converted to a 32-bit float before writing.
+     */
+    std::vector<double> const& getWaveformY() const;
+
+    /**
+     * Return a reference to the underlying waveform data for Z polarization.
+     *
+     * This is used when writing the observer information to disk
+     * and will be converted to a 32-bit float before writing.
+     */
+    std::vector<double> const& getWaveformZ() const;
+
+    /**
+     * Get a reference to the underlying radio implementation.
+     */
+    TObserverImpl& implementation();
+
+  }; // END: class Observer final
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/observers/Observer.inl>

corsika/modules/radio/observers/TimeDomainObserver.hpp

+/*
+ * (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>
+#include <yaml-cpp/yaml.h>
+#include <vector>
+
+namespace corsika {
+
+  /**
+   * An implementation of a time-domain observer that has a customized
+   * start time, sampling rate, and waveform duration.
+   *
+   */
+  class TimeDomainObserver : public Observer<TimeDomainObserver> {
+
+    TimeType const start_time_;         ///< The start time of this waveform.
+    TimeType const duration_;           ///< The duration of this waveform.
+    InverseTimeType const sample_rate_; ///< The sampling rate of this observer.
+    TimeType const ground_hit_time_; ///< The time the primary particle hits the ground.
+    uint64_t const num_bins_;        ///< The number of bins used.
+    std::vector<double> waveformEX_; ///< EX polarization.
+    std::vector<double> waveformEY_; ///< EY polarization.
+    std::vector<double> waveformEZ_; ///< EZ polarization.
+    std::vector<long double> const
+        time_axis_; ///< The time axis corresponding to the electric field.
+
+  public:
+    // import the methods from the observer
+    using Observer<TimeDomainObserver>::getName;
+    using Observer<TimeDomainObserver>::getLocation;
+
+    /**
+     * Construct a new TimeDomainObserver.
+     *
+     * @param name               The name of this observer.
+     * @param location           The location of this observer.
+     * @param coordinateSystem   The coordinate system of this observer.
+     * @param start_time         The starting time of this waveform.
+     * @param duration           The duration of this waveform.
+     * @param sample_rate        The sample rate of this waveform.
+     * @param ground_hit_time    The time the primary particle hits the ground on a
+     * straight vertical line.
+     *
+     */
+    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);
+
+    /**
+     * Receive an electric field at this observer.
+     *
+     * This assumes that the observer will receive
+     *  an *instantaneous* electric field modeled as a delta function (or timebin).
+     *
+     * @param time             The (global) time at which this signal is received.
+     * @param receive_vector   The incident unit vector. (not used at the moment)
+     * @param field            The incident electric field vector.
+     *
+     */
+    // TODO: rethink this method a bit. If the endpoint is at the end of the observer
+    // resolution then you get the startpoint signal but you lose the endpoint signal!
+    void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
+                 ElectricFieldVector const& efield);
+
+    void receive(TimeType const time, Vector<dimensionless_d> const& receive_vector,
+                 VectorPotential const& vectorP);
+
+    /**
+     * Return the time-units of each waveform for X polarization
+     *
+     * This returns them in nanoseconds for ease of use.
+     */
+    auto const& getWaveformX() const;
+
+    /**
+     * Return the time-units of each waveform for Y polarization
+     *
+     * This returns them in nanoseconds for ease of use.
+     */
+    auto const& getWaveformY() const;
+
+    /**
+     * Return the time-units of each waveform for Z polarization
+     *
+     * This returns them in nanoseconds for ease of use.
+     */
+    auto const& getWaveformZ() const;
+
+    /**
+     * Return a label that indicates that this is a time
+     * domain observer
+     *
+     * This returns the string "Time".
+     */
+    std::string const getDomainLabel();
+
+    /**
+     * Creates time-units of each waveform.
+     *
+     * It creates them in nanoseconds for ease of use.
+     */
+    std::vector<long double> createTimeAxis() const;
+
+    /**
+     * Return the time-units of each waveform.
+     *
+     * This returns them in nanoseconds for ease of use.
+     */
+    auto const getAxis() const;
+
+    /**
+     * Returns the sampling rate of the time domain observer.
+     */
+    InverseTimeType const& getSampleRate() const;
+
+    /**
+     * Returns the start time of detection for the time domain observer.
+     */
+    TimeType const& getStartTime() const;
+
+    /**
+     * Reset the observer before starting a new simulation.
+     */
+    void reset();
+
+    /**
+     * Return a YAML configuration for this observer.
+     */
+    YAML::Node getConfig() const;
+
+  }; // END: class TimeDomainObserver
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/observers/TimeDomainObserver.inl>

corsika/modules/radio/propagators/DummyTestPropagator.hpp

+/*
+ * (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/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+  /**
+   * This class implements a dummy propagator that uses
+   * the straight-line (vector) between the particle
+   * location and the observer as the trajectory.
+   * It is intended mainly for fast testing as it only
+   * works with 2 points in a uniform refractive index
+   * atmospheric profile.
+   *
+   */
+  template <typename TEnvironment>
+  class DummyTestPropagator final
+      : public RadioPropagator<DummyTestPropagator<TEnvironment>, TEnvironment> {
+
+    using Base = RadioPropagator<DummyTestPropagator<TEnvironment>, TEnvironment>;
+    using SignalPathCollection = typename Base::SignalPathCollection;
+
+  public:
+    /**
+     * Construct a new SimplePropagator with a given environment.
+     *
+     */
+    DummyTestPropagator(TEnvironment const& env);
+
+    /**
+     * Return the collection of paths from `source` to `destination`.
+     * Hence, the signal propagated from the
+     * emission point to the observer location.
+     *
+     */
+    template <typename Particle>
+    SignalPathCollection propagate(Particle const& particle, Point const& source,
+                                   Point const& destination);
+
+  private:
+    std::deque<Point> points;
+    std::vector<double> rindex;
+
+  }; // End: SimplePropagator
+
+  template <typename TEnvironment>
+  DummyTestPropagator<TEnvironment> make_dummy_test_radio_propagator(
+      TEnvironment const& env) {
+    return DummyTestPropagator<TEnvironment>(env);
+  }
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/DummyTestPropagator.inl>

corsika/modules/radio/propagators/NumericalIntegratingPropagator.hpp

+/*
+ * (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/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+  /**
+   * This class implements a basic propagator that uses
+   * the straight-line (vector) between the particle
+   * location and the observer as the trajectory.
+   * To calculate the time delay of the signal, a basic
+   * numerical integration scheme based on Simpson's rule
+   * takes place. This propagator is slow and not
+   * recommended for big showers simulations.
+   *
+   * This is what is used in ZHAireS and CoREAS in C7.
+   */
+  template <typename TEnvironment>
+  class NumericalIntegratingPropagator final
+      : public RadioPropagator<NumericalIntegratingPropagator<TEnvironment>,
+                               TEnvironment> {
+
+    using Base =
+        RadioPropagator<NumericalIntegratingPropagator<TEnvironment>, TEnvironment>;
+    using SignalPathCollection = typename Base::SignalPathCollection;
+
+  public:
+    /**
+     * Construct a new StraightPropagator with a given environment.
+     *
+     */
+    NumericalIntegratingPropagator(TEnvironment const& env, LengthType const stepsize);
+
+    /**
+     * Return the collection of paths from `start` to `end`.
+     * or from 'source' which is the emission point to 'destination'
+     * which is the location of the observer
+     */
+    template <typename Particle>
+    SignalPathCollection propagate(Particle const& particle, Point const& source,
+                                   Point const& destination) const;
+
+  private:
+    LengthType const stepsize_;
+
+  }; // End: StraightPropagator
+
+  template <typename TEnvironment>
+  NumericalIntegratingPropagator<TEnvironment>
+  make_numerical_integrating_radio_propagator(TEnvironment const& env,
+                                              LengthType const stepsize) {
+    return NumericalIntegratingPropagator<TEnvironment>(env, stepsize);
+  }
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/NumericalIntegratingPropagator.inl>

corsika/modules/radio/propagators/RadioPropagator.hpp

+/*
+ * (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 <vector>
+
+#include <corsika/framework/geometry/Path.hpp>
+#include <corsika/modules/radio/propagators/SignalPath.hpp>
+
+namespace corsika {
+
+  /**
+   * Radio propagators are used to calculate the propagation
+   * paths from particles to observers. Any class that wants
+   * to be used as a RadioPropagator must implement the
+   * following methods:
+   *
+   *  SignalPathCollection Propagate(Particle const& particle,
+   *                                 Point const& source,
+                                     Point const& destination) const
+   */
+  template <typename TImpl, typename TEnvironment>
+  class RadioPropagator {
+
+  protected:
+    // Since we typically know roughly how many paths will
+    // be computed for a given propagator, we use a std::vector here.
+    using SignalPathCollection = std::vector<SignalPath> const;
+
+    TEnvironment const& env_; ///< The environment.
+
+  public:
+    /**
+     * Construct a new RadioPropagator instance.
+     */
+    RadioPropagator(TEnvironment const& env);
+
+  }; // class RadioPropagator
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/RadioPropagator.inl>

corsika/modules/radio/propagators/SignalPath.hpp

+/*
+ * (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/Path.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  /**
+   * Store the photon signal path between two points.
+   *
+   * This is basically a container class
+   */
+  struct SignalPath final : public Path {
+
+    // TODO: discuss if we need average refractivity or average refractive index
+    TimeType const propagation_time_;       ///< The total propagation time.
+    double const average_refractive_index_; ///< The average refractive index.
+    double const refractive_index_source_;  ///< The refractive index at the source.
+    double const
+        refractive_index_destination_; ///< The refractive index at the destination point.
+    Vector<dimensionless_d> const emit_;    ///< The (unit-length) emission vector.
+    Vector<dimensionless_d> const receive_; ///< The (unit-length) receive vector.
+    LengthType const
+        R_distance_; ///< The distance from the point of emission to an observer. TODO:
+                     ///< optical path, not geometrical! (probably)
+
+    /**
+     * Create a new SignalPath instance.
+     */
+    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);
+
+  }; // END: class SignalPath final
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/SignalPath.inl>
\ No newline at end of file

corsika/modules/radio/propagators/TabulatedFlatAtmospherePropagator.hpp

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+#pragma once
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/PhysicalConstants.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/modules/radio/propagators/RadioPropagator.hpp>
+
+namespace corsika {
+
+  /**
+   * This class implements a tabulated propagator that approximates
+   * the Earth's atmosphere as flat. Signal propagation is rectilinear
+   * and this is intended to be used for vertical showers
+   * (<60 degrees zenith angle) for fast simulations. The table needs
+   * to have at least 10 elements.
+   *
+   */
+  template <typename TEnvironment>
+  class TabulatedFlatAtmospherePropagator final
+      : public RadioPropagator<TabulatedFlatAtmospherePropagator<TEnvironment>,
+                               TEnvironment> {
+
+    using Base =
+        RadioPropagator<TabulatedFlatAtmospherePropagator<TEnvironment>, TEnvironment>;
+    using SignalPathCollection = typename Base::SignalPathCollection;
+
+  public:
+    /**
+     * Construct a new FlatEarthPropagator with a given environment.
+     *
+     */
+    TabulatedFlatAtmospherePropagator(TEnvironment const& env, Point const& upperLimit,
+                                      Point const& lowerLimit, LengthType const step);
+
+    /**
+     * Return the collection of paths from `source` to `destination`.
+     * Hence, the signal propagated from the
+     * emission point to the oberserver location.
+     *
+     */
+    template <typename Particle>
+    SignalPathCollection propagate(Particle const& particle, Point const& source,
+                                   Point const& destination);
+
+  private:
+    Point const upperLimit_; ///< the upper point of the table.
+    Point const lowerLimit_; ///< the lowest point of the table (ideally earth's surface).
+    LengthType const step_;  ///< the tabulation step.
+    InverseLengthType const
+        inverseStep_;            ///< inverse of the step used to speed up calculations.
+    LengthType const maxHeight_; ///< z coordinate of upper limit (maximum height).
+    LengthType const minHeight_; ///< z coordinate of lower limit (minimum height) - 1_km
+                                 ///< for safety reasons.
+    std::vector<double> refractivityTable_; ///< the table that stores refractivity.
+    std::vector<LengthType>
+        heightTable_; ///< the table that stores the height using the step above.
+    std::vector<double>
+        integratedRefractivityTable_; ///< the table that stores integrated refractivity.
+    double slopeRefrLower_; ///< used to interpolate refractivity for particles below the
+                            ///< lowest point.
+    double slopeIntRefrLower_; ///< used to interpolate integrated refractivity for
+                               ///< particles below the lowest point.
+    double slopeRefrUpper_; ///< used to interpolate refractivity for particles above the
+                            ///< highest point.
+    double slopeIntRefrUpper_; ///< used to interpolate integrated refractivity for
+                               ///< particles above the highest point.
+    double lastElement_;       ///< last index of tables.
+    std::deque<Point> points;  ///< the points that the signal has propagated through.
+    std::vector<double>
+        rindex; ///< the refractive index values along the signal propagation path.
+
+  }; // End: FlatEarthPropagator
+
+  template <typename TEnvironment>
+  TabulatedFlatAtmospherePropagator<TEnvironment>
+  make_tabulated_flat_atmosphere_radio_propagator(TEnvironment const& env,
+                                                  Point const& upperLimit,
+                                                  Point const& lowerLimit,
+                                                  LengthType const step) {
+    return TabulatedFlatAtmospherePropagator<TEnvironment>(env, upperLimit, lowerLimit,
+                                                           step);
+  }
+
+} // namespace corsika
+
+#include <corsika/detail/modules/radio/propagators/TabulatedFlatAtmospherePropagator.inl>

corsika/modules/sibyll/Decay.hpp

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -75,6 +74,8 @@ namespace corsika::sibyll {
     bool handleAllDecays_ = true;
     bool sibyll_listing_ = false;
     std::set<Code> handledDecays_;
+
+    std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_sibyll_Decay");
   };
 
 } // namespace corsika::sibyll

corsika/modules/sibyll/HadronInteractionModel.hpp

+/*
+ * (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 <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+
+#include <set>
+#include <tuple>
+
+namespace corsika::sibyll {
+
+  /**
+   * @brief Provides the SIBYLL hadron-nucleus interaction model.
+   *
+   * This is a TModel argument for InteractionProcess<TModel>.
+   */
+
+  class HadronInteractionModel {
+
+  public:
+    HadronInteractionModel();
+    /**
+     * @brief construct hadron interaction model with SIBYLL
+     *
+     * @param list: list of particles that should be stable inside SIBYLL. all other
+     * particles should decay (if they can)
+     */
+    HadronInteractionModel(std::set<Code> list);
+    ~HadronInteractionModel();
+
+    /**
+     * @brief Set the Verbose flag.
+     *
+     * If flag is true, SIBYLL will printout additional secondary particle information
+     * lists, etc.
+     *
+     * @param flag to switch.
+     */
+    void setVerbose(bool const flag);
+
+    /**
+     * @brief evaluated validity of collision system.
+     *
+     * sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
+     * neutrons (p,n == nucleon).
+     */
+    bool constexpr isValid(Code const projectileId, Code const targetId,
+                           HEPEnergyType const sqrtSnn) const;
+
+    /**
+     * Returns inelastic AND elastic cross sections.
+     *
+     * These cross sections must correspond to the process described in doInteraction
+     * AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
+     * nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
+     * Sibyll must be useful inside the NuclearInteraction model, which requires that.
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     *
+     * @return a tuple of: inelastic cross section, elastic cross section
+     */
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code const projectile, Code const target, FourMomentum const& projectileP4,
+        FourMomentum const& targetP4) const;
+
+    /**
+     * Returns inelastic (production) cross section.
+     *
+     * This cross section must correspond to the process described in doInteraction.
+     * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+     *
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     *
+     * @return inelastic cross section
+     * elastic cross section
+     */
+    CrossSectionType getCrossSection(Code const projectile, Code const target,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const {
+      return std::get<0>(
+          getCrossSectionInelEla(projectile, target, projectileP4, targetP4));
+    }
+
+    /**
+     * In this function SIBYLL is called to produce one event. The
+     * event is copied (and boosted) into the shower lab frame.
+     *
+     * @param view is the stack object for the secondaries
+     * @param projectile is the Code of the projectile
+     * @param target is the Code of the target
+     * @param projectileP4: four-momentum of projectile
+     * @param targetP4: four-momentum of target
+     */
+
+    template <typename TSecondaries>
+    void doInteraction(TSecondaries& view, Code const projectile, Code const target,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
+
+  private:
+    void setParticleListStable(std::set<Code>);
+    void setAllParticlesUnstable();
+    HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
+    HEPEnergyType constexpr getMaxEnergyCoM() const { return maxEnergyCoM_; }
+
+    // hard model limits
+    static HEPEnergyType constexpr minEnergyCoM_ = 10. * 1e9 * electronvolt;
+    static HEPEnergyType constexpr maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
+    static unsigned int constexpr maxTargetMassNumber_ = 18;
+    static unsigned int constexpr minNuclearTargetA_ = 4;
+
+    std::shared_ptr<spdlog::logger> logger_ =
+        get_logger("corsika_sibyll_HadronInteractionModel");
+
+    // data members
+    int count_ = 0;
+    int nucCount_ = 0;
+    bool sibyll_listing_;
+    bool const internal_decays_;
+    std::set<Code> stable_particles_;
+  };
+
+} // namespace corsika::sibyll
+
+#include <corsika/detail/modules/sibyll/HadronInteractionModel.inl>

corsika/modules/sibyll/Interaction.hpp

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/core/ParticleProperties.hpp>
-#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
-#include <tuple>
-
-namespace corsika::sibyll {
-
-  class Interaction : public InteractionProcess<Interaction> {
-
-  public:
-    Interaction(bool const sibyll_printout_on = false);
-    ~Interaction();
-
-    bool isValidCoMEnergy(HEPEnergyType const ecm) const {
-      return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
-    }
-    //! sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
-    //! neutrons (p,n == nucleon)
-    bool isValidTarget(Code const TargetId) const {
-      return (is_nucleus(TargetId) && (get_nucleus_A(TargetId) >= minNuclearTargetA_) &&
-              (get_nucleus_A(TargetId) < maxTargetMassNumber_)) ||
-             (TargetId == Code::Proton || TargetId == Code::Hydrogen ||
-              TargetId == Code::Neutron);
-    }
-
-    //! returns production and elastic cross section for hadrons in sibyll. Inputs are:
-    //! CorsikaId of beam particle, CorsikaId of target particle and center-of-mass
-    //! energy. Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
-    std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
-        Code const, Code const, HEPEnergyType const) const;
-
-    template <typename TParticle>
-    GrammageType getInteractionLength(TParticle const&) const;
-
-    /**
-       In this function SIBYLL is called to produce one event. The
-       event is copied (and boosted) into the shower lab frame.
-     */
-
-    template <typename TSecondaries>
-    void doInteraction(TSecondaries&);
-
-  private:
-    void setStable(std::vector<Code> const&);
-    void setUnstable(std::vector<Code> const&);
-
-    void setUnstable(Code const);
-    void setStable(Code const);
-    void setAllUnstable();
-    void setAllStable();
-
-    int getMaxTargetMassNumber() const { return maxTargetMassNumber_; }
-    HEPEnergyType getMinEnergyCoM() const { return minEnergyCoM_; }
-    HEPEnergyType getMaxEnergyCoM() const { return maxEnergyCoM_; }
-
-    default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("sibyll");
-    const bool internalDecays_ = true;
-    const HEPEnergyType minEnergyCoM_ = 10. * 1e9 * electronvolt;
-    const HEPEnergyType maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
-    const int maxTargetMassNumber_ = 18;
-    const int minNuclearTargetA_ = 4;
-
-    // data members
-    int count_ = 0;
-    int nucCount_ = 0;
-    bool sibyll_listing_;
-  };
-
-} // namespace corsika::sibyll
-
-#include <corsika/detail/modules/sibyll/Interaction.inl>

corsika/modules/sibyll/InteractionModel.hpp

+/*
+ * (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 {
+  /**
+   * This class combines sibyll::HadronInteractionModel, which can only handle hadron
+   * projectiles, and sibyll::NuclearInteractionModel, which can handle only nucleus
+   * projectile, into a single process. The getCrossSection() and doInteraction() methods
+   * forward to the underlying models, depending on the projectile type.
+   */
+  class InteractionModel {
+  public:
+    using nuclear_model_type = NuclearInteractionModel<HadronInteractionModel>;
+
+    InteractionModel(std::set<Code> const&, std::set<Code> const&);
+
+    CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+                                     FourMomentum const&) const;
+
+    std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+        Code, Code, FourMomentum const&, FourMomentum const&) const;
+
+    template <typename TSecondaries>
+    void doInteraction(TSecondaries&, Code, Code, FourMomentum const&,
+                       FourMomentum const&);
+
+    HadronInteractionModel& getHadronInteractionModel();
+    HadronInteractionModel const& getHadronInteractionModel() const;
+    nuclear_model_type& getNuclearInteractionModel();
+    nuclear_model_type const& getNuclearInteractionModel() const;
+
+  private:
+    HadronInteractionModel hadronSibyll_;
+    nuclear_model_type nuclearSibyll_;
+  };
+} // namespace corsika::sibyll
+
+#include <corsika/detail/modules/sibyll/InteractionModel.inl>

corsika/modules/sibyll/NuclearInteraction.hpp → corsika/modules/sibyll/NuclearInteractionModel.hpp

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
 
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <set>
 
 namespace corsika::sibyll {
 
-  class Interaction; // fwd-decl
-
   /**
+   * The sibyll::NuclearInteractionModel provides the SIBYLL semi superposition model.
+   *
+   * can transform a proton-nucleus interaction model into a nucleus-nucleus interaction
+   * model.
    *
+   * pass list of nuclei in the environment
    *
-   **/
-  template <class TEnvironment>
-  class NuclearInteraction : public InteractionProcess<NuclearInteraction<TEnvironment>> {
+   * @tparam TNucleonModel
+   */
+  template <class TNucleonModel>
+  class NuclearInteractionModel {
 
   public:
-    NuclearInteraction(sibyll::Interaction&, TEnvironment const&);
-    ~NuclearInteraction();
+    NuclearInteractionModel(TNucleonModel&, std::set<Code> const&);
+
+    ~NuclearInteractionModel();
+
+    bool constexpr isValid(Code const projectileId, Code const targetId,
+                           HEPEnergyType const sqrtSnn) const;
 
-    void initializeNuclearCrossSections();
-    void printCrossSectionTable(Code);
-    CrossSectionType readCrossSectionTable(int const, Code const, HEPEnergyType const);
-    HEPEnergyType getMinEnergyPerNucleonCoM() { return gMinEnergyPerNucleonCoM_; }
-    HEPEnergyType getMaxEnergyPerNucleonCoM() { return gMaxEnergyPerNucleonCoM_; }
-    unsigned int constexpr getMaxNucleusAProjectile() { return gMaxNucleusAProjectile_; }
-    unsigned int constexpr getMaxNFragments() { return gMaxNFragments_; }
-    unsigned int constexpr getNEnergyBins() { return gNEnBins_; }
+    void initializeNuclearCrossSections(std::set<Code> const&);
 
-    template <typename Particle>
-    std::tuple<CrossSectionType, CrossSectionType> getCrossSection(Particle const& p,
-                                                                   const Code TargetId);
+    void printCrossSectionTable(Code) const;
+    CrossSectionType readCrossSectionTable(int const, Code const,
+                                           HEPEnergyType const) const;
+    HEPEnergyType getMinEnergyPerNucleonCoM() const { return gMinEnergyPerNucleonCoM_; }
+    HEPEnergyType getMaxEnergyPerNucleonCoM() const { return gMaxEnergyPerNucleonCoM_; }
+    unsigned int constexpr getMaxNucleusAProjectile() const {
+      return gMaxNucleusAProjectile_;
+    }
+    unsigned int constexpr getMaxNFragments() const { return gMaxNFragments_; }
+    unsigned int constexpr getNEnergyBins() const { return gNEnBins_; }
 
-    template <typename Particle>
-    GrammageType getInteractionLength(Particle const&);
+    CrossSectionType getCrossSection(Code const, Code const,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const;
 
     template <typename TSecondaryView>
-    void doInteraction(TSecondaryView&);
+    void doInteraction(TSecondaryView&, Code const, Code const,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4);
 
   private:
     int count_ = 0;
-    int nucCount_ = 0;
 
-    TEnvironment const& environment_;
-    sibyll::Interaction& hadronicInteraction_;
+    TNucleonModel& hadronicInteraction_;
     std::map<Code, int> targetComponentsIndex_;
-    default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("sibyll");
+    std::shared_ptr<spdlog::logger> logger_ =
+        get_logger("corsika_sibyll_NuclearInteractionModel");
     static unsigned int constexpr gNSample_ =
         500; // number of samples in MC estimation of cross section
     static unsigned int constexpr gMaxNucleusAProjectile_ = 56;
@@ -67,4 +76,4 @@ namespace corsika::sibyll {
 
 } // namespace corsika::sibyll
 
-#include <corsika/detail/modules/sibyll/NuclearInteraction.inl>
+#include <corsika/detail/modules/sibyll/NuclearInteractionModel.inl>

corsika/modules/sibyll/ParticleConversion.hpp

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -24,7 +23,7 @@ namespace corsika::sibyll {
      These are the possible projectile for which Sibyll knows the cross section
    */
   enum class SibyllXSClass : int8_t {
-    CannotInteract = 0,
+    CrossSectionUnknown = 0,
     Baryon = 1,
     Pion = 2,
     Kaon = 3,
@@ -33,14 +32,14 @@ namespace corsika::sibyll {
 
 #include <corsika/modules/sibyll/Generated.inc>
 
-  SibyllCode constexpr convertToSibyll(corsika::Code pCode) {
-    return corsika2sibyll[static_cast<corsika::CodeIntType>(pCode)];
+  SibyllCode constexpr convertToSibyll(Code const pCode) {
+    return corsika2sibyll[static_cast<CodeIntType>(pCode)];
   }
 
-  corsika::Code constexpr convertFromSibyll(SibyllCode pCode) {
+  Code constexpr convertFromSibyll(SibyllCode const pCode) {
     auto const s = static_cast<SibyllCodeIntType>(pCode);
     auto const corsikaCode = sibyll2corsika[s - minSibyll];
-    if (corsikaCode == corsika::Code::Unknown) {
+    if (corsikaCode == Code::Unknown) {
       throw std::runtime_error(std::string("SIBYLL/CORSIKA conversion of ")
                                    .append(std::to_string(s))
                                    .append(" impossible"));
@@ -48,18 +47,26 @@ namespace corsika::sibyll {
     return corsikaCode;
   }
 
-  int constexpr convertToSibyllRaw(corsika::Code pCode) {
-    return static_cast<int>(convertToSibyll(pCode));
+  int constexpr convertToSibyllRaw(Code const code) {
+    return static_cast<int>(convertToSibyll(code));
   }
 
-  int constexpr getSibyllXSCode(corsika::Code pCode) {
+  // find which cross section to use for this particle. maps to either proton, pion or
+  // kaon or unknown
+  int constexpr getSibyllXSCode(Code const code) {
+    if (is_nucleus(code))
+      return static_cast<SibyllXSClassIntType>(SibyllXSClass::CrossSectionUnknown);
     return static_cast<SibyllXSClassIntType>(
-        corsika2sibyllXStype[static_cast<corsika::CodeIntType>(pCode)]);
+        corsika2sibyllXStype[static_cast<CodeIntType>(code)]);
   }
 
-  bool constexpr canInteract(corsika::Code pCode) { return getSibyllXSCode(pCode) > 0; }
+  // find if interaction can be generated.
+  bool constexpr canInteract(Code const pCode) {
+    if (is_nucleus(pCode)) return false;
+    return caninteract[static_cast<CodeIntType>(pCode)];
+  }
 
-  HEPMassType getSibyllMass(corsika::Code const);
+  HEPMassType getSibyllMass(Code const);
 
 } // namespace corsika::sibyll
 

corsika/modules/sibyll/Random.hpp

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
-
-#include <corsika/framework/random/RNGManager.hpp>
-#include <random>
-
-/**
- * \file sibyll/Random.hpp
- *
- * This file is an integral part of the sibyll interface. It must be
- * linked to the executable linked to sibyll exactly once
- *
- */
-
-namespace sibyll {
-
-  double rndm_interface() {
-    static corsika::default_prng_type& rng =
-        corsika::RNGManager::getInstance().getRandomStream("sibyll");
-    std::uniform_real_distribution<double> dist;
-    return dist(rng);
-  }
-
-} // namespace sibyll

corsika/modules/sibyll/SibStack.hpp

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -68,15 +67,14 @@ namespace corsika::sibyll {
     }
   };
 
-  template <typename StackIteratorInterface>
-  class ParticleInterface : public corsika::ParticleBase<StackIteratorInterface> {
+  template <typename TStackIterator>
+  class ParticleInterface : public corsika::ParticleBase<TStackIterator> {
 
-    using corsika::ParticleBase<StackIteratorInterface>::getStackData;
-    using corsika::ParticleBase<StackIteratorInterface>::getIndex;
+    using corsika::ParticleBase<TStackIterator>::getStackData;
+    using corsika::ParticleBase<TStackIterator>::getIndex;
 
   public:
-    void setParticleData(const int vID, // corsika::sibyll::SibyllCode vID,
-                         const HEPEnergyType vE, const MomentumVector& vP,
+    void setParticleData(const int vID, const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType vM) {
       setPID(vID);
       setEnergy(vE);
@@ -84,8 +82,7 @@ namespace corsika::sibyll {
       setMass(vM);
     }
 
-    void setParticleData(ParticleInterface<StackIteratorInterface>& /*parent*/,
-                         const int vID, //  corsika::sibyll::SibyllCode vID,
+    void setParticleData(ParticleInterface<TStackIterator>& /*parent*/, const int vID,
                          const HEPEnergyType vE, const MomentumVector& vP,
                          const HEPMassType vM) {
       setPID(vID);