diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl
index 3d7790c29bf80479a275e8a7349953d31c65c557..dbd4180410675b3c61c37438737ccb732db30768 100644
--- a/corsika/detail/media/BaseExponential.inl
+++ b/corsika/detail/media/BaseExponential.inl
@@ -14,23 +14,43 @@
 
 namespace corsika {
 
+  template <typename TDerived>
+  inline BaseExponential<TDerived>::BaseExponential(Point const& point,
+                                                    LengthType const referenceHeight,
+                                                    MassDensityType rho0,
+                                                    LengthType lambda)
+      : rho0_(rho0)
+      , lambda_(lambda)
+      , invLambda_(1 / lambda)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {}
+
   template <typename TDerived>
   inline auto const& BaseExponential<TDerived>::getImplementation() const {
     return *static_cast<TDerived const*>(this);
   }
 
+  template <typename TDerived>
+  inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+      LengthType const height) const {
+    return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+  }
+
   template <typename TDerived>
   inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
-      BaseTrajectory const& traj, LengthType vL, DirectionVector const& axis) const {
-    if (vL == LengthType::zero()) { return GrammageType::zero(); }
+      BaseTrajectory const& traj, DirectionVector const& axis) const {
+    LengthType const length = traj.getLength();
+    if (length == LengthType::zero()) { return GrammageType::zero(); }
 
-    auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
-    auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis).magnitude();
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
 
     if (uDotA == 0) {
-      return vL * rhoStart;
+      return length * rhoStart;
     } else {
-      return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1);
+      return rhoStart * (lambda_ / uDotA) * (exp(uDotA * length * invLambda_) - 1);
     }
   }
 
@@ -38,8 +58,10 @@ namespace corsika {
   inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
       BaseTrajectory const& traj, GrammageType grammage,
       DirectionVector const& axis) const {
-    auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
-    auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+    // this corresponds to height:
+    double const uDotA = traj.getDirection(0).dot(axis).magnitude();
+    MassDensityType const rhoStart =
+        getImplementation().getMassDensity(traj.getPosition(0));
 
     if (uDotA == 0) {
       return grammage / rhoStart;
@@ -54,13 +76,4 @@ namespace corsika {
     }
   }
 
-  template <typename TDerived>
-  inline BaseExponential<TDerived>::BaseExponential(Point const& point,
-                                                    MassDensityType rho0,
-                                                    LengthType lambda)
-      : rho0_(rho0)
-      , lambda_(lambda)
-      , invLambda_(1 / lambda)
-      , point_(point) {}
-
 } // namespace corsika
diff --git a/corsika/detail/media/BaseTabular.inl b/corsika/detail/media/BaseTabular.inl
new file mode 100644
index 0000000000000000000000000000000000000000..3a02f63ad18eb3a89ab1cfdf1c7c134ccc846dd5
--- /dev/null
+++ b/corsika/detail/media/BaseTabular.inl
@@ -0,0 +1,200 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+
+#include <exception>
+
+namespace corsika {
+
+  template <typename TDerived>
+  inline BaseTabular<TDerived>::BaseTabular(
+      Point const& point, LengthType const referenceHeight,
+      std::function<MassDensityType(LengthType)> const& rho, unsigned int const nBins,
+      LengthType const deltaHeight)
+      : nBins_(nBins)
+      , deltaHeight_(deltaHeight)
+      , point_(point)
+      , referenceHeight_(referenceHeight) {
+    density_.resize(nBins_);
+    for (unsigned int bin = 0; bin < nBins; ++bin) {
+      density_[bin] = rho(deltaHeight_ * bin);
+      CORSIKA_LOG_DEBUG("new tabulated atm bin={} h={} rho={}", bin, deltaHeight_ * bin,
+                        density_[bin]);
+    }
+  }
+
+  template <typename TDerived>
+  inline auto const& BaseTabular<TDerived>::getImplementation() const {
+    return *static_cast<TDerived const*>(this);
+  }
+
+  template <typename TDerived>
+  inline MassDensityType BaseTabular<TDerived>::getMassDensity(
+      LengthType const height) const {
+    double const fbin = (height - referenceHeight_) / deltaHeight_;
+    int const bin = int(fbin);
+    if (bin < 0) return MassDensityType::zero();
+    if (bin >= int(nBins_ - 1)) {
+      CORSIKA_LOG_ERROR(
+          "invalid height {} (corrected {}) in BaseTabular atmosphere. Min 0, max {}. If "
+          "max is too low: increase!",
+          height, height - referenceHeight_, nBins_ * deltaHeight_);
+      throw std::runtime_error("invalid height");
+    }
+    return density_[bin] + (fbin - bin) * (density_[bin + 1] - density_[bin]);
+  }
+
+  template <typename TDerived>
+  inline GrammageType BaseTabular<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+
+    Point pCurr = traj.getPosition(0);
+    DirectionVector dCurr = traj.getDirection(0);
+    LengthType height1 = (traj.getPosition(0) - point_).getNorm() - referenceHeight_;
+    LengthType height2 = (traj.getPosition(1) - point_).getNorm() - referenceHeight_;
+
+    LengthType const fullLength = traj.getLength(1);
+
+    int sign = +1; // normal direction
+    if (height1 > height2) {
+      std::swap(height1, height2);
+      pCurr = traj.getPosition(1);
+      dCurr = traj.getDirection(1);
+      sign = -1; // inverted direction
+    }
+
+    double const fbin1 = height1 / deltaHeight_;
+    unsigned int const bin1 = int(fbin1);
+
+    double const fbin2 = height2 / deltaHeight_;
+    unsigned int const bin2 = int(fbin2);
+
+    if (fbin1 == fbin2) { return GrammageType::zero(); }
+
+    if (bin1 >= nBins_ - 1 || bin2 >= nBins_ - 1) {
+      CORSIKA_LOG_ERROR("invalid height {} {} in BaseTabular atmosphere integration",
+                        height1, height2);
+      throw std::runtime_error("invalid height");
+    }
+
+    // interpolated start/end densities
+    MassDensityType const rho1 = getMassDensity(height1 + referenceHeight_);
+    MassDensityType const rho2 = getMassDensity(height2 + referenceHeight_);
+
+    // within first bin
+    if (bin1 == bin2) { return fullLength * (rho2 + rho1) / 2; }
+
+    // inclination of trajectory (local)
+    DirectionVector axis((pCurr - point_).normalized()); // to gravity center
+    double cosTheta = axis.dot(dCurr);
+
+    // distance to next height bin
+    unsigned int bin = bin1;
+    LengthType dD = (deltaHeight_ * (bin + 1) - height1) / cosTheta * sign;
+    LengthType distance = dD;
+
+    GrammageType X = dD * (rho1 + density_[bin + 1]) / 2;
+    double frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength);
+    pCurr = traj.getPosition(frac);
+    dCurr = traj.getDirection(frac);
+
+    for (++bin; bin < bin2; ++bin) {
+      // inclination of trajectory
+      axis = (pCurr - point_).normalized();
+      cosTheta = axis.dot(dCurr);
+      // distance to next height bin
+      dD = deltaHeight_ / cosTheta * sign;
+      distance += dD;
+      GrammageType const dX = dD * (density_[bin] + density_[bin + 1]) / 2;
+      X += dX;
+      frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength);
+      pCurr = traj.getPosition(frac);
+      dCurr = traj.getDirection(frac);
+    }
+
+    // inclination of trajectory
+    axis = ((pCurr - point_).normalized());
+    cosTheta = axis.dot(dCurr);
+    // distance to next height bin
+    dD = (height2 - deltaHeight_ * bin2) / cosTheta * sign;
+    X += dD * (rho2 + density_[bin2]) / 2;
+    distance += dD;
+    return X;
+  }
+
+  template <typename TDerived>
+  inline LengthType BaseTabular<TDerived>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+
+    if (grammage < GrammageType::zero()) {
+      CORSIKA_LOG_ERROR("cannot integrate negative grammage");
+      throw std::runtime_error("negative grammage error");
+    }
+    LengthType const height = (traj.getPosition(0) - point_).getNorm() - referenceHeight_;
+
+    double const fbin = height / deltaHeight_;
+    int bin = int(fbin);
+
+    if (bin >= int(nBins_ - 1)) {
+      CORSIKA_LOG_ERROR("invalid height {} in BaseTabular atmosphere integration",
+                        height);
+      throw std::runtime_error("invalid height");
+    }
+
+    // interpolated start/end densities
+    MassDensityType const rho = getMassDensity(height + referenceHeight_);
+
+    // inclination of trajectory
+    Point pCurr = traj.getPosition(0);
+    DirectionVector dCurr = traj.getDirection(0);
+    DirectionVector axis((pCurr - point_).normalized());
+    double cosTheta = axis.dot(dCurr);
+    int sign = +1; // height increasing along traj
+    if (cosTheta < 0) {
+      cosTheta = -cosTheta; // absolute value only
+      sign = -1;            // height decreasing along traj
+    }
+
+    // height -> distance
+    LengthType const deltaDistance = deltaHeight_ / cosTheta;
+
+    // start with 0 g/cm2
+    GrammageType X(GrammageType::zero());
+    LengthType distance(LengthType::zero());
+
+    // within first bin
+    distance =
+        (sign > 0 ? deltaDistance * (bin + 1 - fbin) : deltaDistance * (fbin - bin));
+    GrammageType binGrammage = (sign > 0 ? distance * (rho + density_[bin + 1]) / 2
+                                         : distance * (rho + density_[bin]) / 2);
+    if (X + binGrammage > grammage) {
+      double const binFraction = (grammage - X) / binGrammage;
+      return distance * binFraction;
+    }
+    X += binGrammage;
+
+    // the following bins (along trajectory)
+    for (bin += sign; bin < int(nBins_) && bin >= 0; bin += sign) {
+
+      binGrammage = deltaDistance * (density_[bin] + density_[bin + 1]) / 2;
+      if (X + binGrammage > grammage) {
+        double const binFraction = (grammage - X) / binGrammage;
+        return distance + deltaDistance * binFraction;
+      }
+      X += binGrammage;
+      distance += deltaDistance;
+    }
+    return std::numeric_limits<double>::infinity() * meter;
+  }
+
+} // namespace corsika
diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl
index b64b2ded2f3546886f490fcfc9e3a37906609d01..f665dda992148f252261327c5a7f931cf5a9416c 100644
--- a/corsika/detail/media/FlatExponential.inl
+++ b/corsika/detail/media/FlatExponential.inl
@@ -18,19 +18,17 @@ namespace corsika {
 
   template <typename T>
   inline FlatExponential<T>::FlatExponential(Point const& point,
-                                             Vector<dimensionless_d> const& axis,
+                                             DirectionVector const& axis,
                                              MassDensityType rho, LengthType lambda,
                                              NuclearComposition const& nuclComp)
-      : BaseExponential<FlatExponential<T>>(point, rho, lambda)
+      : BaseExponential<FlatExponential<T>>(point, 0_m, rho, lambda)
       , axis_(axis)
       , nuclComp_(nuclComp) {}
 
   template <typename T>
   inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
-    return BaseExponential<FlatExponential<T>>::getRho0() *
-           exp(BaseExponential<FlatExponential<T>>::getInvLambda() *
-               (point - BaseExponential<FlatExponential<T>>::getAnchorPoint())
-                   .dot(axis_));
+    return BaseExponential<FlatExponential<T>>::getMassDensity(
+        (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm());
   }
 
   template <typename T>
@@ -40,8 +38,8 @@ namespace corsika {
 
   template <typename T>
   inline GrammageType FlatExponential<T>::getIntegratedGrammage(
-      BaseTrajectory const& line, LengthType to) const {
-    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
+      BaseTrajectory const& line) const {
+    return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
   }
 
   template <typename T>
diff --git a/corsika/detail/media/HomogeneousMedium.inl b/corsika/detail/media/HomogeneousMedium.inl
index 368073c9f61f2bed7f05f907128940cb16ef1d7f..032209bd3739ab5fe3d15ec3df23a9a5767f3ff9 100644
--- a/corsika/detail/media/HomogeneousMedium.inl
+++ b/corsika/detail/media/HomogeneousMedium.inl
@@ -32,9 +32,9 @@ namespace corsika {
   }
 
   template <typename T>
-  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(BaseTrajectory const&,
-                                                                  LengthType to) const {
-    return to * density_;
+  inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+      BaseTrajectory const& track) const {
+    return track.getLength() * density_;
   }
 
   template <typename T>
diff --git a/corsika/detail/media/InhomogeneousMedium.inl b/corsika/detail/media/InhomogeneousMedium.inl
index 8706765addf54494958788dd5709a7bc2ab9094b..f53f2c3be92471c0652d4921dfb2f14b457e6c90 100644
--- a/corsika/detail/media/InhomogeneousMedium.inl
+++ b/corsika/detail/media/InhomogeneousMedium.inl
@@ -36,8 +36,8 @@ namespace corsika {
 
   template <typename T, typename TDensityFunction>
   inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
-      BaseTrajectory const& line, LengthType to) const {
-    return densityFunction_.getIntegrateGrammage(line, to);
+      BaseTrajectory const& line) const {
+    return densityFunction_.getIntegrateGrammage(line);
   }
 
   template <typename T, typename TDensityFunction>
diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index 74f00f101bd284fc79f06e8f0e7c02a655b197c7..17371ff3b775986099e6436894c1213a945e4733 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
@@ -13,6 +13,7 @@
 #include <corsika/media/FlatExponential.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
 #include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
 
 namespace corsika {
 
@@ -41,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;
 
@@ -55,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
@@ -63,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));
@@ -74,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) {
@@ -93,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>>(
@@ -103,6 +103,40 @@ namespace corsika {
     layers_.push(std::move(node));
   }
 
+  template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+            typename... TModelArgs>
+  inline void
+  LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+      addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+                      unsigned int const nBins, LengthType const deltaHeight,
+                      LengthType const 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_, 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));
+  }
+
   template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
             typename... TModelArgs>
   inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
@@ -132,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};
     }
   };
 
diff --git a/corsika/detail/media/LinearApproximationIntegrator.inl b/corsika/detail/media/LinearApproximationIntegrator.inl
index ced2dee38b7f477c7c81d9bc80c1fe27024094b6..8ee60d9844b200b9d46ebb71af938698a2b6d2f2 100644
--- a/corsika/detail/media/LinearApproximationIntegrator.inl
+++ b/corsika/detail/media/LinearApproximationIntegrator.inl
@@ -19,10 +19,14 @@ namespace corsika {
 
   template <typename TDerived>
   inline auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
-      BaseTrajectory const& line, LengthType length) const {
+      BaseTrajectory const& line) const {
+    LengthType const length = line.getLength();
     auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
     auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
                                                                 line.getDirection(0));
+    CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1,
+                     line.getPosition(0), line.getDirection(0),
+                     (c0 + 0.5 * c1 * length) * length);
     return (c0 + 0.5 * c1 * length) * length;
   }
 
diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl
index 61911772a432e5d4d1c5b6cec3fc13f9ea2c8cdf..9b921872811e157c9d0d232286b3714444983727 100644
--- a/corsika/detail/media/SlidingPlanarExponential.inl
+++ b/corsika/detail/media/SlidingPlanarExponential.inl
@@ -8,52 +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, 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 height =
-        (point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-            .getNorm() -
-        referenceHeight_;
-    return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
-           exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * height);
+    auto const heightFull =
+        (point - BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+            .getNorm();
+    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(
-      BaseTrajectory const& traj, LengthType l) const {
-    auto const axis = (traj.getPosition(0) -
-                       BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
-                          .normalized();
-    return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, l,
-                                                                               axis);
+  template <typename TDerived>
+  inline GrammageType SlidingPlanarExponential<TDerived>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    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);
   }
 
diff --git a/corsika/detail/media/SlidingPlanarTabular.inl b/corsika/detail/media/SlidingPlanarTabular.inl
new file mode 100644
index 0000000000000000000000000000000000000000..c8e15792cafea62f7ba7bb05b8388fdca2fa0608
--- /dev/null
+++ b/corsika/detail/media/SlidingPlanarTabular.inl
@@ -0,0 +1,47 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  template <typename T>
+  inline SlidingPlanarTabular<T>::SlidingPlanarTabular(
+      Point const& p0, std::function<MassDensityType(LengthType)> const& rho,
+      unsigned int const nBins, LengthType const deltaHeight,
+      NuclearComposition const& nuclComp, LengthType referenceHeight)
+      : BaseTabular<SlidingPlanarTabular<T>>(p0, referenceHeight, rho, nBins, deltaHeight)
+      , nuclComp_(nuclComp) {}
+
+  template <typename T>
+  inline MassDensityType SlidingPlanarTabular<T>::getMassDensity(
+      Point const& point) const {
+    auto const heightFull =
+        (point - BaseTabular<SlidingPlanarTabular<T>>::getAnchorPoint()).getNorm();
+    return BaseTabular<SlidingPlanarTabular<T>>::getMassDensity(heightFull);
+  }
+
+  template <typename T>
+  inline NuclearComposition const& SlidingPlanarTabular<T>::getNuclearComposition()
+      const {
+    return nuclComp_;
+  }
+
+  template <typename T>
+  inline GrammageType SlidingPlanarTabular<T>::getIntegratedGrammage(
+      BaseTrajectory const& traj) const {
+    return BaseTabular<SlidingPlanarTabular<T>>::getIntegratedGrammage(traj);
+  }
+
+  template <typename T>
+  inline LengthType SlidingPlanarTabular<T>::getArclengthFromGrammage(
+      BaseTrajectory const& traj, GrammageType const grammage) const {
+    return BaseTabular<SlidingPlanarTabular<T>>::getArclengthFromGrammage(traj, grammage);
+  }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 7772f26609fdd65a3a1216235f4788ba5eda31a8..537a3d85fbf347303d00a1f519fa92ed8317e4d7 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -154,8 +154,7 @@ namespace corsika {
 
     if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
 
-    GrammageType const dX =
-        p.getNode()->getModelProperties().getIntegratedGrammage(t, t.getLength());
+    GrammageType const dX = p.getNode()->getModelProperties().getIntegratedGrammage(t);
     CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
                       p.getChargeNumber(), dX / 1_g * square(1_cm));
     HEPEnergyType dE = getTotalEnergyLoss(p, dX);
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index aa50f5866f008f1611489d33b7ebcb2a1aa736ab..30f0b10a1925d78fcdd7ad178b12accf31ffb358 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -97,8 +97,7 @@ namespace corsika::proposal {
     if (vT.getLength() == 0_m) return ProcessReturn::Ok;
 
     // calculate passed grammage
-    auto dX =
-        vP.getNode()->getModelProperties().getIntegratedGrammage(vT, vT.getLength());
+    auto dX = vP.getNode()->getModelProperties().getIntegratedGrammage(vT);
 
     // get or build corresponding track integral calculator and solve the
     // integral
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 4b40f19f40da8eba02ba6dac1735d913c641cd91..9c9bcfe1ac5f84824c446f4ee0dc9f0f201c5808 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -88,7 +88,8 @@ namespace corsika::urqmd {
         break;
       default: { // LCOV_EXCL_START since this can never happen due to canInteract
         CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
-        return CrossSectionType::zero(); // LCOV_EXCL_STOP
+        return CrossSectionType::zero();
+        // LCOV_EXCL_STOP
       }
     }
 
@@ -395,9 +396,9 @@ namespace corsika::urqmd {
     boost::filesystem::ifstream file(filename, std::ios::in);
 
     if (!file.is_open()) {
-      throw std::runtime_error(
-          filename.native() +
-          " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
+      // LCOV_EXCL_START since this is pointless to test
+      throw std::runtime_error(filename.native() + " could not be opened.");
+      // LCOV_EXCL_STOP
     }
 
     std::string line;
diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp
index 66f907943386d6f0d56f24577d39c5671cb00ab0..507ed640fa130eda007a5e71a3dec3437a530d9b 100644
--- a/corsika/media/BaseExponential.hpp
+++ b/corsika/media/BaseExponential.hpp
@@ -23,9 +23,18 @@ namespace corsika {
    */
   template <typename TDerived>
   class BaseExponential {
+
+  public:
+    BaseExponential(Point const& point, LengthType const referenceHeight,
+                    MassDensityType rho0, LengthType lambda);
+
+    Point const& getAnchorPoint() const { return point_; }
+
   protected:
     auto const& getImplementation() const;
 
+    MassDensityType getMassDensity(LengthType const height) const;
+
     // clang-format off
     /**
      * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
@@ -40,7 +49,7 @@ namespace corsika {
      * \f]
      */
     // clang-format on
-    GrammageType getIntegratedGrammage(BaseTrajectory const& line, LengthType vL,
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line,
                                        DirectionVector const& axis) const;
 
     // clang-format off
@@ -64,18 +73,12 @@ namespace corsika {
     LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage,
                                         DirectionVector const& axis) const;
 
-  public:
-    BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda);
-
-    Point const& getAnchorPoint() const { return point_; }
-    MassDensityType getRho0() const { return rho0_; }
-    InverseLengthType getInvLambda() const { return invLambda_; }
-
   private:
     MassDensityType const rho0_;
     LengthType const lambda_;
     InverseLengthType const invLambda_;
     Point const point_;
+    LengthType const referenceHeight_;
 
   }; // class BaseExponential
 
diff --git a/corsika/media/BaseTabular.hpp b/corsika/media/BaseTabular.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..cecf7cc8a03cb68905217ac2300500e6f69e7e48
--- /dev/null
+++ b/corsika/media/BaseTabular.hpp
@@ -0,0 +1,57 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
+
+#include <limits>
+#include <vector>
+
+namespace corsika {
+
+  /**
+   * This class provides the grammage/length conversion functionality for
+   * (locally) flat tabulated atmospheres.
+   */
+  template <typename TDerived>
+  class BaseTabular {
+
+  public:
+    BaseTabular(Point const& point, LengthType const referenceHeight,
+                std::function<MassDensityType(LengthType)> const& rho,
+                unsigned int const nBins, LengthType const deltaHeight);
+
+    Point const& getAnchorPoint() const { return point_; }
+
+  protected:
+    auto const& getImplementation() const;
+
+    MassDensityType getMassDensity(LengthType const height) const;
+
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const;
+
+    LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+                                        GrammageType const grammage) const;
+
+  private:
+    unsigned int const nBins_;
+    LengthType deltaHeight_;
+    std::vector<MassDensityType> density_;
+    Point const point_;
+    LengthType const referenceHeight_;
+
+  }; // class BaseTabular
+
+} // namespace corsika
+
+#include <corsika/detail/media/BaseTabular.inl>
diff --git a/corsika/media/FlatExponential.hpp b/corsika/media/FlatExponential.hpp
index 71767d705a030db91179ef70bb3441cfaf5c9d2d..6f64295e83efcc36a86ce4e43cc94a37e6593d78 100644
--- a/corsika/media/FlatExponential.hpp
+++ b/corsika/media/FlatExponential.hpp
@@ -41,8 +41,7 @@ namespace corsika {
 
     NuclearComposition const& getNuclearComposition() const override;
 
-    GrammageType getIntegratedGrammage(BaseTrajectory const& line,
-                                       LengthType to) const override;
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
 
     LengthType getArclengthFromGrammage(BaseTrajectory const& line,
                                         GrammageType grammage) const override;
diff --git a/corsika/media/HomogeneousMedium.hpp b/corsika/media/HomogeneousMedium.hpp
index 95a9fb0ae91ee41aa8b8f6036afbda3f105c8207..16331957db5334d042e0f4f4ff6e8f67d1fd471b 100644
--- a/corsika/media/HomogeneousMedium.hpp
+++ b/corsika/media/HomogeneousMedium.hpp
@@ -30,8 +30,7 @@ namespace corsika {
 
     NuclearComposition const& getNuclearComposition() const override;
 
-    GrammageType getIntegratedGrammage(BaseTrajectory const&,
-                                       LengthType to) const override;
+    GrammageType getIntegratedGrammage(BaseTrajectory const&) const override;
 
     LengthType getArclengthFromGrammage(BaseTrajectory const&,
                                         GrammageType grammage) const override;
diff --git a/corsika/media/IMediumModel.hpp b/corsika/media/IMediumModel.hpp
index 87e702c5dba8cf7aafdbc0ac3fe0f5f9ba090e75..e3e030392ff949ce41e3758b64eb975e06727ee3 100644
--- a/corsika/media/IMediumModel.hpp
+++ b/corsika/media/IMediumModel.hpp
@@ -21,11 +21,25 @@ namespace corsika {
 
     virtual MassDensityType getMassDensity(Point const&) const = 0;
 
-    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
-    // approach; for now, only lines are supported
-    virtual GrammageType getIntegratedGrammage(BaseTrajectory const&,
-                                               LengthType) const = 0;
-
+    /**
+     * Integrate the matter density along trajectory.
+     *
+     * @return GrammageType as integrated matter density along the BaseTrajectory
+     *
+     * @todo think about the mixin inheritance of the trajectory vs the BaseTrajectory
+     *       approach; for now, only lines are supported (?).
+     */
+    virtual GrammageType getIntegratedGrammage(BaseTrajectory const&) const = 0;
+
+    /**
+     * Calculates the length along the trajectory.
+     *
+     * The length along the trajectory is determined at which the integrated matter
+     * density is reached. If the specified matter density cannot be reached (is too
+     * large) the result becomes meaningless and could be "infinity" (discuss this).
+     *
+     * @return LengthType the length corresponding to grammage.
+     */
     virtual LengthType getArclengthFromGrammage(BaseTrajectory const&,
                                                 GrammageType) const = 0;
 
diff --git a/corsika/media/InhomogeneousMedium.hpp b/corsika/media/InhomogeneousMedium.hpp
index a257bc508969c4620172b4dae45e0cf52fa06ae3..66cf00f731af3a3ef9d4940ce6de9c10f3b63287 100644
--- a/corsika/media/InhomogeneousMedium.hpp
+++ b/corsika/media/InhomogeneousMedium.hpp
@@ -32,8 +32,7 @@ namespace corsika {
 
     NuclearComposition const& getNuclearComposition() const override;
 
-    GrammageType getIntegratedGrammage(BaseTrajectory const& line,
-                                       LengthType to) const override;
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
 
     LengthType getArclengthFromGrammage(BaseTrajectory const& pLine,
                                         GrammageType grammage) const override;
diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
index 5de160d136e54f61537d9c49943cc863137bb7dd..6c178f4589537fff2a79c7e166662791c30a1306 100644
--- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+++ b/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,15 +79,19 @@ namespace corsika {
     void addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary);
     void addLinearLayer(LengthType c, 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(); }
 
     void assemble(Environment<TMediumInterface>& env);
     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;
@@ -94,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>
diff --git a/corsika/media/LinearApproximationIntegrator.hpp b/corsika/media/LinearApproximationIntegrator.hpp
index 0e2eb0a295cc373c5749acb444b13a206c05622d..ef442ebdb10882d82b9c8ed02323f76c5de20148 100644
--- a/corsika/media/LinearApproximationIntegrator.hpp
+++ b/corsika/media/LinearApproximationIntegrator.hpp
@@ -20,7 +20,7 @@ namespace corsika {
     auto const& getImplementation() const;
 
   public:
-    auto getIntegrateGrammage(BaseTrajectory const& line, LengthType length) const;
+    auto getIntegrateGrammage(BaseTrajectory const& line) const;
 
     auto getArclengthFromGrammage(BaseTrajectory const& line,
                                   GrammageType grammage) const;
diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp
index a06c293a4d1056e74b40f52bb7fac7f3cf3a0070..4ef91bb4e1ca338241c61d805f1fa3dd68c47ba4 100644
--- a/corsika/media/SlidingPlanarExponential.hpp
+++ b/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>
 
@@ -46,15 +45,13 @@ namespace corsika {
 
     NuclearComposition const& getNuclearComposition() const override;
 
-    GrammageType getIntegratedGrammage(BaseTrajectory const& line,
-                                       LengthType l) const override;
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
 
     LengthType getArclengthFromGrammage(BaseTrajectory const& line,
                                         GrammageType grammage) const override;
 
   private:
     NuclearComposition const nuclComp_;
-    LengthType const referenceHeight_;
   };
 
 } // namespace corsika
diff --git a/corsika/media/SlidingPlanarTabular.hpp b/corsika/media/SlidingPlanarTabular.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..afd0d88e259723dcb6f77abe7471a75701f78298
--- /dev/null
+++ b/corsika/media/SlidingPlanarTabular.hpp
@@ -0,0 +1,61 @@
+/*
+ * (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.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/BaseTrajectory.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/BaseTabular.hpp>
+
+namespace corsika {
+
+  // clang-format off
+  /**
+   * The SlidingPlanarTabular models mass density as
+   * \f[
+   *   \varrho(r) = \varrho_0 \rho\left( |p_0 - r| \right).
+   * \f]
+   * For grammage/length conversion, the density distribution is approximated as
+   * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the 
+   * axis pointing rom \f$ p_0 \f$ to \f$ r_0 \f$ defining the local height.
+   */
+  // clang-format on
+
+  template <typename TDerived>
+  class SlidingPlanarTabular : public BaseTabular<SlidingPlanarTabular<TDerived>>,
+                               public TDerived {
+
+    using Base = BaseTabular<SlidingPlanarTabular<TDerived>>;
+
+  public:
+    SlidingPlanarTabular(Point const& p0,
+                         std::function<MassDensityType(LengthType)> const& rho,
+                         unsigned int const nBins, LengthType const deltaHeight,
+                         NuclearComposition const& nuclComp,
+                         LengthType referenceHeight = LengthType::zero());
+
+    MassDensityType getMassDensity(Point const& point) const override;
+
+    NuclearComposition const& getNuclearComposition() const override;
+
+    GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
+
+    LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+                                        GrammageType grammage) const override;
+
+  private:
+    NuclearComposition const nuclComp_;
+  };
+
+} // namespace corsika
+
+#include <corsika/detail/media/SlidingPlanarTabular.inl>
diff --git a/corsika/stack/DummyStack.hpp b/corsika/stack/DummyStack.hpp
index 9bbd51403f98fd53b9d64c92460bcbaeafece400..bcba30f1d55d2ca5197ca804d22ff672096c9d16 100644
--- a/corsika/stack/DummyStack.hpp
+++ b/corsika/stack/DummyStack.hpp
@@ -71,6 +71,8 @@ namespace corsika::dummy_stack {
      */
     void copy(const int /*i1*/, const int /*i2*/) {}
 
+    void swap(const int, const int) {}
+
     void incrementSize() { entries_++; }
     void decrementSize() { entries_--; }
 
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index e3d056e22202e7af786c3f1eec00fcc36e077a61..70816c0d30968b67caa65f98b204c4ab0b09c338 100644
--- a/examples/CMakeLists.txt
+++ b/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)
 
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index fd048470df3c3d4c598ab35c0be67b3db0ebf0e1..275ca88e26daebf256a57ce540152e3cb21957d1 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -97,14 +97,8 @@ void registerRandomStreams(int seed) {
 template <typename T>
 using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
 
-// argv : 1.number of nucleons, 2.number of protons,
-//        3.total energy in GeV, 4.number of showers,
-//        5.seed (0 by default to generate random values for all)
-
 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 +148,30 @@ 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")
+      ->default_str("info")
+      ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
+      ->group("Misc.");
 
   // parse the command line options into the variables
   CLI11_PARSE(app, argc, argv);
 
+  string const loglevel =
+      (app.count("--verbosity") ? app["--verbosity"]->as<string>() : "info");
+  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 +199,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 += 10_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 +257,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 +307,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);
@@ -311,7 +337,7 @@ int main(int argc, char** argv) {
 
   // observation plane
   Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
-  ObservationPlane<setup::Tracking, NoOutput> observationLevel(
+  ObservationPlane<setup::Tracking> observationLevel(
       obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
   // register the observation plane with the output
   output.add("particles", observationLevel);
@@ -350,9 +376,10 @@ int main(int argc, char** argv) {
     output.startOfShower();
 
     // directory for outputs
-    string const labHist_file = "inthist_lab_verticalEAS_" + to_string(i_shower) + ".npz";
-    string const cMSHist_file = "inthist_cms_verticalEAS_" + to_string(i_shower) + ".npz";
-    string const longprof_file = "longprof_verticalEAS_" + to_string(i_shower) + ".txt";
+    string const outdir(app["--filename"]->as<std::string>());
+    string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
+    string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
+    string const longprof_file = outdir + "/longprof_" + to_string(i_shower) + ".txt";
 
     // setup particle stack, and add primary particle
     stack.clear();
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 0ff4d28b7425d1fd988af0718ccecbe84b4b8ba1..d8b8385f3696ff1fee4f48518253bd24e97fbdb0 100644
--- a/examples/em_shower.cpp
+++ b/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));
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index e065c4ee8ddc787e2bd56d2f8502f37cfb7d0bd1..4674512cc4d0f204538999c6bf6a3eb384704ea2 100644
--- a/examples/hybrid_MC.cpp
+++ b/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));
diff --git a/examples/mars.cpp b/examples/mars.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..92ed62ad9cd76578ebbbaeddf4c42702e1b64416
--- /dev/null
+++ b/examples/mars.cpp
@@ -0,0 +1,457 @@
+/*
+ * (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.
+ */
+
+/* clang-format off */
+// InteractionCounter used boost/histogram, which
+// fails if boost/type_traits have been included before. Thus, we have
+// to include it first...
+#include <corsika/framework/process/InteractionCounter.hpp>
+/* clang-format on */
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/SaveBoostHistogram.hpp>
+#include <corsika/framework/process/ProcessSequence.hpp>
+#include <corsika/framework/process/SwitchProcessSequence.hpp>
+#include <corsika/framework/process/InteractionCounter.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/CorsikaFenv.hpp>
+#include <corsika/framework/core/Cascade.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+#include <corsika/output/OutputManager.hpp>
+#include <corsika/output/NoOutput.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/SlidingPlanarExponential.hpp>
+
+#include <corsika/modules/BetheBlochPDG.hpp>
+#include <corsika/modules/LongitudinalProfile.hpp>
+#include <corsika/modules/ObservationPlane.hpp>
+#include <corsika/modules/OnShellCheck.hpp>
+#include <corsika/modules/StackInspector.hpp>
+#include <corsika/modules/TrackWriter.hpp>
+#include <corsika/modules/ParticleCut.hpp>
+#include <corsika/modules/Pythia8.hpp>
+#include <corsika/modules/Sibyll.hpp>
+#include <corsika/modules/UrQMD.hpp>
+#include <corsika/modules/PROPOSAL.hpp>
+#include <corsika/modules/QGSJetII.hpp>
+
+#include <corsika/setup/SetupStack.hpp>
+#include <corsika/setup/SetupTrajectory.hpp>
+
+#include <CLI/App.hpp>
+#include <CLI/Formatter.hpp>
+#include <CLI/Config.hpp>
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <string>
+
+/*
+  NOTE, WARNING, ATTENTION
+
+  The .../Random.hpp implement the hooks of external modules to the C8 random
+  number generator. It has to occur excatly ONCE per linked
+  executable. If you include the header below multiple times and
+  link this togehter, it will fail.
+ */
+#include <corsika/modules/Random.hpp>
+
+using namespace corsika;
+using namespace std;
+
+using Particle = setup::Stack::particle_type;
+
+typedef decltype(1 * pascal) PressureType;
+typedef decltype(1 * degree_celsius) TemperatureType;
+
+class MarsAtmModel {
+public:
+  MarsAtmModel() = delete;
+  MarsAtmModel(PressureType a, InverseLengthType b, TemperatureType c,
+               decltype(1 * degree_celsius / 1_m) d)
+      : a_(a)
+      , b_(b)
+      , c_(c)
+      , d_(d) {}
+
+  MassDensityType operator()(LengthType height) const {
+    PressureType const pressure = a_ * exp(-b_ * height);
+    TemperatureType const temperature = -c_ - d_ * height + 273.1_K; // in K
+    constexpr decltype(square(1_m) / (square(1_s) * 1_K)) constant =
+        1000 * 0.1921 * square(1_m) / (square(1_s) * 1_K);
+    return pressure / (constant * temperature);
+  }
+
+private:
+  PressureType a_;
+  InverseLengthType b_;
+  TemperatureType c_;
+  decltype(1_K / 1_m) d_;
+};
+
+void registerRandomStreams(int seed) {
+  RNGManager<>::getInstance().registerRandomStream("cascade");
+  RNGManager<>::getInstance().registerRandomStream("qgsjet");
+  RNGManager<>::getInstance().registerRandomStream("sibyll");
+  RNGManager<>::getInstance().registerRandomStream("pythia");
+  RNGManager<>::getInstance().registerRandomStream("urqmd");
+  RNGManager<>::getInstance().registerRandomStream("proposal");
+  if (seed == 0) {
+    std::random_device rd;
+    seed = rd();
+    cout << "new random seed (auto) " << seed << endl;
+  }
+  RNGManager<>::getInstance().setSeed(seed);
+}
+
+template <typename T>
+using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
+
+int main(int argc, char** argv) {
+
+  // the main command line description
+  CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};
+
+  // some options that we want to fill in
+  int A, Z, nevent = 0;
+
+  // the following section adds the options to the parser
+
+  // we start by definining a sub-group for the primary ID
+  auto opt_Z = app.add_option("-Z", Z, "Atomic number for primary")
+                   ->check(CLI::Range(0, 26))
+                   ->group("Primary");
+  auto opt_A = app.add_option("-A", A, "Atomic mass number for primary")
+                   ->needs(opt_Z)
+                   ->check(CLI::Range(1, 58))
+                   ->group("Primary");
+  app.add_option("-p,--pdg", "PDG code for primary.")
+      ->excludes(opt_A)
+      ->excludes(opt_Z)
+      ->group("Primary");
+  // the remainding options
+  app.add_option("-E,--energy", "Primary energy in GeV")
+      ->required()
+      ->check(CLI::PositiveNumber)
+      ->group("Primary");
+  app.add_option("-z,--zenith", "Primary zenith angle (deg)")
+      ->required()
+      ->default_val(0.)
+      ->check(CLI::Range(0, 90))
+      ->group("Primary");
+  app.add_option("-a,--azimuth", "Primary azimuth angle (deg)")
+      ->default_val(0.)
+      ->check(CLI::Range(0, 360))
+      ->group("Primary");
+  app.add_option("-N,--nevent", nevent, "The number of events/showers to run.")
+      ->required()
+      ->check(CLI::PositiveNumber)
+      ->group("Library/Output");
+  app.add_option("-f,--filename", "Filename for output library.")
+      ->required()
+      ->default_val("corsika_library")
+      ->check(CLI::NonexistentPath)
+      ->group("Library/Output");
+  app.add_option("-s,--seed", "The random number seed.")
+      ->default_val(12351739)
+      ->check(CLI::NonNegativeNumber)
+      ->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);
+
+  string const loglevel =
+      (app.count("--verbosity") ? app["verbosity"]->as<string>() : "info");
+  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
+  if (app.count("--pdg") == 0) {
+    if ((app.count("-A") == 0) || (app.count("-Z") == 0)) {
+      std::cerr << "If --pdg is not provided, then both -A and -Z are required."
+                << std::endl;
+      return 1;
+    }
+  }
+
+  // initialize random number sequence(s)
+  registerRandomStreams(app["--seed"]->as<int>());
+
+  /* === START: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+  using EnvType = setup::Environment;
+  EnvType env;
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+  Point const center{rootCS, 0_m, 0_m, 0_m};
+  LengthType const radiusMars = 3389.5_km;
+  auto builder =
+      make_layered_spherical_atmosphere_builder<setup::EnvironmentInterface, MyExtraEnv>::
+          create(center,
+                 radiusMars,                                   // Mars
+                 Medium::AirDry1Atm,                           // Mars, close enough
+                 MagneticFieldVector{rootCS, 0_T, 0_uT, 0_T}); // Mars
+
+  builder.setNuclearComposition(                             // Mars
+      {{Code::Nitrogen, Code::Oxygen}, {1. / 3., 2. / 3.}}); // simplified
+  //{{Code::Carbon, Code::Oxygen, // 95.97 CO2
+  //          Code::Nitrogen},            // 1.89 N2 + 1.93 Argon + 0.146 O2
+  //       {0.9597 / 3, 0.9597 * 2 / 3,
+  //      1 - 0.9597}}); // values taken from AIRES manual, Ar removed for now
+
+  MarsAtmModel layer1(0.699e3 * pascal, 0.00009 / 1_m, 31.0 * degree_celsius,
+                      0.000998 * 1 * degree_celsius / 1_m);
+  MarsAtmModel layer2(0.699e3 * pascal, 0.00009 / 1_m, 23.4 * degree_celsius,
+                      0.00222 * 1 * degree_celsius / 1_m);
+
+  builder.addTabularLayer(layer1, 100, 100_m, 7_km);
+  builder.addTabularLayer(layer2, 300, 500_m, 100_km);
+  builder.addLinearLayer(1e9_cm, 112.8_km);
+  builder.assemble(env);
+  /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+
+  ofstream atmout("mars.dat");
+  for (LengthType h = 0_m; h < 110_km; h += 10_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 === */
+
+  // parse the primary ID as a PDG or A/Z code
+  Code beamCode;
+  HEPEnergyType mass;
+
+  // check if we want to use a PDG code instead
+  if (app.count("--pdg") > 0) {
+    beamCode = convert_from_PDG(PDGCode(app["--pdg"]->as<int>()));
+    mass = get_mass(beamCode);
+  } else {
+    // check manually for proton and neutrons
+    if ((A == 0) && (Z == 1)) beamCode = Code::Proton;
+    if ((A == 1) && (Z == 1)) beamCode = Code::Neutron;
+    mass = get_nucleus_mass(A, Z);
+  }
+
+  // particle energy
+  HEPEnergyType const E0 = 1_GeV * app["--energy"]->as<float>();
+
+  // direction of the shower in (theta, phi) space
+  auto const thetaRad = app["--zenith"]->as<float>() / 180. * M_PI;
+  auto const phiRad = app["--azimuth"]->as<float>() / 180. * M_PI;
+
+  // convert Elab to Plab
+  HEPMomentumType P0 = sqrt((E0 - mass) * (E0 + mass));
+
+  // convert the momentum to the zenith and azimuth angle of the primary
+  auto const [px, py, pz] =
+      std::make_tuple(P0 * sin(thetaRad) * cos(phiRad), P0 * sin(thetaRad) * sin(phiRad),
+                      -P0 * cos(thetaRad));
+  auto plab = MomentumVector(rootCS, {px, py, pz});
+  /* === END: CONSTRUCT PRIMARY PARTICLE === */
+
+  /* === START: CONSTRUCT GEOMETRY === */
+  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));
+  Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+  Point const injectionPos =
+      showerCore + DirectionVector{rootCS,
+                                   {-sin(thetaRad) * cos(phiRad),
+                                    -sin(thetaRad) * sin(phiRad), cos(thetaRad)}} *
+                       t;
+
+  // we make the axis much longer than the inj-core distance since the
+  // profile will go beyond the core, depending on zenith angle
+  ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos) * 1.2, env};
+  /* === END: CONSTRUCT GEOMETRY === */
+
+  // create the output manager that we then register outputs with
+  OutputManager output(app["--filename"]->as<std::string>());
+
+  /* === START: SETUP PROCESS LIST === */
+  corsika::sibyll::Interaction sibyll;
+  InteractionCounter sibyllCounted(sibyll);
+
+  corsika::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
+  InteractionCounter sibyllNucCounted(sibyllNuc);
+
+  corsika::pythia8::Decay decayPythia;
+
+  // use sibyll decay routine for decays of particles unknown to pythia
+  corsika::sibyll::Decay decaySibyll{{
+      Code::N1440Plus,
+      Code::N1440MinusBar,
+      Code::N1440_0,
+      Code::N1440_0Bar,
+      Code::N1710Plus,
+      Code::N1710MinusBar,
+      Code::N1710_0,
+      Code::N1710_0Bar,
+
+      Code::Pi1300Plus,
+      Code::Pi1300Minus,
+      Code::Pi1300_0,
+
+      Code::KStar0_1430_0,
+      Code::KStar0_1430_0Bar,
+      Code::KStar0_1430_Plus,
+      Code::KStar0_1430_MinusBar,
+  }};
+
+  // decaySibyll.printDecayConfig();
+
+  ParticleCut cut{1_GeV, 1_GeV, 1_GeV, 1_GeV, false};
+  corsika::proposal::Interaction emCascade(env);
+  corsika::proposal::ContinuousProcess emContinuous(env);
+  InteractionCounter emCascadeCounted(emCascade);
+
+  LongitudinalProfile longprof{showerAxis, 1_g / square(1_cm)};
+
+  corsika::urqmd::UrQMD urqmd;
+  InteractionCounter urqmdCounted{urqmd};
+  StackInspector<setup::Stack> stackInspect(5000, false, E0);
+
+  // assemble all processes into an ordered process list
+  struct EnergySwitch {
+    HEPEnergyType cutE_;
+    EnergySwitch(HEPEnergyType cutE)
+        : cutE_(cutE) {}
+    bool operator()(const Particle& p) { return (p.getKineticEnergy() < cutE_); }
+  };
+  auto hadronSequence = make_select(EnergySwitch(80_GeV), urqmdCounted,
+                                    make_sequence(sibyllNucCounted, sibyllCounted));
+  auto decaySequence = make_sequence(decayPythia, decaySibyll);
+
+  // track writer
+  TrackWriter trackWriter;
+  output.add("tracks", trackWriter); // register TrackWriter
+
+  // observation plane
+  Plane const obsPlane(showerCore, DirectionVector(rootCS, {0., 0., 1.}));
+  ObservationPlane<setup::Tracking> observationLevel(
+      obsPlane, DirectionVector(rootCS, {1., 0., 0.}));
+  // register the observation plane with the output
+  output.add("particles", observationLevel);
+
+  // assemble the final process sequence
+  auto sequence =
+      make_sequence(stackInspect, hadronSequence, decaySequence, emCascadeCounted,
+                    emContinuous, cut, trackWriter, observationLevel, longprof);
+  /* === END: SETUP PROCESS LIST === */
+
+  // create the cascade object using the default stack and tracking implementation
+  setup::Tracking tracking;
+  setup::Stack stack;
+  Cascade EAS(env, tracking, sequence, output, stack);
+
+  // print our primary parameters all in one place
+  if (app["--pdg"]->count() > 0) {
+    CORSIKA_LOG_INFO("Primary PDG ID: {}", app["--pdg"]->as<int>());
+  } else {
+    CORSIKA_LOG_INFO("Primary Z/A: {}/{}", Z, A);
+  }
+  CORSIKA_LOG_INFO("Primary Energy: {}", E0);
+  CORSIKA_LOG_INFO("Primary Momentum: {}", P0);
+  CORSIKA_LOG_INFO("Point of Injection: {}", injectionPos.getCoordinates());
+  CORSIKA_LOG_INFO("Shower Axis Length: {}", (showerCore - injectionPos).getNorm() * 1.2);
+
+  // trigger the output manager to open the library for writing
+  output.startOfLibrary();
+
+  // loop over each shower
+  for (int i_shower = 1; i_shower < nevent + 1; i_shower++) {
+    CORSIKA_LOG_INFO("Shower {} / {} ", i_shower, nevent);
+
+    // trigger the start of the outputs for this shower
+    output.startOfShower();
+
+    // directory for outputs
+    string const outdir(app["--filename"]->as<std::string>());
+    string const labHist_file = outdir + "/inthist_lab_" + to_string(i_shower) + ".npz";
+    string const cMSHist_file = outdir + "/inthist_cms_" + to_string(i_shower) + ".npz";
+    string const longprof_file = outdir + "/longprof" + to_string(i_shower) + ".txt";
+
+    // setup particle stack, and add primary particle
+    stack.clear();
+
+    // add the desired particle to the stack
+    if (A > 1) {
+      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns, A, Z));
+    } else {
+      stack.addParticle(std::make_tuple(beamCode, plab, injectionPos, 0_ns));
+    }
+
+    // if we want to fix the first location of the shower
+    if (app["--force-interaction"]) EAS.forceInteraction();
+
+    // run the shower
+    EAS.run();
+
+    cut.showResults();
+    // emContinuous.showResults();
+    observationLevel.showResults();
+    const HEPEnergyType Efinal = cut.getCutEnergy() + cut.getInvEnergy() +
+                                 cut.getEmEnergy() + // emContinuous.getEnergyLost() +
+                                 observationLevel.getEnergyGround();
+    cout << "total cut energy (GeV): " << Efinal / 1_GeV << endl
+         << "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
+    observationLevel.reset();
+    cut.reset();
+    // emContinuous.reset();
+
+    auto const hists = sibyllCounted.getHistogram() + sibyllNucCounted.getHistogram() +
+                       urqmdCounted.getHistogram();
+
+    save_hist(hists.labHist(), labHist_file, true);
+    save_hist(hists.CMSHist(), cMSHist_file, true);
+    longprof.save(longprof_file);
+
+    // trigger the output manager to save this shower to disk
+    output.endOfShower();
+  }
+
+  // and finalize the output on disk
+  output.endOfLibrary();
+}
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 6ae766bb197e2e6610aa503442cea05704652b3b..6800a84be8858d31ed01d36d213ef56754a265b1 100644
--- a/examples/vertical_EAS.cpp
+++ b/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));
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index b39fbb357757c28d51b6ff13b0a0f7dfceba480f..372cfdb9a578e321354517d3772c812920022fb5 100644
--- a/tests/media/testEnvironment.cpp
+++ b/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>
@@ -81,6 +83,9 @@ TEST_CASE("HomogeneousMedium") {
                                              std::vector<float>{1.f});
   HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
 
+  CHECK(protonComposition.getFractions() == std::vector<float>{1.});
+  CHECK(protonComposition.getComponents() == std::vector<Code>{Code::Proton});
+
   CHECK_THROWS(NuclearComposition({Code::Proton}, {1.1}));
   CHECK_THROWS(NuclearComposition({Code::Proton}, {0.99}));
 }
@@ -94,36 +99,42 @@ 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);
-  auto const tEnd = 5_s;
+  SpeedType const speed = 20_m / second;
+  LengthType const length = 2_m;
+  TimeType const tEnd = length / speed;
+
+  CHECK(medium.getNuclearComposition().getFractions() == std::vector<float>{1.});
+  CHECK(medium.getNuclearComposition().getComponents() ==
+        std::vector<Code>{Code::Proton});
 
   SECTION("horizontal") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {20_cm / second, 0_m / second, 0_m / second}));
+                                 gCS, {speed, 0_m / second, 0_m / second}));
     setup::Trajectory const trajectory =
         setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
-    CHECK((medium.getIntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
-    CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
+    CHECK((medium.getIntegratedGrammage(trajectory) / (rho0 * length)) == Approx(1));
+    CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * length) / length) ==
+          Approx(1));
   }
 
   SECTION("vertical") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 0_m / second, 5_m / second}));
+                                 gCS, {0_m / second, 0_m / second, speed}));
     setup::Trajectory const trajectory =
         setup::testing::make_track<setup::Trajectory>(line, tEnd);
-    LengthType const length = 2 * lambda;
     GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
 
-    CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1));
+    CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1));
     CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
   }
 
   SECTION("escape grammage") {
     Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 0_m / second, -5_m / second}));
+                                 gCS, {SpeedType::zero(), SpeedType::zero(), -speed}));
     setup::Trajectory const trajectory =
         setup::testing::make_track<setup::Trajectory>(line, tEnd);
     GrammageType const escapeGrammage = rho0 * lambda;
@@ -134,15 +145,15 @@ TEST_CASE("FlatExponential") {
   }
 
   SECTION("inclined") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 5_m / second, 5_m / second}));
+    Line const line(gOrigin,
+                    Vector<SpeedType::dimension_type>(
+                        gCS, {0_m / second, speed / sqrt(2.), speed / sqrt(2.)}));
     setup::Trajectory const trajectory =
         setup::testing::make_track<setup::Trajectory>(line, tEnd);
     double const cosTheta = M_SQRT1_2;
-    LengthType const length = 2 * lambda;
     GrammageType const exact =
         rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
-    CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1));
+    CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1));
     CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
   }
 }
@@ -179,16 +190,286 @@ TEST_CASE("SlidingPlanarExponential") {
 
     CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
           flat.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
-    CHECK(medium.getIntegratedGrammage(trajectory, 2_m).magnitude() ==
-          flat.getIntegratedGrammage(trajectory, 2_m).magnitude());
+    CHECK(medium.getIntegratedGrammage(trajectory).magnitude() ==
+          flat.getIntegratedGrammage(trajectory).magnitude());
     CHECK(medium.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
           flat.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
   }
 }
 
-auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
+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") {
 
-struct Exponential {
+  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("not possible") {
+    CHECK_THROWS(medium.getMassDensity({gCS, {0_m, 1e10_m, 0_m}}));
+
+    SpeedType const speed = 5_m / second;
+    TimeType const tEnd = 1e10_s;
+    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 trajectory =
+        setup::testing::make_track<setup::Trajectory>(line, tEnd);
+    CHECK_THROWS(medium.getIntegratedGrammage(trajectory));
+
+    Line const line2(
+        {gCS, {0_m, 0_m, 1e9_m}},
+        Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, speed}));
+    setup::Trajectory const trajectory2 =
+        setup::testing::make_track<setup::Trajectory>(line2, tEnd);
+    CHECK_THROWS(medium.getArclengthFromGrammage(trajectory2, 1e3_g / square(1_cm)));
+  }
+
+  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 {
   auto operator()(Point const& p) const {
     return exp(p.getCoordinates()[0] / 1_m) * rho0;
   }
@@ -213,41 +494,50 @@ TEST_CASE("InhomogeneousMedium") {
 
   Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
 
+  SpeedType const speed = 20_m / second;
   Line line(gOrigin, Vector<SpeedType::dimension_type>(
-                         gCS, {20_m / second, 0_m / second, 0_m / second}));
+                         gCS, {speed, SpeedType::zero(), SpeedType::zero()}));
 
-  auto const tEnd = 5_s;
+  // the tested LinearApproximationIntegrator really does a single step only. It is very
+  // poor for exponentials with a bit larger step-width.
+  TimeType const tEnd = 0.001_s;
   setup::Trajectory const trajectory =
       setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
-  Exponential const e;
-  DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
+  ExponentialTest const expTest;
+  DensityFunction<ExponentialTest, LinearApproximationIntegrator> const rho(expTest);
 
   SECTION("DensityFunction") {
-    CHECK(e.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+    CHECK(expTest.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
           Approx(1));
-    CHECK(rho.evaluateAt(gOrigin) == e(gOrigin));
+    CHECK(rho.evaluateAt(gOrigin) == expTest(gOrigin));
   }
 
   auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
   auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); };
 
-  auto constexpr l = 15_cm;
+  LengthType const length = tEnd * speed;
 
   NuclearComposition const composition{{Code::Proton}, {1.f}};
   InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
 
+  CORSIKA_LOG_INFO("test={} l={} {} {}", rho.getIntegrateGrammage(trajectory), length,
+                   exactGrammage(length), 1_m * rho0 * (exp(length / 1_m) - 1));
+
   SECTION("Integration") {
-    CHECK(rho.getIntegrateGrammage(trajectory, l) / exactGrammage(l) ==
+    CORSIKA_LOG_INFO("test={} {} {}", rho.getIntegrateGrammage(trajectory),
+                     exactGrammage(length),
+                     rho.getIntegrateGrammage(trajectory) / exactGrammage(length));
+    CHECK(rho.getIntegrateGrammage(trajectory) / exactGrammage(length) ==
           Approx(1).epsilon(1e-2));
-    CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(l)) /
-              exactLength(exactGrammage(l)) ==
+    CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(length)) /
+              exactLength(exactGrammage(length)) ==
           Approx(1).epsilon(1e-2));
     CHECK(rho.getMaximumLength(trajectory, 1e-2) >
-          l); // todo: write reasonable test when implementation is working
+          length); // todo: write reasonable test when implementation is working
 
-    CHECK(rho.getIntegrateGrammage(trajectory, l) ==
-          inhMedium.getIntegratedGrammage(trajectory, l));
+    CHECK(rho.getIntegrateGrammage(trajectory) ==
+          inhMedium.getIntegratedGrammage(trajectory));
     CHECK(rho.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
           inhMedium.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
     CHECK(inhMedium.getNuclearComposition() == composition);
@@ -278,7 +568,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
 
   CHECK(builder.getSize() == 0);
 
-  auto const R = builder.getEarthRadius();
+  auto const R = builder.getPlanetRadius();
 
   CHECK(univ->getChildNodes().size() == 1);
 
@@ -323,7 +613,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);
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index b3ccf4d75c4ef8f98aff716b2be559775c3a0a3f..6e9f347b8c7447cf227afdd191acc5362c7e1283 100644
--- a/tests/media/testMagneticField.cpp
+++ b/tests/media/testMagneticField.cpp
@@ -25,7 +25,6 @@ using namespace corsika;
 TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
 
   logging::set_level(logging::level::info);
-  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
 
   CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
   Point const gOrigin(gCS, {0_m, 0_m, 0_m});
@@ -85,8 +84,4 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
   // and the associated trajectory
   setup::Trajectory const track =
       setup::testing::make_track<setup::Trajectory>(line, tEnd);
-
-  // and check the integrated grammage
-  CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
-  CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
 }
diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp
index a10348717921ded38235b69a81328c694b0b87fe..bb5a9f0e86bb7398df1480844829fba1bf26f83c 100644
--- a/tests/media/testMedium.cpp
+++ b/tests/media/testMedium.cpp
@@ -92,18 +92,21 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
   CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
   medium.getNuclearComposition();
 
+  SpeedType const speed = 1_m / second;
+
   // create a line of length 1 m
-  Line const line(gOrigin,
-                  VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second}));
+  Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second}));
 
   // the end time of our line
   auto const tEnd = 1_s;
 
+  LengthType const length = tEnd * speed;
+
   // and the associated trajectory
   setup::Trajectory const trajectory =
       corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
   // and check the integrated grammage
-  CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+  CHECK((medium.getIntegratedGrammage(trajectory) / (density * length)) == Approx(1));
   CHECK((medium.getArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
 }
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index acef32c4f4fa5260ac841cf1db53a780f9f80900..62422f7a8c8d76da32f25a0972b7391ee7fb3ccf 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -76,19 +76,22 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
   CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
   medium.getNuclearComposition();
 
+  SpeedType const speed = 1_m / second;
+
   // create a line of length 1 m
-  Line const line(gOrigin,
-                  VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second}));
+  Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second}));
 
   // the end time of our line
   auto const tEnd = 1_s;
 
+  LengthType const length = tEnd * speed;
+
   // and the associated trajectory
   setup::Trajectory const track =
       setup::testing::make_track<setup::Trajectory>(line, tEnd);
 
   // and check the integrated grammage
-  CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  CHECK((medium.getIntegratedGrammage(track) / (density * length)) == Approx(1));
   CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
 }
 
@@ -148,24 +151,26 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
   REQUIRE(density == medium__.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
   medium__.getNuclearComposition();
 
-  // create a line of length 1 m
-  Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                               gCS, {1_m / second, 0_m / second, 0_m / second}));
+  SpeedType const velocity = 1_m / second;
 
   // the end time of our line
-  auto const tEnd = 1_s;
+  TimeType const tEnd = 1_s;
+
+  LengthType const length = tEnd * velocity;
+
+  // create a line of length 1 m
+  Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+                               gCS, {velocity, 0_m / second, 0_m / second}));
 
   // and the associated trajectory
   setup::Trajectory const track =
       setup::testing::make_track<setup::Trajectory>(line, tEnd);
-  //  // and the associated trajectory
-  //  Trajectory<Line> const trajectory(line, tEnd);
 
   // and check the integrated grammage
-  REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium.getIntegratedGrammage(track) / (density * length)) == Approx(1));
   REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
-  REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium_.getIntegratedGrammage(track) / (density * length)) == Approx(1));
   REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
-  REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium__.getIntegratedGrammage(track) / (density * length)) == Approx(1));
   REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
 }
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index a4c4ef5933bafa426e9e0266545368ba811dff64..aa2a006f29602efe6b66e841d304df29a357921c 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -38,6 +38,7 @@ TEST_CASE("epos", "modules") {
           corsika::epos::convertFromEpos(corsika::epos::EposCode::Electron));
     CHECK(Code::Proton ==
           corsika::epos::convertFromEpos(corsika::epos::EposCode::Proton));
+    CHECK_THROWS(corsika::epos::convertFromEpos(corsika::epos::EposCode::Unknown));
   }
 
   SECTION("corsika -> epos") {
diff --git a/tests/stack/testGeometryNodeStackExtension.cpp b/tests/stack/testGeometryNodeStackExtension.cpp
index 9cc29edcce53750f3ab1b4c5ba4b5f2def5885ad..2bc2e1160836ba85a39b1d93d7916c8459f93af5 100644
--- a/tests/stack/testGeometryNodeStackExtension.cpp
+++ b/tests/stack/testGeometryNodeStackExtension.cpp
@@ -86,4 +86,30 @@ TEST_CASE("GeometryNodeStackExtension", "stack") {
     CHECK(v == 99 * data);
     CHECK(s.getEntries() == 0);
   }
+
+  SECTION("copy and swap") {
+
+    const int data1 = 16;
+    const int data2 = 17;
+
+    TestStack s;
+    // add 99 particles, each 10th particle is a nucleus with A=i and Z=A/2!
+    for (int i = 0; i < 4; ++i) {
+      auto p = s.addParticle(std::tuple<dummy_stack::NoData>{noData});
+      p.setNode(i % 2 == 0 ? &data1 : &data2);
+    }
+
+    CHECK(*((s.first() + 0)).getNode() == 16);
+    CHECK(*((s.first() + 1)).getNode() == 17);
+    CHECK(*((s.first() + 2)).getNode() == 16);
+    CHECK(*((s.first() + 3)).getNode() == 17);
+
+    s.copy(s.first() + 0, s.first() + 1);
+    CHECK(*((s.first() + 0)).getNode() == 16);
+    CHECK(*((s.first() + 1)).getNode() == 16);
+
+    s.swap(s.first() + 3, s.first() + 1);
+    CHECK(*((s.first() + 1)).getNode() == 17);
+    CHECK(*((s.first() + 3)).getNode() == 16);
+  }
 }