README.md

@@ -82,7 +82,7 @@ If you do not have Conan installed, it can be
 installed with:
 
 ``` shell
-pip install --user conan~=1.55.0
+pip install --user conan~=1.57.0
 ```
 
 ### Compiling
@@ -94,7 +94,7 @@ git clone --recursive git@gitlab.iap.kit.edu:AirShowerPhysics/corsika.git
 mkdir corsika-build
 cd corsika-build
 ../corsika/conan-install.sh
-cmake ../corsika -DCMAKE_INSTALL_PREFIX=../corsika-install
+cmake ../corsika -DCMAKE_BUILD_TYPE="RelWithDebInfo" -DCMAKE_INSTALL_PREFIX=../corsika-install
 make -j8
 make install
 ```

README.rst

@@ -86,7 +86,7 @@ If you do not have Conan installed, it can be
 installed with:
 ::
 
-  pip install --user conan~=1.55.0
+  pip install --user conan~=1.57.0
 
   
 Compiling
@@ -99,7 +99,7 @@ Once Conan is installed, follow these steps to download and install CORSIKA 8:
   mkdir corsika-build
   cd corsika-build
   ../corsika/conan-install.sh
-  cmake ../corsika -DCMAKE_INSTALL_PREFIX=../corsika-install
+  cmake ../corsika -DCMAKE_BUILD_TYPE="RelWithDebInfo" -DCMAKE_INSTALL_PREFIX=../corsika-install
   make -j8
   make install
 

corsika/detail/media/CORSIKA7Atmospheres.inl

@@ -20,13 +20,7 @@ namespace corsika {
         TEnvironmentInterface, TExtraEnv>::create(center, constants::EarthRadius::Mean,
                                                   std::forward<TArgs>(args)...);
 
-    // composition values from AIRES manual
-    builder.setNuclearComposition({{
-                                       Code::Nitrogen,
-                                       Code::Argon,
-                                       Code::Oxygen,
-                                   },
-                                   {0.7847, 0.0047, 1. - 0.7847 - 0.0047}});
+    builder.setNuclearComposition(standardAirComposition);
 
     // add the standard atmosphere layers
     auto const params = atmosphereParameterList[static_cast<uint8_t>(atmId)];

corsika/detail/media/ExponentialRefractiveIndex.inl

@@ -24,7 +24,8 @@ namespace corsika {
       , radius_(radius) {}
 
   template <typename T>
-  double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
+  inline double ExponentialRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
     return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
   }
 

corsika/detail/media/GeomagneticModel.inl

@@ -8,6 +8,8 @@
 
 #include <corsika/framework/core/Logging.hpp>
 #include <boost/math/tr1.hpp>
+#include <boost/filesystem.hpp>
+#include <boost/filesystem/fstream.hpp>
 
 #include <stdexcept>
 #include <string>
@@ -216,4 +218,4 @@ namespace corsika {
                                magneticfield_geo[2] * -1_nT};
   }
 
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika

corsika/detail/media/GladstoneDaleRefractiveIndex.inl

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

corsika/detail/media/NuclearComposition.inl

@@ -46,7 +46,7 @@ namespace corsika {
   }
 
   template <typename TFunction>
-  inline auto NuclearComposition::getWeighted(TFunction const& func) const {
+  inline auto NuclearComposition::getWeighted(TFunction func) const {
     using ResultQuantity = decltype(func(std::declval<Code>()));
     auto const product = [&](auto const compID, auto const fraction) {
       return func(compID) * fraction;
@@ -66,7 +66,7 @@ namespace corsika {
   } // namespace corsika
 
   template <typename TFunction>
-  inline auto NuclearComposition::getWeightedSum(TFunction const& func) const
+  inline auto NuclearComposition::getWeightedSum(TFunction func) const
       -> decltype(func(std::declval<Code>())) {
     using ResultQuantity = decltype(func(std::declval<Code>()));
 

corsika/detail/modules/conex/CONEXhybrid.inl

@@ -7,6 +7,7 @@
  */
 
 #include <corsika/framework/core/Logging.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
 #include <corsika/modules/conex/CONEXhybrid.hpp>
 #include <corsika/modules/conex/CONEX_f.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
@@ -17,7 +18,7 @@
 
 #include <algorithm>
 #include <fstream>
-#include <iomanip>
+#include <numeric>
 #include <utility>
 
 namespace corsika {
@@ -83,6 +84,28 @@ namespace corsika {
     CORSIKA_LOG_DEBUG("showerCore (C8): {}",
                       showerCore_.getCoordinates(center.getCoordinateSystem()));
 
+    auto const& components = ::corsika::standardAirComposition.getComponents();
+    auto const& fractions = ::corsika::standardAirComposition.getFractions();
+    if (::corsika::standardAirComposition.getSize() != 3) {
+      throw std::runtime_error{"CONEXhybrid only usable with standard 3-component air"};
+    }
+
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.aira.begin(),
+                   get_nucleus_A);
+    std::transform(components.cbegin(), components.cend(), ::conex::cxair_.airz.begin(),
+                   get_nucleus_Z);
+    std::copy(fractions.cbegin(), fractions.cend(), ::conex::cxair_.airw.begin());
+
+    ::conex::cxair_.airava =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.aira.cbegin(), 0.);
+    ::conex::cxair_.airavz =
+        std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(),
+                           ::conex::cxair_.airz.cbegin(), 0.);
+
+    // this is the CONEX default but actually unused there
+    ::conex::cxair_.airi = {82.0e-09, 95.0e-09, 188.e-09};
+
     int randomSeeds[3] = {1234, 0,
                           0}; // SEEDS ARE NOT USED. All random numbers are obtained from
                               // the CORSIKA 8 stream "conex" and "epos"!

corsika/detail/modules/proposal/InteractionModel.inl

@@ -133,7 +133,10 @@ namespace corsika::proposal {
           this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
         } else {
           auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
-          view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
+          // use mass provided by PROPOSAL to ensure correct conversion to kinetic energy
+          auto massProposal =
+              PROPOSAL::ParticleDef::GetParticleDefForType(s.type).mass * 1_MeV;
+          view.addSecondary(std::make_tuple(sec_code, E - massProposal, dir));
         }
       }
     }

corsika/detail/modules/radio/RadioProcess.inl

@@ -36,13 +36,13 @@ namespace corsika {
                                     TPropagator>::doContinuous(const Step<Particle>& step,
                                                                const bool) {
     // we want the following particles:
-    // Code::Electron & Code::Positron & Code::Gamma
+    // Code::Electron & Code::Positron
 
     // we wrap Simulate() in doContinuous as the plan is to add particle level
     // filtering or thinning for calculation of the radio emission. This is
     // important for controlling the runtime of radio (by ignoring particles
     // that aren't going to contribute i.e. heavy hadrons)
-    // if (valid(particle, track)) {
+    // if (valid(step)) {
     auto const particleID_{step.getParticlePre().getPID()};
     if ((particleID_ == Code::Electron) || (particleID_ == Code::Positron)) {
       CORSIKA_LOG_DEBUG("Particle for radio calculation: {} ", particleID_);
@@ -62,17 +62,28 @@ namespace corsika {
     return meter * std::numeric_limits<double>::infinity();
   }
 
-  // this should all be moved at a separate radio output function
-  // LCOV_EXCL_START
   template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
   inline void RadioProcess<TAntennaCollection, TRadioImpl, TPropagator>::startOfLibrary(
       const boost::filesystem::path& directory) {
 
-    // loop over every antenna and set the initial path
-    // this also writes the time-bins to disk.
-    for (auto& antenna : antennas_.getAntennas()) {
-      antenna.startOfLibrary(directory, this->implementation().algorithm);
-    }
+    // setup the streamer
+    output_.initStreamer((directory / ("antennas.parquet")).string());
+    // LCOV_EXCL_START
+    // build the schema
+    output_.addField("Time", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ex", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ey", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+
+    output_.addField("Ez", parquet::Repetition::REQUIRED, parquet::Type::DOUBLE,
+                     parquet::ConvertedType::NONE);
+    // LCOV_EXCL_STOP
+    // and build the streamer
+    output_.buildStreamer();
   }
 
   template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
@@ -83,13 +94,55 @@ namespace corsika {
     // flush data to disk, and then reset the antenna
     // before the next event
     for (auto& antenna : antennas_.getAntennas()) {
-      antenna.endOfShower(event_, this->implementation().algorithm,
-                          antenna.getSampleRate() * 1_s);
+
+      auto const sampleRate = antenna.getSampleRate() * 1_s;
+      auto const radioImplementation =
+          static_cast<std::string>(this->implementation().algorithm);
+
+      // get the axis labels for this antenna and write the first row.
+      axistype axis = antenna.implementation().getAxis();
+
+      // get the copy of the waveform data for this event
+      std::vector<double> const& dataX = antenna.implementation().getWaveformX();
+      std::vector<double> const& dataY = antenna.implementation().getWaveformY();
+      std::vector<double> const& dataZ = antenna.implementation().getWaveformZ();
+
+      // check for the axis name
+      std::string label = "Unknown";
+      if (antenna.getDomainLabel() == "Time") {
+        label = "Time";
+      }
+      // LCOV_EXCL_START
+      else if (antenna.getDomainLabel() == "Frequency") {
+        label = "Frequency";
+      }
+      // LCOV_EXCL_STOP
+      if (radioImplementation == "ZHS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          auto time = (axis.at(i + 1) + axis.at(i)) / 2.;
+          auto Ex = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
+          auto Ey = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
+          auto Ez = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
+
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(time) << static_cast<double>(Ex)
+              << static_cast<double>(Ey) << static_cast<double>(Ez) << parquet::EndRow;
+        }
+      } else if (radioImplementation == "CoREAS" && label == "Time") {
+        for (size_t i = 0; i < axis.size() - 1; i++) {
+          *(output_.getWriter())
+              << showerId_ << static_cast<double>(axis[i])
+              << static_cast<double>(dataX[i]) << static_cast<double>(dataY[i])
+              << static_cast<double>(dataZ[i]) << parquet::EndRow;
+        }
+      }
+
       antenna.reset();
     }
+    output_.closeStreamer();
 
     // increment our event counter
-    event_++;
+    showerId_++;
   }
 
   template <typename TAntennaCollection, typename TRadioImpl, typename TPropagator>
@@ -123,6 +176,5 @@ namespace corsika {
 
     return config;
   }
-  // LCOV_EXCL_STOP
 
 } // namespace corsika
\ No newline at end of file

corsika/detail/modules/radio/ZHS.inl

@@ -214,7 +214,7 @@ namespace corsika {
               VectorPotential Vp =
                   betaPerp * sign * constants * f / denominator / midPaths[i].R_distance_;
               antenna.receive(detectionTime1, betaPerp, Vp);
-              // intermidiate contributions
+              // intermediate contributions
               for (int it{1}; it < numberOfBins; ++it) {
                 Vp = betaPerp * sign * constants / denominator / midPaths[i].R_distance_;
                 antenna.receive(

corsika/detail/modules/radio/antennas/Antenna.inl

@@ -28,85 +28,6 @@ namespace corsika {
     return name_;
   }
 
-  template <typename TAntennaImpl>
-  inline void Antenna<TAntennaImpl>::startOfLibrary(
-      const boost::filesystem::path& directory, const std::string radioImplementation) {
-
-    // calculate and save our filename
-    filename_ = (directory / this->getName()).string() + ".npz";
-
-    // get the axis labels for this antenna and write the first row.
-    axistype axis = this->implementation().getAxis();
-
-    // check for the axis name
-    std::string label = "Unknown";
-    if constexpr (TAntennaImpl::is_time_domain) {
-      label = "Time";
-    } else if constexpr (TAntennaImpl::is_freq_domain) {
-      label = "Frequency";
-    }
-
-    if (radioImplementation == "ZHS" && TAntennaImpl::is_time_domain) {
-      for (size_t i = 0; i < axis.size() - 1; i++) {
-        axis.at(i) = (axis.at(i + 1) + axis.at(i)) / 2.;
-      }
-      // explicitly convert the arrays to the needed type for cnpy
-      axis.pop_back();
-      long double const* raw_data = axis.data();
-      std::vector<size_t> N = {
-          axis.size()}; // cnpy needs a vector here -- this should be axis.size() - 1 --
-      // write the labels to the first row of the NumPy file
-      cnpy::npz_save(filename_, label, raw_data, N, "w");
-    } else {
-      // explicitly convert the arrays to the needed type for cnpy
-      long double const* raw_data = axis.data();
-      std::vector<size_t> N = {axis.size()}; // cnpy needs a vector here
-      // write the labels to the first row of the NumPy file
-      cnpy::npz_save(filename_, label, raw_data, N, "w");
-    }
-  }
-
-  template <typename TAntennaImpl>
-  inline void Antenna<TAntennaImpl>::endOfShower(const int event,
-                                                 std::string const& radioImplementation,
-                                                 const double sampleRate) {
-
-    // get the copy of the waveform data for this event
-    // we transpose it so that we can match dimensions with the
-    // time array that is already in the output file
-    std::vector<double> const& dataX = this->implementation().getWaveformX();
-    std::vector<double> const& dataY = this->implementation().getWaveformY();
-    std::vector<double> const& dataZ = this->implementation().getWaveformZ();
-
-    if (radioImplementation == "ZHS") {
-      std::vector<double> electricFieldX(dataX.size() - 1,
-                                         0); // num_bins_, std::vector<double>(3, 0)
-      std::vector<double> electricFieldY(dataY.size() - 1, 0);
-      std::vector<double> electricFieldZ(dataZ.size() - 1, 0);
-      for (size_t i = 0; i < electricFieldX.size(); i++) {
-        electricFieldX.at(i) = -(dataX.at(i + 1) - dataX.at(i)) * sampleRate;
-        electricFieldY.at(i) = -(dataY.at(i + 1) - dataY.at(i)) * sampleRate;
-        electricFieldZ.at(i) = -(dataZ.at(i + 1) - dataZ.at(i)) * sampleRate;
-      }
-      // cnpy needs a vector for the shape
-      cnpy::npz_save(filename_, std::to_string(event) + "X", electricFieldX.data(),
-                     {electricFieldX.size()}, "a");
-      cnpy::npz_save(filename_, std::to_string(event) + "Y", electricFieldY.data(),
-                     {electricFieldY.size()}, "a");
-      cnpy::npz_save(filename_, std::to_string(event) + "Z", electricFieldZ.data(),
-                     {electricFieldZ.size()}, "a");
-    } else {
-      // cnpy needs a vector for the shape
-      // and write this event to the .npz archive
-      cnpy::npz_save(filename_, std::to_string(event) + "X", dataX.data(), {dataX.size()},
-                     "a");
-      cnpy::npz_save(filename_, std::to_string(event) + "Y", dataY.data(), {dataY.size()},
-                     "a");
-      cnpy::npz_save(filename_, std::to_string(event) + "Z", dataZ.data(), {dataZ.size()},
-                     "a");
-    }
-  }
-
   template <typename TAntennaImpl>
   inline TAntennaImpl& Antenna<TAntennaImpl>::implementation() {
     return static_cast<TAntennaImpl&>(*this);

corsika/detail/modules/radio/antennas/TimeDomainAntenna.inl

@@ -97,6 +97,8 @@ namespace corsika {
 
   inline auto const& TimeDomainAntenna::getWaveformZ() const { return waveformEZ_; }
 
+  inline std::string const TimeDomainAntenna::getDomainLabel() { return "Time"; }
+
   inline std::vector<long double> TimeDomainAntenna::createTimeAxis() const {
 
     // create a 1-D xtensor to store time values so we can print them later.
@@ -115,7 +117,7 @@ namespace corsika {
     return times;
   }
 
-  inline auto const& TimeDomainAntenna::getAxis() const { return time_axis_; }
+  inline auto const TimeDomainAntenna::getAxis() const { return time_axis_; }
 
   inline InverseTimeType const& TimeDomainAntenna::getSampleRate() const {
     return sample_rate_;

corsika/detail/modules/urqmd/UrQMD.inl

@@ -19,6 +19,7 @@
 #include <corsika/framework/utility/COMBoost.hpp>
 
 #include <boost/filesystem.hpp>
+#include <boost/filesystem/fstream.hpp>
 #include <boost/multi_array.hpp>
 
 #include <algorithm>

corsika/detail/output/YAMLStreamer.inl

@@ -8,6 +8,8 @@
 
 #pragma once
 
+#include <boost/filesystem/fstream.hpp>
+
 namespace corsika {
 
   inline void YAMLStreamer::writeYAML(YAML::Node const& node,

corsika/framework/core/PhysicalUnits.hpp

@@ -92,6 +92,8 @@ namespace corsika::units::si {
       phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
   using InverseTimeType =
       phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
+  using InverseMassDensityType =
+      phys::units::quantity<phys::units::dimensions<3, -1, 0>, double>;
   using InverseGrammageType =
       phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
   using MagneticFluxType =

corsika/framework/process/DynamicInteractionProcess.hpp

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/geometry/FourVector.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/process/InteractionProcess.hpp>
+#include <corsika/framework/process/ProcessTraits.hpp>
+
+#include <boost/type_index.hpp>
+
+#include <memory>
+
+namespace corsika {
+
+  /**
+   * This class allows selecting/using different InteractionProcesses at runtime without
+   * recompiling the process sequence. The implementation is based on the "type-erasure"
+   * technique.
+   *
+   * @tparam TStack the stack type; has to match the one used by Cascade together with
+   * the ProcessSequence containing the DynamicInteractionProcess.
+   */
+  template <typename TStack>
+  class DynamicInteractionProcess
+      : InteractionProcess<DynamicInteractionProcess<TStack>> {
+  public:
+    using stack_view_type = typename TStack::stack_view_type;
+
+    DynamicInteractionProcess() = default;
+
+    /**
+     * Create new DynamicInteractionProcess. Calls to this instance will be forwared
+     * to the process referred to by obj. The newly created DynamicInteractionProcess
+     * shares ownership of the underlying process, so that you do not need to worry
+     * about it going out of scope.
+     */
+    template <typename TInteractionProcess>
+    DynamicInteractionProcess(std::shared_ptr<TInteractionProcess> obj)
+        : concept_{std::make_unique<ConcreteModel<TInteractionProcess>>(std::move(obj))} {
+      CORSIKA_LOG_DEBUG("creating DynamicInteractionProcess from {}",
+                        boost::typeindex::type_id<TInteractionProcess>().pretty_name());
+
+      // poor man's check while we don't have C++20 concepts
+      static_assert(is_interaction_process_v<TInteractionProcess>);
+
+      // since we only forward interaction process API, all other types of corsika
+      // processes are prohibited
+      static_assert(!is_continuous_process_v<TInteractionProcess>);
+      static_assert(!is_decay_process_v<TInteractionProcess>);
+      static_assert(!is_stack_process_v<TInteractionProcess>);
+      static_assert(!is_cascade_equations_process_v<TInteractionProcess>);
+      static_assert(!is_secondaries_process_v<TInteractionProcess>);
+      static_assert(!is_boundary_process_v<TInteractionProcess>);
+      static_assert(!is_cascade_equations_process_v<TInteractionProcess>);
+    }
+
+    /// forwards arguments to doInteraction() of wrapped instance
+    void doInteraction(stack_view_type& view, Code projectileId, Code targetId,
+                       FourMomentum const& projectileP4, FourMomentum const& targetP4) {
+      return concept_->doInteraction(view, projectileId, targetId, projectileP4,
+                                     targetP4);
+    }
+
+    /// forwards arguments to getCrossSection() of wrapped instance
+    CrossSectionType getCrossSection(Code projectileId, Code targetId,
+                                     FourMomentum const& projectileP4,
+                                     FourMomentum const& targetP4) const {
+      return concept_->getCrossSection(projectileId, targetId, projectileP4, targetP4);
+    }
+
+  private:
+    struct IInteractionModel {
+      virtual ~IInteractionModel() = default;
+      virtual void doInteraction(stack_view_type&, Code, Code, FourMomentum const&,
+                                 FourMomentum const&) = 0;
+      virtual CrossSectionType getCrossSection(Code, Code, FourMomentum const&,
+                                               FourMomentum const&) const = 0;
+    };
+
+    template <typename TModel>
+    struct ConcreteModel final : IInteractionModel {
+      ConcreteModel(std::shared_ptr<TModel> obj) noexcept
+          : model_{std::move(obj)} {}
+
+      void doInteraction(stack_view_type& view, Code projectileId, Code targetId,
+                         FourMomentum const& projectileP4,
+                         FourMomentum const& targetP4) override {
+        return model_->doInteraction(view, projectileId, targetId, projectileP4,
+                                     targetP4);
+      }
+
+      CrossSectionType getCrossSection(Code projectileId, Code targetId,
+                                       FourMomentum const& projectileP4,
+                                       FourMomentum const& targetP4) const override {
+        return model_->getCrossSection(projectileId, targetId, projectileP4, targetP4);
+      }
+
+    private:
+      std::shared_ptr<TModel> model_;
+    };
+
+    std::unique_ptr<IInteractionModel> concept_;
+  };
+} // namespace corsika

corsika/media/CORSIKA7Atmospheres.hpp

@@ -11,6 +11,7 @@
 #include <corsika/media/IRefractiveIndexModel.hpp>
 #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
 #include <corsika/framework/utility/ImplementsMixin.hpp>
+#include <corsika/media/NuclearComposition.hpp>
 
 // for detail namespace, NoExtraModelInner, NoExtraModel and traits
 #include <corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -204,6 +205,10 @@ namespace corsika {
   void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
                                 Point const& center, TArgs... args);
 
+  //! The standard/default air composition with fraction values based on CORSIKA 7
+  static inline NuclearComposition const standardAirComposition{
+      {Code::Nitrogen, Code::Oxygen, Code::Argon}, {0.78479, .21052, 0.00469}};
+
 } // namespace corsika
 
-#include <corsika/detail/media/CORSIKA7Atmospheres.inl>
\ No newline at end of file
+#include <corsika/detail/media/CORSIKA7Atmospheres.inl>

corsika/media/GladstoneDaleRefractiveIndex.hpp

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

corsika/media/NuclearComposition.hpp

@@ -51,7 +51,7 @@ namespace corsika {
      *  @retval returns the vector with weighted return types of func.
      */
     template <typename TFunction>
-    auto getWeighted(TFunction const& func) const;
+    auto getWeighted(TFunction func) const;
 
     /**
      * Sum all all relative composition weighted by func(element)
@@ -65,8 +65,7 @@ namespace corsika {
      *  @retval returns the weighted sum with the type defined by the return type of func.
      */
     template <typename TFunction>
-    auto getWeightedSum(TFunction const& func) const
-        -> decltype(func(std::declval<Code>()));
+    auto getWeightedSum(TFunction func) const -> decltype(func(std::declval<Code>()));
 
     /**
      * Number of elements in the composition array