corsika/detail/framework/stack/CombinedStack.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/stack/SecondaryView.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/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
@@ -15,14 +14,18 @@ namespace corsika {
   template <typename T>
   template <typename... Args>
   ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
-      double const n0, InverseLengthType const lambda, Args&&... args)
+      double const n0, InverseLengthType const lambda, Point const center,
+      LengthType const radius, Args&&... args)
       : T(std::forward<Args>(args)...)
-      , n_0(n0)
-      , lambda_(lambda) {}
+      , n0_(n0)
+      , lambda_(lambda)
+      , center_(center)
+      , radius_(radius) {}
 
   template <typename T>
-  double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
-    return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+  inline double ExponentialRefractiveIndex<T>::getRefractiveIndex(
+      Point const& point) const {
+    return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
   }
 
 } // namespace corsika

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

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

corsika/detail/media/GladstoneDaleRefractiveIndex.inl

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