From aa1c3660719c5ae9260455ae4d7bb468e645a018 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 20 Jan 2020 13:56:03 +0100
Subject: [PATCH] fixed numerical issues in SlidingPlanarAtmosphere
---
Environment/Convenience.cpp | 19 ++++++++------
Environment/Convenience.h | 1 +
Environment/SlidingPlanarExponential.h | 34 ++++++++++++++------------
3 files changed, 31 insertions(+), 23 deletions(-)
diff --git a/Environment/Convenience.cpp b/Environment/Convenience.cpp
index fb2361c03..f3b052d6d 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 b8298134a..ab2b8a25d 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 71f8e27d9..fdbde0f4e 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);
}
};
--
GitLab