diff --git a/Environment/BaseExponential.h b/Environment/BaseExponential.h
new file mode 100644
index 0000000000000000000000000000000000000000..6d499b366359e2ca606505dc076050bbe5fe9cbc
--- /dev/null
+++ b/Environment/BaseExponential.h
@@ -0,0 +1,83 @@
+/*
+ * (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 <cassert>
+#include <limits>
+
+/**
+ *
+ */
+
+namespace corsika::environment {
+
+ template <class TDerived>
+ class BaseExponential {
+ protected:
+ corsika::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); }
+
+ corsika::units::si::GrammageType IntegratedGrammage(
+ corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::units::si::LengthType pTo,
+ geometry::Vector<units::si::dimensionless_d> const& axis) const {
+ auto const vDotA = line.NormalizedDirection().dot(axis).magnitude();
+
+ if (vDotA == 0) {
+ return pTo * GetImplementation().GetMassDensity(line.GetR0());
+ } else {
+ return GetImplementation().GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
+ (exp(vDotA * pTo / fLambda) - 1);
+ }
+ }
+
+ corsika::units::si::LengthType ArclengthFromGrammage(
+ corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::units::si::GrammageType pGrammage,
+ geometry::Vector<units::si::dimensionless_d> const& axis) const {
+ auto const vDotA = line.NormalizedDirection().dot(axis).magnitude();
+
+ if (vDotA == 0) {
+ return pGrammage / GetImplementation().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;
+ }
+ }
+ }
+
+ public:
+ BaseExponential(geometry::Point const& p0, units::si::MassDensityType rho,
+ units::si::LengthType lambda)
+ : fRho0(rho)
+ , fLambda(lambda)
+ , fInvLambda(1 / lambda)
+ , fP0(p0) {}
+ };
+
+} // namespace corsika::environment
+#endif
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 0190d6dc8dadd18dc2b1bd6d0345be23545a3e47..13fd4f1e736d70a00b8380fdb6de6d906349b4a5 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -10,6 +10,7 @@ set (
DensityFunction.h
Environment.h
NameModel.h
+ BaseExponential.h
FlatExponential.h
SlidingPlanarExponential.h
)
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index 19e42c725c4b5204dc9c6fc166f9be00f44bdb7a..97fba098e4ee0f09698d24c0b90bbb2829df414f 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -11,12 +11,12 @@
#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>
@@ -29,64 +29,39 @@
namespace corsika::environment {
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)
+ NuclearComposition nuclComp)
+ : Base(p0, rho, lambda)
, fAxis(axis)
- , fP0(p0) {}
+ , fNuclComp(nuclComp) {}
corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const& p) const override {
- return fRho0 * exp(fInvLambda * (p - fP0).dot(fAxis));
+ return Base::fRho0 * exp(Base::fInvLambda * (p - 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();
-
- if (vDotA == 0) {
- return pTo * GetMassDensity(line.GetR0());
- } else {
- return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
- (exp(vDotA * pTo / fLambda) - 1);
- }
+ return Base::IntegratedGrammage(line, pTo, 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;
- }
- }
+ return Base::ArclengthFromGrammage(line, pGrammage, fAxis);
}
};
-
} // namespace corsika::environment
#endif
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
index 3dc05917361eaf39d7fbffd1c0b5ecb20475a3ca..c04c084ec5775f4307004a907c065c9761bb57b5 100644
--- a/Environment/SlidingPlanarExponential.h
+++ b/Environment/SlidingPlanarExponential.h
@@ -11,15 +11,16 @@
#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/environment/FlatExponential.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 <fenv.h>
#include <cassert>
#include <limits>
@@ -30,30 +31,22 @@
namespace corsika::environment {
template <class T>
- class SlidingPlanarExponential : public T {
- corsika::units::si::MassDensityType const fRho0;
- units::si::LengthType const fLambda;
- units::si::InverseLengthType const fInvLambda;
+ class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
+ public T {
NuclearComposition const fNuclComp;
- geometry::Vector<units::si::dimensionless_d> const fAxis;
- geometry::Point const fP0;
+
+ using Base = BaseExponential<SlidingPlanarExponential<T>>;
public:
- SlidingPlanarExponential(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) {}
+ SlidingPlanarExponential(geometry::Point const& p0, units::si::MassDensityType rho,
+ units::si::LengthType lambda, NuclearComposition nuclComp)
+ : Base(p0, rho, lambda)
+ , fNuclComp(nuclComp) {}
corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const& p) const override {
- auto const height = (p - fP0).norm();
- return fRho0 * exp(fInvLambda * height);
+ auto const height = (p - Base::fP0).norm();
+ return Base::fRho0 * exp(Base::fInvLambda * height);
}
NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
@@ -61,36 +54,17 @@ namespace corsika::environment {
corsika::units::si::GrammageType IntegratedGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
corsika::units::si::LengthType pTo) const override {
-
- auto const axis = (fP0 - line.GetR0()).normalized();
- auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
-
- if (vDotA == 0) {
- return pTo * GetMassDensity(line.GetR0());
- } else {
- return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
- (exp(vDotA * pTo / fLambda) - 1);
- }
+ feenableexcept(FE_INVALID);
+ auto const axis = (Base::fP0 - line.GetR0()).normalized();
+ return Base::IntegratedGrammage(line, pTo, axis);
}
corsika::units::si::LengthType ArclengthFromGrammage(
corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
corsika::units::si::GrammageType pGrammage) const override {
- auto const axis = (fP0 - line.GetR0()).normalized();
- 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;
- }
- }
+ feenableexcept(FE_INVALID);
+ auto const axis = (Base::fP0 - line.GetR0()).normalized();
+ return Base::ArclengthFromGrammage(line, pGrammage, axis);
}
};
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 {