updates and fixes

corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl

@@ -13,7 +13,7 @@
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
 #include <corsika/media/SlidingPlanarExponential.hpp>
-//#include <corsika/media/SlidingPlanarTabular.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
 
 namespace corsika {
 
@@ -42,7 +42,7 @@ namespace corsika {
       TModelArgs...>::addExponentialLayer(GrammageType b, LengthType c,
                                           LengthType upperBoundary) {
 
-    auto const radius = earthRadius_ + upperBoundary;
+    auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
 
@@ -56,7 +56,7 @@ namespace corsika {
       auto lastBound = [&](auto... argPack) {
         return std::make_shared<
             TMediumModelExtra<SlidingPlanarExponential<TMediumInterface>>>(
-            argPack..., center_, rho0, -c, *composition_, earthRadius_);
+            argPack..., center_, rho0, -c, *composition_, planetRadius_);
       };
 
       // now unpack the additional arguments
@@ -64,7 +64,7 @@ namespace corsika {
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
-          center_, rho0, -c, *composition_, earthRadius_);
+          center_, rho0, -c, *composition_, planetRadius_);
     }
 
     layers_.push(std::move(node));
@@ -75,14 +75,14 @@ namespace corsika {
   inline void LayeredSphericalAtmosphereBuilder<
       TMediumInterface, TMediumModelExtra,
       TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
-    auto const radius = earthRadius_ + upperBoundary;
+    auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
 
     auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
         std::make_unique<Sphere>(center_, radius));
 
-    units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
+    units::si::GrammageType constexpr b = 1_g / (1_cm * 1_cm);
     auto const rho0 = b / c;
 
     if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
@@ -94,7 +94,6 @@ namespace corsika {
 
       // now unpack the additional arguments
       auto model = std::apply(lastBound, additionalModelArgs_);
-
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
@@ -108,33 +107,33 @@ namespace corsika {
             typename... TModelArgs>
   inline void
   LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
-      addTabularLayer(GrammageType b,
-                      std::function<MassDensityType(LengthType)> const& funcRho,
-                      LengthType upperBoundary) {
+      addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+                      unsigned int const nBins, LengthType const deltaHeight,
+                      LengthType const upperBoundary) {
 
-    auto const radius = earthRadius_ + upperBoundary;
+    auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
 
     auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
         std::make_unique<Sphere>(center_, radius));
-    /*
-        if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
-          // helper lambda in which the last 5 arguments to make_shared<...> are bound
-          auto lastBound = [&](auto... argPack) {
-            return std::make_shared<
-                TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
-                argPack..., center_, rho0, -c, *composition_, earthRadius_);
-          };
-
-          // now unpack the additional arguments
-          auto model = std::apply(lastBound, additionalModelArgs_);
-          node->setModelProperties(std::move(model));
-        } else {
-          node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
-              center_, rho0, -c, *composition_, earthRadius_);
-        }
-    */
+
+    if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+      // helper lambda in which the last 5 arguments to make_shared<...> are bound
+      auto lastBound = [&](auto... argPack) {
+        return std::make_shared<
+            TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
+            argPack..., center_, funcRho, nBins, deltaHeight, *composition_,
+            planetRadius_);
+      };
+
+      // now unpack the additional arguments
+      auto model = std::apply(lastBound, additionalModelArgs_);
+      node->setModelProperties(std::move(model));
+    } else {
+      node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
+          center_, funcRho, nBins, deltaHeight, *composition_, planetRadius_);
+    }
     layers_.push(std::move(node));
   }
 
@@ -167,10 +166,10 @@ namespace corsika {
   template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
   struct make_layered_spherical_atmosphere_builder {
     template <typename... TArgs>
-    static auto create(Point const& center, LengthType earthRadius, TArgs... args) {
+    static auto create(Point const& center, LengthType planetRadius, TArgs... args) {
       return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
                                                TArgs...>{std::forward<TArgs>(args)...,
-                                                         center, earthRadius};
+                                                         center, planetRadius};
     }
   };
 

corsika/detail/media/SlidingPlanarExponential.inl

@@ -8,50 +8,51 @@
 
 #pragma once
 
-#include <corsika/media/SlidingPlanarExponential.hpp>
-
 namespace corsika {
 
-  template <typename T>
-  inline SlidingPlanarExponential<T>::SlidingPlanarExponential(
+  template <typename TDerived>
+  inline SlidingPlanarExponential<TDerived>::SlidingPlanarExponential(
       Point const& p0, MassDensityType rho0, LengthType lambda,
       NuclearComposition const& nuclComp, LengthType referenceHeight)
-      : BaseExponential<SlidingPlanarExponential<T>>(p0, referenceHeight, rho0, lambda)
-      , nuclComp_(nuclComp)
-      , referenceHeight_(referenceHeight) {}
+      : BaseExponential<SlidingPlanarExponential<TDerived>>(p0, referenceHeight, rho0,
+                                                            lambda)
+      , nuclComp_(nuclComp) {}
 
-  template <typename T>
-  inline MassDensityType SlidingPlanarExponential<T>::getMassDensity(
+  template <typename TDerived>
+  inline MassDensityType SlidingPlanarExponential<TDerived>::getMassDensity(
       Point const& point) const {
     auto const heightFull =
-        (point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
+        (point - BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
             .getNorm();
-    return BaseExponential<SlidingPlanarExponential<T>>::getMassDensity(heightFull);
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getMassDensity(
+        heightFull);
   }
 
-  template <typename T>
-  inline NuclearComposition const& SlidingPlanarExponential<T>::getNuclearComposition()
-      const {
+  template <typename TDerived>
+  inline NuclearComposition const&
+  SlidingPlanarExponential<TDerived>::getNuclearComposition() const {
     return nuclComp_;
   }
 
-  template <typename T>
-  inline GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
+  template <typename TDerived>
+  inline GrammageType SlidingPlanarExponential<TDerived>::getIntegratedGrammage(
       BaseTrajectory const& traj) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj,
-                                                                               axis);
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getIntegratedGrammage(
+        traj, axis);
   }
 
-  template <typename T>
-  inline LengthType SlidingPlanarExponential<T>::getArclengthFromGrammage(
+  template <typename TDerived>
+  inline LengthType SlidingPlanarExponential<TDerived>::getArclengthFromGrammage(
       BaseTrajectory const& traj, GrammageType const grammage) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getArclengthFromGrammage(
+    auto const axis =
+        (traj.getPosition(0) -
+         BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .normalized();
+    return BaseExponential<SlidingPlanarExponential<TDerived>>::getArclengthFromGrammage(
         traj, grammage, axis);
   }
 

corsika/media/BaseExponential.hpp

@@ -29,8 +29,6 @@ namespace corsika {
                     MassDensityType rho0, LengthType lambda);
 
     Point const& getAnchorPoint() const { return point_; }
-    MassDensityType getRho0() const { return rho0_; }
-    InverseLengthType getInvLambda() const { return invLambda_; }
 
   protected:
     auto const& getImplementation() const;

corsika/media/LayeredSphericalAtmosphereBuilder.hpp

@@ -22,6 +22,7 @@
 #include <stack>
 #include <tuple>
 #include <type_traits>
+#include <functional>
 
 namespace corsika {
 
@@ -68,9 +69,9 @@ namespace corsika {
 
   protected:
     LayeredSphericalAtmosphereBuilder(TModelArgs... args, Point const& center,
-                                      LengthType earthRadius)
+                                      LengthType planetRadius)
         : center_(center)
-        , earthRadius_(earthRadius)
+        , planetRadius_(planetRadius)
         , additionalModelArgs_{args...} {}
 
   public:
@@ -78,9 +79,9 @@ namespace corsika {
     void addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary);
     void addLinearLayer(LengthType c, LengthType upperBoundary);
 
-    void addTabularLayer(GrammageType b,
-                         std::function<MassDensityType(LengthType)> const& funcRho,
-                         LengthType upperBoundary);
+    void addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+                         unsigned int const nBins, LengthType const deltaHeight,
+                         LengthType const upperBoundary);
 
     int getSize() const { return layers_.size(); }
 
@@ -88,9 +89,9 @@ namespace corsika {
     Environment<TMediumInterface> assemble();
 
     /**
-     * Get the current Earth radius.
+     * Get the current planet radius.
      */
-    LengthType getEarthRadius() const { return earthRadius_; }
+    LengthType getPlanetRadius() const { return planetRadius_; }
 
   private:
     void checkRadius(LengthType r) const;
@@ -98,7 +99,7 @@ namespace corsika {
     std::unique_ptr<NuclearComposition> composition_;
     Point center_;
     LengthType previousRadius_{LengthType::zero()};
-    LengthType earthRadius_;
+    LengthType planetRadius_;
     std::tuple<TModelArgs...> const additionalModelArgs_;
 
     std::stack<typename VolumeTreeNode<TMediumInterface>::VTNUPtr>

corsika/media/SlidingPlanarExponential.hpp

@@ -13,7 +13,6 @@
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/framework/geometry/BaseTrajectory.hpp>
 
@@ -53,7 +52,6 @@ namespace corsika {
 
   private:
     NuclearComposition const nuclComp_;
-    LengthType const referenceHeight_;
   };
 
 } // namespace corsika

examples/CMakeLists.txt

@@ -60,4 +60,6 @@ CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
 add_executable (corsika corsika.cpp)
 target_link_libraries (corsika CORSIKA8::CORSIKA8)
 
+add_executable (mars mars.cpp)
+target_link_libraries (mars CORSIKA8::CORSIKA8)
 

examples/corsika.cpp

@@ -103,8 +103,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
 
 int main(int argc, char** argv) {
 
-  logging::set_level(logging::level::info);
-
   // the main command line description
   CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};
 
@@ -154,10 +152,31 @@ int main(int argc, char** argv) {
       ->group("Misc.");
   app.add_flag("--force-interaction", "Force the location of the first interaction.")
       ->group("Misc.");
+  app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.")
+      ->default_val("info")
+      ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
+      ->group("Misc.");
 
   // parse the command line options into the variables
   CLI11_PARSE(app, argc, argv);
 
+  if (app.count("--verbosity")) {
+    string const loglevel = app["verbosity"]->as<string>();
+    if (loglevel == "warn") {
+      logging::set_level(logging::level::warn);
+    } else if (loglevel == "info") {
+      logging::set_level(logging::level::info);
+    } else if (loglevel == "debug") {
+      logging::set_level(logging::level::debug);
+    } else if (loglevel == "trace") {
+#ifndef DEBUG
+      CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
+      return 1;
+#endif
+      logging::set_level(logging::level::trace);
+    }
+  }
+
   // check that we got either PDG or A/Z
   // this can be done with option_groups but the ordering
   // gets all messed up
@@ -185,17 +204,29 @@ int main(int argc, char** argv) {
                                                                            50_uT, 0_T});
   builder.setNuclearComposition(
       {{Code::Nitrogen, Code::Oxygen},
-       {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+       {0.7847, 1. - 0.7847}}); // values taken from AIRES manual, Ar removed for now
 
   builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
   builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
   builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
   builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
   builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
-  builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
+  builder.addLinearLayer(1e9_cm, 112.8_km);
   builder.assemble(env);
   /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
 
+  ofstream atmout("earth.dat");
+  for (LengthType h = 0_m; h < 110_km; h += 100_m) {
+    Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h};
+    auto rho =
+        env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity(
+            ptest);
+    atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n";
+  }
+  atmout.close();
+
+  /* === START: CONSTRUCT PRIMARY PARTICLE === */
+
   /* === START: CONSTRUCT PRIMARY PARTICLE === */
 
   // parse the primary ID as a PDG or A/Z code
@@ -231,8 +262,8 @@ int main(int argc, char** argv) {
   /* === END: CONSTRUCT PRIMARY PARTICLE === */
 
   /* === START: CONSTRUCT GEOMETRY === */
-  auto const observationHeight = 0_km + builder.getEarthRadius();
-  auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+  auto const observationHeight = 0_km + builder.getPlanetRadius();
+  auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
                       static_pow<2>(injectionHeight));
@@ -281,9 +312,9 @@ int main(int argc, char** argv) {
       Code::KStar0_1430_MinusBar,
   }};
 
-  decaySibyll.printDecayConfig();
+  // decaySibyll.printDecayConfig();
 
-  ParticleCut cut{50_GeV, 50_GeV, 50_GeV, 50_GeV, false};
+  ParticleCut cut{1_GeV, 1_GeV, 1_GeV, 1_GeV, false};
   corsika::proposal::Interaction emCascade(env);
   corsika::proposal::ContinuousProcess emContinuous(env);
   InteractionCounter emCascadeCounted(emCascade);

examples/em_shower.cpp

@@ -128,8 +128,8 @@ int main(int argc, char** argv) {
   cout << "input momentum: " << plab.getComponents() / 1_GeV
        << ", norm = " << plab.getNorm() << endl;
 
-  auto const observationHeight = 1.4_km + builder.getEarthRadius();
-  auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+  auto const observationHeight = 1.4_km + builder.getPlanetRadius();
+  auto const injectionHeight = 112.75_km + builder.getPlanetRadius();
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
                       static_pow<2>(injectionHeight));

examples/hybrid_MC.cpp

@@ -178,8 +178,8 @@ int main(int argc, char** argv) {
   cout << "input momentum: " << plab.getComponents() / 1_GeV
        << ", norm = " << plab.getNorm() << endl;
 
-  auto const observationHeight = 0_km + builder.getEarthRadius();
-  auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+  auto const observationHeight = 0_km + builder.getPlanetRadius();
+  auto const injectionHeight = 112.75_km + builder.getPlanetRadius();
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
                       static_pow<2>(injectionHeight));

examples/vertical_EAS.cpp

@@ -177,8 +177,8 @@ int main(int argc, char** argv) {
   cout << "input momentum: " << plab.getComponents() / 1_GeV
        << ", norm = " << plab.getNorm() << endl;
 
-  auto const observationHeight = 0_km + builder.getEarthRadius();
-  auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+  auto const observationHeight = 0_km + builder.getPlanetRadius();
+  auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
   auto const t = -observationHeight * cos(thetaRad) +
                  sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
                       static_pow<2>(injectionHeight));

tests/media/testEnvironment.cpp

@@ -11,6 +11,7 @@
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/RootCoordinateSystem.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
+
 #include <corsika/media/DensityFunction.hpp>
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
@@ -26,6 +27,7 @@
 #include <corsika/media/LinearApproximationIntegrator.hpp>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
 #include <corsika/media/VolumeTreeNode.hpp>
 
 #include <SetupTestTrajectory.hpp>
@@ -94,7 +96,7 @@ TEST_CASE("FlatExponential") {
 
   Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
   LengthType const lambda = 3_m;
-  auto const rho0 = 1_g / static_pow<3>(1_cm);
+  auto const rho0 = 1_g / cube(1_cm);
   FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
                                              protonComposition);
   SpeedType const speed = 20_m / second;
@@ -188,6 +190,256 @@ TEST_CASE("SlidingPlanarExponential") {
   }
 }
 
+struct RhoFuncConst {
+  MassDensityType operator()(LengthType) const { return 1_g / cube(1_cm); }
+  static GrammageType integrate(LengthType dL) { return dL * 1_g / cube(1_cm); }
+};
+
+struct RhoFuncExp {
+  MassDensityType operator()(LengthType height) const {
+    return 1_g / cube(1_cm) * exp(-height / 1000_m);
+  }
+  static GrammageType integrate(BaseTrajectory const& traj, Point const& origin,
+                                LengthType const& refH) {
+    LengthType height1 = (traj.getPosition(0) - origin).getNorm() - refH;
+    LengthType height2 = (traj.getPosition(1) - origin).getNorm() - refH;
+    if (height1 > height2) { std::swap(height1, height2); }
+
+    DirectionVector const axis(
+        (traj.getPosition(0) - origin).normalized()); // to gravity center
+    double const cosTheta = axis.dot(traj.getDirection(0));
+
+    CORSIKA_LOG_INFO("h1={} h2={} cT={} rho1={}, rho2={}", height1, height2, cosTheta,
+                     1_g / cube(1_cm) * exp(-height1 / 1000_m),
+                     1_g / cube(1_cm) * exp(-height2 / 1000_m));
+    return (1_km * 1_g / cube(1_cm) * exp(-height1 / 1000_m) -
+            1_km * 1_g / cube(1_cm) * exp(-height2 / 1000_m)) /
+           cosTheta;
+  }
+  static GrammageType integrate(BaseTrajectory const& traj, LengthType const& length,
+                                Point const& origin, LengthType const& refH) {
+    LengthType height1 = (traj.getPosition(0) - origin).getNorm() - refH;
+    LengthType height2 =
+        (traj.getPosition(0) + traj.getDirection(0) * length - origin).getNorm() - refH;
+    if (height1 > height2) { std::swap(height1, height2); }
+
+    DirectionVector const axis(
+        (traj.getPosition(0) - origin).normalized()); // to gravity center
+    double const cosTheta = axis.dot(traj.getDirection(0));
+
+    CORSIKA_LOG_INFO("h1={} h2={} cT={}", height1, height2, cosTheta);
+    return (1_km * 1_g / cube(1_cm) * exp(-height1 / 1000_m) -
+            1_km * 1_g / cube(1_cm) * exp(-height2 / 1000_m)) /
+           cosTheta;
+  }
+};
+
+TEST_CASE("SlidingPlanarTabular") {
+
+  logging::set_level(logging::level::info);
+
+  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+                                             std::vector<float>{1.f});
+
+  RhoFuncConst rhoFunc;
+  SlidingPlanarTabular<IMediumModel> const medium(gOrigin, rhoFunc, 1000, 10_m,
+                                                  protonComposition);
+
+  SECTION("density") {
+    CHECK(medium.getMassDensity({gCS, {0_m, 0_m, 3_m}}) /
+              medium.getMassDensity({gCS, {0_m, 3_m, 0_m}}) ==
+          Approx(1));
+    CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}) == 1_g / cube(1_cm));
+    CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 300_m}}) == 1_g / cube(1_cm));
+  }
+
+  SECTION("vertical") {
+    SpeedType const speed = 5_m / second;
+    TimeType const tEnd1 = 1_s;
+    LengthType const length1 = speed * tEnd1;
+    TimeType const tEnd2 = 300_s;
+    LengthType const length2 = speed * tEnd2;
+    Line const line(
+        {gCS, {0_m, 0_m, 1_m}},
+        Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, speed}));
+    setup::Trajectory const trajectory1 =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd1);
+    Line const line1Reverse(
+        trajectory1.getPosition(1),
+        Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, -speed}));
+    setup::Trajectory const trajectory1Reverse =
+        setup::testing::make_track<setup::Trajectory>(line1Reverse, tEnd1);
+
+    setup::Trajectory const trajectory2 =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd2);
+    Line const line2Reverse(
+        trajectory2.getPosition(0),
+        Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, -speed}));
+    setup::Trajectory const trajectory2Reverse =
+        setup::testing::make_track<setup::Trajectory>(line2Reverse, tEnd2);
+
+    // failures
+    CHECK_THROWS(medium.getArclengthFromGrammage(trajectory1, -1_kg / square(1_cm)));
+
+    MassDensityType const rho0 = 1_g / cube(1_cm);
+
+    // short track
+    CHECK(medium.getIntegratedGrammage(trajectory1) == length1 * rho0);
+    LengthType const testD1 = length1 / 200; // within bin
+    CHECK(medium.getArclengthFromGrammage(trajectory1, rho0 * testD1) / testD1 ==
+          Approx(1));
+    // short track, reverse
+    CHECK(medium.getIntegratedGrammage(trajectory1Reverse) == length1 * rho0);
+    CHECK(medium.getArclengthFromGrammage(trajectory1Reverse, rho0 * testD1) / testD1 ==
+          Approx(1));
+
+    // long track
+    CHECK(medium.getIntegratedGrammage(trajectory2) == length2 * 1_g / cube(1_cm));
+    LengthType const testD2 = length2 / 25; // multi bin
+    CHECK(medium.getArclengthFromGrammage(trajectory2, rho0 * testD2) == testD2);
+  }
+
+  SECTION("inclined") {
+    SpeedType const speed = 5_m / second;
+    TimeType const tEnd1 = 1_s;
+    LengthType const length1 = speed * tEnd1;
+    TimeType const tEnd2 = 300_s;
+    LengthType const length2 = speed * tEnd2;
+    Line const line({gCS, {0_m, 0_m, 1_m}},
+                    Vector<SpeedType::dimension_type>(
+                        gCS, {speed / sqrt(2.), 0_m / second, speed / sqrt(2.)}));
+    setup::Trajectory const trajectory1 =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd1);
+    Line const line1Reverse(
+        trajectory1.getPosition(1),
+        Vector<SpeedType::dimension_type>(
+            gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+    setup::Trajectory const trajectory1Reverse =
+        setup::testing::make_track<setup::Trajectory>(line1Reverse, tEnd1);
+
+    setup::Trajectory const trajectory2 =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd2);
+    Line const line2Reverse(
+        trajectory2.getPosition(1),
+        Vector<SpeedType::dimension_type>(
+            gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+    setup::Trajectory const trajectory2Reverse =
+        setup::testing::make_track<setup::Trajectory>(line2Reverse, tEnd2);
+
+    MassDensityType const rho0 = 1_g / cube(1_cm);
+
+    // short track
+    CHECK(medium.getIntegratedGrammage(trajectory1) / (length1 * rho0) == Approx(1));
+    LengthType const testD1 = length1 / 200; // within bin
+    CHECK(medium.getArclengthFromGrammage(trajectory1, RhoFuncConst::integrate(testD1)) /
+              testD1 ==
+          Approx(1));
+    // short track, reverse
+    CHECK(medium.getIntegratedGrammage(trajectory1Reverse) / (length1 * rho0) ==
+          Approx(1));
+    CHECK(medium.getArclengthFromGrammage(trajectory1Reverse, rho0 * testD1) / testD1 ==
+          Approx(1));
+
+    // long track
+    CHECK(medium.getIntegratedGrammage(trajectory2) / (length2 * rho0) ==
+          Approx(1).epsilon(0.01));
+    LengthType const testD2 = length2 / 25; // multi bin
+    CHECK(medium.getArclengthFromGrammage(trajectory2, rho0 * testD2) / testD2 ==
+          Approx(1).epsilon(0.01));
+    // long track reverse
+    CORSIKA_LOG_INFO("length2={}", length2);
+    CHECK(medium.getIntegratedGrammage(trajectory2Reverse) / (length2 * rho0) ==
+          Approx(1).epsilon(0.01));
+    CHECK(medium.getArclengthFromGrammage(trajectory2Reverse, rho0 * testD2) / testD2 ==
+          Approx(1).epsilon(0.01));
+  }
+
+  /*The exponential test is taken over phase-space where the exponential is not so steep
+   * and is samples in sufficient substeps. An reference-height offset of 1000_km is used.
+   * Thus, density is given from 1000 to 1010 km. And curvature effects are small.
+   */
+
+  RhoFuncExp rhoFuncExp;
+  SlidingPlanarTabular<IMediumModel> const mediumExp(gOrigin, rhoFuncExp, 1000, 10_m,
+                                                     protonComposition, 1000_km);
+
+  SECTION("exponential") {
+
+    SpeedType const speed = 5_m / second;
+    TimeType const tEnd1 = 1_s;
+    LengthType const length1 = speed * tEnd1;
+    TimeType const tEnd2 = 300_s;
+    LengthType const length2 = speed * tEnd2;
+    Line const line({gCS, {0_m, 0_m, 1000.005_km}},
+                    Vector<SpeedType::dimension_type>(
+                        gCS, {speed / sqrt(2.), 0_m / second, speed / sqrt(2.)}));
+    setup::Trajectory const trajectory1 =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd1);
+    Line const line1Reverse(
+        trajectory1.getPosition(1),
+        Vector<SpeedType::dimension_type>(
+            gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+    setup::Trajectory const trajectory1Reverse =
+        setup::testing::make_track<setup::Trajectory>(line1Reverse, tEnd1);
+
+    setup::Trajectory const trajectory2 =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd2);
+
+    CORSIKA_LOG_INFO("{} {}", RhoFuncExp::integrate(trajectory1, gOrigin, 1000_km),
+                     length1);
+
+    // short track
+    GrammageType const testShortX = RhoFuncExp::integrate(trajectory1, gOrigin, 1000_km);
+    CHECK(mediumExp.getIntegratedGrammage(trajectory1) / testShortX ==
+          Approx(1).epsilon(0.01));
+    LengthType const testD1 = length1 / 200; // within bin
+    GrammageType const testD1X =
+        RhoFuncExp::integrate(trajectory1, testD1, gOrigin, 1000_km);
+    CHECK(mediumExp.getArclengthFromGrammage(trajectory1, testD1X) / testD1 ==
+          Approx(1).epsilon(0.01));
+    // short track, reverse
+    CHECK(mediumExp.getIntegratedGrammage(trajectory1Reverse) / testShortX ==
+          Approx(1).epsilon(0.01));
+    CHECK(mediumExp.getArclengthFromGrammage(trajectory1Reverse, testD1X) / testD1 ==
+          Approx(1).epsilon(0.01));
+
+    // long track
+    GrammageType const testLongX = RhoFuncExp::integrate(trajectory2, gOrigin, 1000_km);
+    CORSIKA_LOG_INFO("testLongX={}", testLongX);
+    CHECK(mediumExp.getIntegratedGrammage(trajectory2) / testLongX ==
+          Approx(1).epsilon(0.01));
+    LengthType const testD2 = length2 / 25; // multi bin
+    GrammageType const testD2X =
+        RhoFuncExp::integrate(trajectory2, testD2, gOrigin, 1000_km);
+    CHECK(mediumExp.getArclengthFromGrammage(trajectory2, testD2X) / testD2 ==
+          Approx(1).epsilon(0.01));
+    // long track, reverse
+
+    // first full trajectory2 reverse
+    Line line2Reverse(trajectory2.getPosition(1),
+                      Vector<SpeedType::dimension_type>(
+                          gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+    setup::Trajectory trajectory2Reverse =
+        setup::testing::make_track<setup::Trajectory>(line2Reverse, tEnd2);
+
+    CHECK(mediumExp.getIntegratedGrammage(trajectory2Reverse) / testLongX ==
+          Approx(1).epsilon(0.01));
+
+    // but now shorter trajectory2 reversed to correspond 100% to testD2
+
+    line2Reverse = Line(trajectory2.getPosition(0) + trajectory2.getDirection(0) * testD2,
+                        Vector<SpeedType::dimension_type>(
+                            gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+    auto const trajectory2ReverseShort =
+        setup::testing::make_track<setup::Trajectory>(line2Reverse, testD2 / speed);
+
+    CORSIKA_LOG_INFO("here {} {} {}", trajectory2ReverseShort.getLength(), testD2,
+                     testD2X / 1_g * square(1_cm));
+    CHECK(mediumExp.getArclengthFromGrammage(trajectory2ReverseShort, testD2X) / testD2 ==
+          Approx(1).epsilon(0.01));
+  }
+}
+
 MassDensityType constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
 
 struct ExponentialTest {
@@ -289,7 +541,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
 
   CHECK(builder.getSize() == 0);
 
-  auto const R = builder.getEarthRadius();
+  auto const R = builder.getPlanetRadius();
 
   CHECK(univ->getChildNodes().size() == 1);
 
@@ -334,7 +586,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
 
   CHECK(builder.getSize() == 0);
   CHECK(univ->getChildNodes().size() == 1);
-  auto const R = builder.getEarthRadius();
+  auto const R = builder.getPlanetRadius();
 
   // check magnetic field at several locations
   const Point pTest(gCS, -10_m, 4_m, R + 35_m);