diff --git a/Environment/Convenience.cpp b/Environment/Convenience.cpp
index fb2361c03570b0349f0f524677a23e33698fa170..f3b052d6d1101a7fbb2284dda15e78ac94a2e347 100644
--- a/Environment/Convenience.cpp
+++ b/Environment/Convenience.cpp
@@ -36,10 +36,11 @@ void LayeredSphericalAtmosphereBuilder::addExponentialLayer(
   auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
       std::make_unique<geometry::Sphere>(center_, radius));
 
-  auto const rho0 = b / c * exp(seaLevel_ / c);
+  auto const rho0 = b / c;
+  std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
 
-  node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(center_, rho0, -c,
-                                                                   *composition_);
+  node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(
+      center_, rho0, -c, *composition_, seaLevel_);
 
   layers_.push(std::move(node));
 }
@@ -55,6 +56,8 @@ void LayeredSphericalAtmosphereBuilder::addLinearLayer(
   units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
   auto const rho0 = b / c;
 
+  std::cout << "rho0 = " << rho0;
+
   auto node = std::make_unique<VolumeTreeNode<environment::IMediumModel>>(
       std::make_unique<geometry::Sphere>(center_, radius));
   node->SetModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
@@ -63,9 +66,13 @@ void LayeredSphericalAtmosphereBuilder::addLinearLayer(
 }
 
 Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
-  Environment<IMediumModel> environment_;
+  Environment<IMediumModel> env;
+  assemble(env);
+  return env;
+}
 
-  auto& universe = environment_.GetUniverse();
+void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
+  auto& universe = env.GetUniverse();
   auto* outmost = universe.get();
 
   while (!layers_.empty()) {
@@ -75,6 +82,4 @@ Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
     layers_.pop();
     outmost = tmp;
   }
-
-  return environment_;
 }
diff --git a/Environment/Convenience.h b/Environment/Convenience.h
index b8298134a3aa10d91d9e6e0dfa3b229a20b62e94..ab2b8a25d95887d198dc568536fec81c3d40b461 100644
--- a/Environment/Convenience.h
+++ b/Environment/Convenience.h
@@ -48,6 +48,7 @@ namespace corsika::environment {
 
     void addLinearLayer(units::si::LengthType, units::si::LengthType);
 
+    void assemble(Environment<IMediumModel>&);
     Environment<IMediumModel> assemble();
   };
 
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
index 71f8e27d917e3e34fc5d8c2fc663d5142d07aa59..fdbde0f4e1d3257ba73768ea1e4d123e6ee6fbbc 100644
--- a/Environment/SlidingPlanarExponential.h
+++ b/Environment/SlidingPlanarExponential.h
@@ -38,36 +38,38 @@ namespace corsika::environment {
   template <class T>
   class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
                                    public T {
-    NuclearComposition const fNuclComp;
+    NuclearComposition const nuclComp_;
+    units::si::LengthType const referenceHeight_;
 
     using Base = BaseExponential<SlidingPlanarExponential<T>>;
 
   public:
-    SlidingPlanarExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
-                             units::si::LengthType vLambda, NuclearComposition vNuclComp)
-        : Base(vP0, vRho, vLambda)
-        , fNuclComp(vNuclComp) {}
+    SlidingPlanarExponential(geometry::Point const& p0, units::si::MassDensityType rho0,
+                             units::si::LengthType lambda, NuclearComposition nuclComp, units::si::LengthType referenceHeight = units::si::LengthType::zero())
+        : Base(p0, rho0, lambda)
+        , nuclComp_(nuclComp),
+        referenceHeight_(referenceHeight) {}
 
     units::si::MassDensityType GetMassDensity(
-        geometry::Point const& vP) const override {
-      auto const height = (vP - Base::fP0).norm();
+        geometry::Point const& p) const override {
+      auto const height = (p - Base::fP0).norm() - referenceHeight_;
       return Base::fRho0 * exp(Base::fInvLambda * height);
     }
 
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+    NuclearComposition const& GetNuclearComposition() const override { return nuclComp_; }
 
     units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::LengthType vL) const override {
-      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
-      return Base::IntegratedGrammage(vLine, vL, axis);
+        geometry::Trajectory<geometry::Line> const& line,
+        units::si::LengthType l) const override {
+      auto const axis = (line.GetR0() - Base::fP0).normalized();
+      return Base::IntegratedGrammage(line, l, axis);
     }
 
     units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::GrammageType vGrammage) const override {
-      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
-      return Base::ArclengthFromGrammage(vLine, vGrammage, axis);
+        geometry::Trajectory<geometry::Line> const& line,
+        units::si::GrammageType grammage) const override {
+      auto const axis = (line.GetR0() - Base::fP0).normalized();
+      return Base::ArclengthFromGrammage(line, grammage, axis);
     }
   };