diff --git a/Environment/BaseExponential.h b/Environment/BaseExponential.h
new file mode 100644
index 0000000000000000000000000000000000000000..15585b6af48cdb33f4484d35da228ece4c8831f0
--- /dev/null
+++ b/Environment/BaseExponential.h
@@ -0,0 +1,114 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_Environment_BaseExponential_h_
+#define _include_Environment_BaseExponential_h_
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <limits>
+
+namespace corsika::environment {
+
+ /**
+ * This class provides the grammage/length conversion functionality for
+ * (locally) flat exponential atmospheres.
+ */
+ template <class TDerived>
+ class BaseExponential {
+ protected:
+ units::si::MassDensityType const fRho0;
+ units::si::LengthType const fLambda;
+ units::si::InverseLengthType const fInvLambda;
+ geometry::Point const fP0;
+
+ auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
+
+ // clang-format off
+ /**
+ * For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
+ * direction \f$ \vec{u} \f$ is given by
+ * \f[
+ * X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right)
+ * \f], where \f$ \varrho_0 \f$ is the density at the starting point.
+ *
+ * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
+ * \f[
+ * X = \varrho_0 l;
+ * \f]
+ */
+ // clang-format on
+ units::si::GrammageType IntegratedGrammage(
+ geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
+ geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
+ auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
+ auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+
+ if (uDotA == 0) {
+ return vL * rhoStart;
+ } else {
+ return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1);
+ }
+ }
+
+ // clang-format off
+ /**
+ * For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized)
+ * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by
+ * \f[
+ * l = \begin{cases}
+ * \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 +
+ * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda}
+ * \infty & \text{else,}
+ * \end{cases}
+ * \f] where \f$ \varrho_0 \f$ is the density at the starting point.
+ *
+ * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
+ * \f[
+ * l = \frac{X}{\varrho_0}
+ * \f]
+ */
+ // clang-format on
+ units::si::LengthType ArclengthFromGrammage(
+ geometry::Trajectory<geometry::Line> const& vLine,
+ units::si::GrammageType vGrammage,
+ geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
+ auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
+ auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+
+ if (uDotA == 0) {
+ return vGrammage / rhoStart;
+ } else {
+ auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1;
+ if (logArg > 0) {
+ return fLambda / uDotA * log(logArg);
+ } else {
+ return std::numeric_limits<typename decltype(
+ vGrammage)::value_type>::infinity() *
+ units::si::meter;
+ }
+ }
+ }
+
+ public:
+ BaseExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
+ units::si::LengthType vLambda)
+ : fRho0(vRho)
+ , fLambda(vLambda)
+ , fInvLambda(1 / vLambda)
+ , fP0(vP0) {}
+ };
+
+} // namespace corsika::environment
+#endif
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 0a94c410aeebbf353d225e9e7ea2822d814b1728..13fd4f1e736d70a00b8380fdb6de6d906349b4a5 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -10,7 +10,9 @@ set (
DensityFunction.h
Environment.h
NameModel.h
+ BaseExponential.h
FlatExponential.h
+ SlidingPlanarExponential.h
)
set (
diff --git a/Environment/FlatAtmosphere/FlatAtmosphere.h b/Environment/FlatAtmosphere/FlatAtmosphere.h
deleted file mode 100644
index 6839def48c2f6fddfafd0beb5fd4ca46058ebd09..0000000000000000000000000000000000000000
--- a/Environment/FlatAtmosphere/FlatAtmosphere.h
+++ /dev/null
@@ -1,10 +0,0 @@
-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * 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.
- */
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index 61fa2dced76357946bccfc9fc2951752361c9a42..8428508594f6c4286655103da523d77426a46cdb 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -12,82 +12,60 @@
#ifndef _include_Environment_FlatExponential_h_
#define _include_Environment_FlatExponential_h_
+#include <corsika/environment/BaseExponential.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
-#include <cassert>
-#include <limits>
-
-/**
- *
- */
-
namespace corsika::environment {
+ //clang-format off
+ /**
+ * flat exponential density distribution with
+ * \f[
+ * \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
+ * \vec{a} \right).
+ * \f]
+ * \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
+ * with the scale parameter \f$ \lambda \f$.
+ */
+ //clang-format on
template <class T>
- class FlatExponential : public T {
- corsika::units::si::MassDensityType const fRho0;
- units::si::LengthType const fLambda;
- units::si::InverseLengthType const fInvLambda;
- NuclearComposition const fNuclComp;
+ class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
geometry::Vector<units::si::dimensionless_d> const fAxis;
- geometry::Point const fP0;
+ NuclearComposition const fNuclComp;
+
+ using Base = BaseExponential<FlatExponential<T>>;
public:
- FlatExponential(geometry::Point const& p0,
- geometry::Vector<units::si::dimensionless_d> const& axis,
- units::si::MassDensityType rho, units::si::LengthType lambda,
- NuclearComposition pNuclComp)
- : fRho0(rho)
- , fLambda(lambda)
- , fInvLambda(1 / lambda)
- , fNuclComp(pNuclComp)
- , fAxis(axis)
- , fP0(p0) {}
+ FlatExponential(geometry::Point const& vP0,
+ geometry::Vector<units::si::dimensionless_d> const& vAxis,
+ units::si::MassDensityType vRho, units::si::LengthType vLambda,
+ NuclearComposition vNuclComp)
+ : Base(vP0, vRho, vLambda)
+ , fAxis(vAxis)
+ , fNuclComp(vNuclComp) {}
- corsika::units::si::MassDensityType GetMassDensity(
- corsika::geometry::Point const& p) const override {
- return fRho0 * exp(fInvLambda * (p - fP0).dot(fAxis));
+ units::si::MassDensityType GetMassDensity(geometry::Point const& vP) const override {
+ return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
}
- NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
- corsika::units::si::GrammageType IntegratedGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
- corsika::units::si::LengthType pTo) const override {
- auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
+ NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
- if (vDotA == 0) {
- return pTo * GetMassDensity(line.GetR0());
- } else {
- return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
- (exp(vDotA * pTo / fLambda) - 1);
- }
+ units::si::GrammageType IntegratedGrammage(
+ geometry::Trajectory<geometry::Line> const& vLine,
+ units::si::LengthType vTo) const override {
+ return Base::IntegratedGrammage(vLine, vTo, fAxis);
}
- corsika::units::si::LengthType ArclengthFromGrammage(
- corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
- corsika::units::si::GrammageType pGrammage) const override {
- auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
-
- if (vDotA == 0) {
- return pGrammage / GetMassDensity(line.GetR0());
- } else {
- auto const logArg = pGrammage * vDotA / (fRho0 * fLambda) + 1;
- if (logArg > 0) {
- return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
- } else {
- return std::numeric_limits<typename decltype(
- pGrammage)::value_type>::infinity() *
- corsika::units::si::meter;
- }
- }
+ units::si::LengthType ArclengthFromGrammage(
+ geometry::Trajectory<geometry::Line> const& vLine,
+ units::si::GrammageType vGrammage) const override {
+ return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
}
};
-
} // namespace corsika::environment
#endif
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
new file mode 100644
index 0000000000000000000000000000000000000000..3cffa3afb7b4cb0ce803b4f0ed4edfd1a2392445
--- /dev/null
+++ b/Environment/SlidingPlanarExponential.h
@@ -0,0 +1,74 @@
+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_Environment_SlidingPlanarExponential_h_
+#define _include_Environment_SlidingPlanarExponential_h_
+
+#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::environment {
+
+ // clang-format off
+ /**
+ * The SlidingPlanarExponential models mass density as
+ * \f[
+ * \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \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
+ * from \f$ p_0 \f$ to \f$ r_0 \f$.
+ */
+ //clang-format on
+
+ template <class T>
+ class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
+ public T {
+ NuclearComposition const fNuclComp;
+
+ 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) {}
+
+ units::si::MassDensityType GetMassDensity(
+ geometry::Point const& vP) const override {
+ auto const height = (vP - Base::fP0).norm();
+ return Base::fRho0 * exp(Base::fInvLambda * height);
+ }
+
+ NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+
+ 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);
+ }
+
+ 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);
+ }
+ };
+
+} // namespace corsika::environment
+#endif
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 206849e7eb17b6d33a23fffddb4abefbb16900e8..0821a24271c951456a4a7364e6aad241d6fdc6d3 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -19,6 +19,7 @@
#include <corsika/environment/InhomogeneousMedium.h>
#include <corsika/environment/LinearApproximationIntegrator.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/RootCoordinateSystem.h>
@@ -100,6 +101,41 @@ TEST_CASE("FlatExponential") {
}
}
+TEST_CASE("SlidingPlanarExponential") {
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ LengthType const lambda = 3_m;
+ auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+ auto const tEnd = 5_s;
+
+ SlidingPlanarExponential<IMediumModel> const medium(gOrigin, rho0, lambda,
+ protonComposition);
+
+ SECTION("density") {
+ CHECK(medium.GetMassDensity({gCS, {0_m, 0_m, 3_m}}) /
+ medium.GetMassDensity({gCS, {0_m, 3_m, 0_m}}) ==
+ Approx(1));
+ }
+
+ SECTION("vertical") {
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+ FlatExponential<IMediumModel> const flat(gOrigin, axis, rho0, lambda,
+ protonComposition);
+ Line const line({gCS, {0_m, 0_m, 1_m}},
+ Vector<SpeedType::dimension_type>(
+ gCS, {0_m / second, 0_m / second, 5_m / second}));
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ CHECK(medium.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
+ flat.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
+ CHECK(medium.IntegratedGrammage(trajectory, 2_m).magnitude() ==
+ flat.IntegratedGrammage(trajectory, 2_m).magnitude());
+ CHECK(medium.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
+ flat.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
+ }
+}
+
auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
struct Exponential {