fixed numerical issues in SlidingPlanarAtmosphere

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_;
 }

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();
   };
 

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);
     }
   };