corsika/detail/framework/stack/Stack.inl

 /*
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
@@ -348,7 +347,7 @@ namespace corsika {
   inline void Stack<StackData, MParticleInterface, MSecondaryProducer>::swap(
       unsigned int const a, unsigned int const b) {
     data_.swap(a, b);
-    std::swap(deleted_[a], deleted_[b]);
+    std::vector<bool>::swap(deleted_[a], deleted_[b]);
   }
 
   template <typename StackData, template <typename> typename MParticleInterface,

corsika/detail/framework/utility/COMBoost.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
@@ -23,7 +22,7 @@ namespace corsika {
                             HEPMassType const massTarget)
       : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
       , rotatedCS_{make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents())} {
-    auto const pProjectile = P4projectile.getSpaceLikeComponents();
+    auto const& pProjectile = P4projectile.getSpaceLikeComponents();
     auto const pProjNormSquared = pProjectile.getSquaredNorm();
     auto const pProjNorm = sqrt(pProjNormSquared);
 
@@ -82,7 +81,7 @@ namespace corsika {
 
   template <typename FourVector>
   inline FourVector COMBoost::toCoM(FourVector const& p4) const {
-    auto pComponents = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
+    auto const pComponents = p4.getSpaceLikeComponents().getComponents(rotatedCS_);
     Eigen::Vector3d eVecRotated = pComponents.getEigenVector();
     Eigen::Vector2d lab;
 
@@ -134,8 +133,10 @@ namespace corsika {
     inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
   }
 
-  inline CoordinateSystemPtr COMBoost::getRotatedCS() const { return rotatedCS_; }
+  inline CoordinateSystemPtr const& COMBoost::getRotatedCS() const { return rotatedCS_; }
 
-  inline CoordinateSystemPtr COMBoost::getOriginalCS() const { return originalCS_; }
+  inline CoordinateSystemPtr const& COMBoost::getOriginalCS() const {
+    return originalCS_;
+  }
 
 } // namespace corsika

corsika/detail/framework/utility/CorsikaData.inl

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

corsika/detail/framework/utility/CorsikaFenvFallback.inl

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

corsika/detail/framework/utility/CorsikaFenvOSX.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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.
  */
 
 /**

corsika/detail/framework/utility/CubicSolver.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
@@ -254,7 +253,7 @@ namespace corsika {
 
     if (pre_opt.size()) {
       x1 = pre_opt[0]; //*std::max_element(pre_opt.begin(), pre_opt.end());
-#ifdef DEBUG
+#ifdef _C8_DEBUG_
       for (long double test_v : pre_opt) {
         CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
                           cubic_function(test_v, a, b, c, d));

corsika/detail/framework/utility/LinearSolver.inl

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

corsika/detail/framework/utility/QuadraticSolver.inl

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

corsika/detail/framework/utility/QuarticSolver.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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
@@ -41,8 +40,8 @@ namespace corsika {
       long double y = x3[0]; // there is always at least one solution
       // The essence - choosing Y with maximal absolute value.
       if (x3.size() == 3) {
-        if (fabs(x3[1]) > fabs(y)) y = x3[1];
-        if (fabs(x3[2]) > fabs(y)) y = x3[2];
+        if (std::abs(x3[1]) > std::abs(y)) y = x3[1];
+        if (std::abs(x3[2]) > std::abs(y)) y = x3[2];
       }
 
       long double q1, q2, p1, p2;
@@ -50,13 +49,13 @@ namespace corsika {
 
       long double Det = y * y - 4 * e;
       CORSIKA_LOG_TRACE("Det={}", Det);
-      if (fabs(Det) < epsilon) // in other words - D==0
+      if (std::abs(Det) < epsilon) // in other words - D==0
       {
         q1 = q2 = y * 0.5;
         // g1+g2 = b && g1+g2 = c-y   <=>   g^2 - b*g + c-y = 0    (p === g)
         Det = b * b - 4 * (c - y);
         CORSIKA_LOG_TRACE("Det={}", Det);
-        if (fabs(Det) < epsilon) { // in other words - D==0
+        if (std::abs(Det) < epsilon) { // in other words - D==0
           p1 = p2 = b * 0.5;
         } else {
           if (Det < 0) return {};
@@ -80,7 +79,7 @@ namespace corsika {
       std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
       std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
       if (quad2.size() > 0) {
-        for (auto val : quad2) quad1.push_back(val);
+        for (auto const val : quad2) quad1.push_back(val);
       }
       return quad1;
     }

corsika/detail/framework/utility/SaveBoostHistogram.inl

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

corsika/detail/media/BaseExponential.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -50,7 +49,7 @@ namespace corsika {
     if (uDotA == 0) {
       return length * rhoStart;
     } else {
-      return rhoStart * (lambda_ / uDotA) * (exp(uDotA * length * invLambda_) - 1);
+      return rhoStart * (lambda_ / uDotA) * expm1(uDotA * length * invLambda_);
     }
   }
 
@@ -66,9 +65,9 @@ namespace corsika {
     if (uDotA == 0) {
       return grammage / rhoStart;
     } else {
-      auto const logArg = grammage * invLambda_ * uDotA / rhoStart + 1;
-      if (logArg > 0) {
-        return lambda_ / uDotA * log(logArg);
+      auto const logArg = grammage * invLambda_ * uDotA / rhoStart;
+      if (logArg > -1) {
+        return lambda_ / uDotA * log1p(logArg);
       } else {
         return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
                meter;

corsika/detail/media/BaseTabular.inl

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

corsika/detail/media/CORSIKA7Atmospheres.inl

 /*
  * (c) Copyright 2021 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
@@ -20,13 +19,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/Environment.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -49,4 +48,24 @@ namespace corsika {
         std::make_unique<TVolumeType>(std::forward<TVolumeArgs>(args)...));
   }
 
+  template <typename IEnvironmentModel>
+  std::set<Code> const get_all_elements_in_universe(
+      Environment<IEnvironmentModel> const& env) {
+
+    auto const& universe = *(env.getUniverse());
+    auto const allElementsInUniverse = std::invoke([&]() {
+      std::set<Code> allElementsInUniverse;
+      auto collectElements = [&](auto& vtn) {
+        if (vtn.hasModelProperties()) {
+          auto const& comp =
+              vtn.getModelProperties().getNuclearComposition().getComponents();
+          for (auto const c : comp) allElementsInUniverse.insert(c);
+        }
+      };
+      universe.walk(collectElements);
+      return allElementsInUniverse;
+    });
+    return allElementsInUniverse;
+  }
+
 } // namespace corsika

corsika/detail/media/ExponentialRefractiveIndex.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the 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,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/FlatExponential.inl

 /*
  * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
  */
 
 #pragma once
@@ -29,7 +28,7 @@ namespace corsika {
   template <typename T>
   inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
     return BaseExponential<FlatExponential<T>>::getMassDensity(
-        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm());
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
   }
 
   template <typename T>

corsika/detail/media/GeomagneticModel.inl

-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
 
+#include <corsika/framework/core/Logging.hpp>
+#include <boost/math/tr1.hpp>
 #include <boost/filesystem.hpp>
+#include <boost/filesystem/fstream.hpp>
 
 #include <stdexcept>
 #include <string>
@@ -10,16 +17,15 @@
 namespace corsika {
 
   inline GeomagneticModel::GeomagneticModel(Point const& center,
-                                                std::string const& dataFile)
+                                            boost::filesystem::path const path)
       : center_(center) {
 
     // Read in coefficients
-    boost::filesystem::path const path = corsika::corsika_data(dataFile);
     boost::filesystem::ifstream file(path, std::ios::in);
 
     // Exit if file opening failed
     if (!file.is_open()) {
-      CORSIKA_LOG_ERROR("Failed opening data file {}", dataFile);
+      CORSIKA_LOG_ERROR("Failed opening data file {}", path);
       throw std::runtime_error("Cannot load GeomagneticModel data.");
     }
 
@@ -73,8 +79,7 @@ namespace corsika {
       int iEpoch = int(epoch);
 
       if (parameters_.count(iEpoch) != 0) {
-        throw std::runtime_error(
-            "GeomagneticModel input file has duplicate Epoch. Fix.");
+        throw std::runtime_error("GeomagneticModel input file has duplicate Epoch. Fix.");
       }
       parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
 
@@ -86,28 +91,30 @@ namespace corsika {
       }
     }
     file.close();
+
+    if (parameters_.size() == 0) {
+      CORSIKA_LOG_ERROR("No input data read!");
+      throw std::runtime_error("No input data read");
+    }
   }
 
   inline MagneticFieldVector GeomagneticModel::getField(double const year,
-                                                          LengthType const altitude,
-                                                          double const latitude,
-                                                          double const longitude) {
+                                                        LengthType const altitude,
+                                                        double const latitude,
+                                                        double const longitude) {
 
     int iYear = int(year);
-    int iEpoch = 0;
-    for (auto parIt = parameters_.rbegin(); parIt != parameters_.rend(); ++parIt) {
-      if (parIt->first <= iYear) {
-        iEpoch = parIt->first;
-        break;
-      }
+    auto iEpoch = parameters_.rbegin();
+    for (; iEpoch != parameters_.rend(); ++iEpoch) {
+      if (iEpoch->first <= iYear) { break; }
     }
-    double epoch = double(iEpoch);
-    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch, year);
-    if (iEpoch == 0) {
+    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch->first, year);
+    if (iEpoch == parameters_.rend()) {
       CORSIKA_LOG_WARN("Year {} is before first EPOCH. Results unclear.", year);
+      iEpoch--; // move one epoch back
     }
-    if (altitude < -1_km || altitude > 850_km) {
-      CORSIKA_LOG_WARN("Altitude should be between -1_km and 850_km.");
+    if (altitude < -1_km || altitude > 600_km) {
+      CORSIKA_LOG_WARN("Altitude should be between -1_km and 600_km.");
     }
     if (latitude <= -90 || latitude >= 90) {
       CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
@@ -119,12 +126,21 @@ namespace corsika {
       CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
     }
 
-    const double lat_geo = latitude * constants::pi / 180;
-    const double lon = longitude * constants::pi / 180;
+    double const epoch = double(iEpoch->first);
+    auto iNextEpoch = iEpoch; // next epoch for interpolation
+    --iNextEpoch;
+    bool const lastEpoch = (iEpoch == parameters_.rbegin());
+    auto const delta_t = year - epoch;
+    CORSIKA_LOG_DEBUG(
+        "identified: t_epoch={}, delta_t={}, lastEpoch={} (false->interpolate)", epoch,
+        delta_t, lastEpoch);
+
+    double const lat_geo = latitude * constants::pi / 180;
+    double const lon = longitude * constants::pi / 180;
 
     // Transform into spherical coordinates
-    const double f = 1 / 298.257223563;
-    const double e_squared = f * (2 - f);
+    double constexpr f = 1 / 298.257223563;
+    double constexpr e_squared = f * (2 - f);
     LengthType R_c =
         constants::EarthRadius::Equatorial / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
     LengthType p = (R_c + altitude) * cos(lat_geo);
@@ -135,16 +151,33 @@ namespace corsika {
     double legendre, next_legendre, derivate_legendre;
     double magneticfield[3] = {0, 0, 0};
 
-    for (size_t j = 0; j < parameters_[iEpoch].size(); j++) {
+    for (size_t j = 0; j < iEpoch->second.size(); j++) {
 
-      ParameterLine p = parameters_[iEpoch][j];
+      ParameterLine p = iEpoch->second[j];
 
       // Time interpolation
-      p.g = p.g + (year - epoch) * p.dg;
-      p.h = p.h + (year - epoch) * p.dh;
+      if (iEpoch == parameters_.rbegin()) {
+        // this is the latest epoch in time, or time-dependence (dg/dh) was specified
+        // we use the extrapolation factors dg/dh:
+        p.g = p.g + delta_t * p.dg;
+        p.h = p.h + delta_t * p.dh;
+      } else {
+        // we linearly interpolate between two epochs
+        ParameterLine const next_p = iNextEpoch->second[j];
+        double const length = iNextEpoch->first - epoch;
+        double p_g = p.g + (next_p.g - p.g) * delta_t / length;
+        double p_h = p.h + (next_p.h - p.h) * delta_t / length;
+        CORSIKA_LOG_TRACE(
+            "interpolation: delta-g={}, delta-h={}, delta-t={}, length={} g1={} g2={} "
+            "g={} h={} ",
+            next_p.g - p.g, next_p.h - p.h, year - epoch, length, next_p.g, p.g, p_g,
+            p_h);
+        p.g = p_g;
+        p.h = p_h;
+      }
 
-      legendre = pow(-1, p.m) * std::assoc_legendre(p.n, p.m, sin(lat_sph));
-      next_legendre = pow(-1, p.m) * std::assoc_legendre(p.n + 1, p.m, sin(lat_sph));
+      legendre = boost::math::tr1::assoc_legendre(p.n, p.m, sin(lat_sph));
+      next_legendre = boost::math::tr1::assoc_legendre(p.n + 1, p.m, sin(lat_sph));
 
       // Schmidt semi-normalization
       if (p.m > 0) {
@@ -180,7 +213,8 @@ namespace corsika {
                            magneticfield[2] * cos(lat_sph - lat_geo);
 
     return MagneticFieldVector{center_.getCoordinateSystem(), magneticfield_geo[0] * 1_nT,
-                               magneticfield_geo[1] * -1_nT, magneticfield_geo[2] * -1_nT};
+                               magneticfield_geo[1] * -1_nT,
+                               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 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  template <typename T>
+  template <typename... Args>
+  inline GladstoneDaleRefractiveIndex<T>::GladstoneDaleRefractiveIndex(
+      double const referenceRefractiveIndex, Point const point, Args&&... args)
+      : T(std::forward<Args>(args)...)
+      , referenceRefractivity_(referenceRefractiveIndex - 1)
+      , referenceInvDensity_(1 / this->getMassDensity(point)) {}
+
+  template <typename T>
+  inline double GladstoneDaleRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return referenceRefractivity_ * (this->getMassDensity(point) * referenceInvDensity_) +
+           1.;
+  }
+
+} // namespace corsika

corsika/detail/media/HomogeneousMedium.inl

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

corsika/detail/media/InhomogeneousMedium.inl

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