From b7deecfabe0b063b6e511f9d7bf72fd7f1fb360b Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 15 Jun 2021 18:14:35 +0200
Subject: [PATCH] some density modelling refactoring
---
corsika/detail/media/BaseExponential.inl | 47 ++++++++-----
corsika/detail/media/FlatExponential.inl | 14 ++--
corsika/detail/media/HomogeneousMedium.inl | 6 +-
corsika/detail/media/InhomogeneousMedium.inl | 4 +-
.../LayeredSphericalAtmosphereBuilder.inl | 35 ++++++++++
.../media/LinearApproximationIntegrator.inl | 6 +-
.../detail/media/SlidingPlanarExponential.inl | 14 ++--
.../modules/energy_loss/BetheBlochPDG.inl | 3 +-
.../modules/proposal/ContinuousProcess.inl | 3 +-
corsika/media/BaseExponential.hpp | 21 +++---
corsika/media/FlatExponential.hpp | 3 +-
corsika/media/HomogeneousMedium.hpp | 3 +-
corsika/media/IMediumModel.hpp | 24 +++++--
corsika/media/InhomogeneousMedium.hpp | 3 +-
.../LayeredSphericalAtmosphereBuilder.hpp | 4 ++
.../media/LinearApproximationIntegrator.hpp | 2 +-
corsika/media/SlidingPlanarExponential.hpp | 3 +-
tests/media/testEnvironment.cpp | 69 +++++++++++--------
tests/media/testMagneticField.cpp | 5 --
tests/media/testMedium.cpp | 9 ++-
tests/media/testRefractiveIndex.cpp | 29 ++++----
21 files changed, 193 insertions(+), 114 deletions(-)
diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl
index 3d7790c29..dbd418041 100644
--- a/corsika/detail/media/BaseExponential.inl
+++ b/corsika/detail/media/BaseExponential.inl
@@ -14,23 +14,43 @@
namespace corsika {
+ template <typename TDerived>
+ inline BaseExponential<TDerived>::BaseExponential(Point const& point,
+ LengthType const referenceHeight,
+ MassDensityType rho0,
+ LengthType lambda)
+ : rho0_(rho0)
+ , lambda_(lambda)
+ , invLambda_(1 / lambda)
+ , point_(point)
+ , referenceHeight_(referenceHeight) {}
+
template <typename TDerived>
inline auto const& BaseExponential<TDerived>::getImplementation() const {
return *static_cast<TDerived const*>(this);
}
+ template <typename TDerived>
+ inline MassDensityType BaseExponential<TDerived>::getMassDensity(
+ LengthType const height) const {
+ return rho0_ * exp(invLambda_ * (height - referenceHeight_));
+ }
+
template <typename TDerived>
inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
- BaseTrajectory const& traj, LengthType vL, DirectionVector const& axis) const {
- if (vL == LengthType::zero()) { return GrammageType::zero(); }
+ BaseTrajectory const& traj, DirectionVector const& axis) const {
+ LengthType const length = traj.getLength();
+ if (length == LengthType::zero()) { return GrammageType::zero(); }
- auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
- auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+ // this corresponds to height:
+ double const uDotA = traj.getDirection(0).dot(axis).magnitude();
+ MassDensityType const rhoStart =
+ getImplementation().getMassDensity(traj.getPosition(0));
if (uDotA == 0) {
- return vL * rhoStart;
+ return length * rhoStart;
} else {
- return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1);
+ return rhoStart * (lambda_ / uDotA) * (exp(uDotA * length * invLambda_) - 1);
}
}
@@ -38,8 +58,10 @@ namespace corsika {
inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
BaseTrajectory const& traj, GrammageType grammage,
DirectionVector const& axis) const {
- auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
- auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
+ // this corresponds to height:
+ double const uDotA = traj.getDirection(0).dot(axis).magnitude();
+ MassDensityType const rhoStart =
+ getImplementation().getMassDensity(traj.getPosition(0));
if (uDotA == 0) {
return grammage / rhoStart;
@@ -54,13 +76,4 @@ namespace corsika {
}
}
- template <typename TDerived>
- inline BaseExponential<TDerived>::BaseExponential(Point const& point,
- MassDensityType rho0,
- LengthType lambda)
- : rho0_(rho0)
- , lambda_(lambda)
- , invLambda_(1 / lambda)
- , point_(point) {}
-
} // namespace corsika
diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl
index b64b2ded2..f665dda99 100644
--- a/corsika/detail/media/FlatExponential.inl
+++ b/corsika/detail/media/FlatExponential.inl
@@ -18,19 +18,17 @@ namespace corsika {
template <typename T>
inline FlatExponential<T>::FlatExponential(Point const& point,
- Vector<dimensionless_d> const& axis,
+ DirectionVector const& axis,
MassDensityType rho, LengthType lambda,
NuclearComposition const& nuclComp)
- : BaseExponential<FlatExponential<T>>(point, rho, lambda)
+ : BaseExponential<FlatExponential<T>>(point, 0_m, rho, lambda)
, axis_(axis)
, nuclComp_(nuclComp) {}
template <typename T>
inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
- return BaseExponential<FlatExponential<T>>::getRho0() *
- exp(BaseExponential<FlatExponential<T>>::getInvLambda() *
- (point - BaseExponential<FlatExponential<T>>::getAnchorPoint())
- .dot(axis_));
+ return BaseExponential<FlatExponential<T>>::getMassDensity(
+ (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm());
}
template <typename T>
@@ -40,8 +38,8 @@ namespace corsika {
template <typename T>
inline GrammageType FlatExponential<T>::getIntegratedGrammage(
- BaseTrajectory const& line, LengthType to) const {
- return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
+ BaseTrajectory const& line) const {
+ return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
}
template <typename T>
diff --git a/corsika/detail/media/HomogeneousMedium.inl b/corsika/detail/media/HomogeneousMedium.inl
index 368073c9f..032209bd3 100644
--- a/corsika/detail/media/HomogeneousMedium.inl
+++ b/corsika/detail/media/HomogeneousMedium.inl
@@ -32,9 +32,9 @@ namespace corsika {
}
template <typename T>
- inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(BaseTrajectory const&,
- LengthType to) const {
- return to * density_;
+ inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
+ BaseTrajectory const& track) const {
+ return track.getLength() * density_;
}
template <typename T>
diff --git a/corsika/detail/media/InhomogeneousMedium.inl b/corsika/detail/media/InhomogeneousMedium.inl
index 8706765ad..f53f2c3be 100644
--- a/corsika/detail/media/InhomogeneousMedium.inl
+++ b/corsika/detail/media/InhomogeneousMedium.inl
@@ -36,8 +36,8 @@ namespace corsika {
template <typename T, typename TDensityFunction>
inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
- BaseTrajectory const& line, LengthType to) const {
- return densityFunction_.getIntegrateGrammage(line, to);
+ BaseTrajectory const& line) const {
+ return densityFunction_.getIntegrateGrammage(line);
}
template <typename T, typename TDensityFunction>
diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index 74f00f101..791d57972 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
@@ -13,6 +13,7 @@
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
+//#include <corsika/media/SlidingPlanarTabular.hpp>
namespace corsika {
@@ -103,6 +104,40 @@ namespace corsika {
layers_.push(std::move(node));
}
+ template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
+ typename... TModelArgs>
+ inline void
+ LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
+ addTabularLayer(GrammageType b,
+ std::function<MassDensityType(LengthType)> const& funcRho,
+ LengthType upperBoundary) {
+
+ auto const radius = earthRadius_ + upperBoundary;
+ checkRadius(radius);
+ previousRadius_ = radius;
+
+ auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+ std::make_unique<Sphere>(center_, radius));
+ /*
+ if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+ // helper lambda in which the last 5 arguments to make_shared<...> are bound
+ auto lastBound = [&](auto... argPack) {
+ return std::make_shared<
+ TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
+ argPack..., center_, rho0, -c, *composition_, earthRadius_);
+ };
+
+ // now unpack the additional arguments
+ auto model = std::apply(lastBound, additionalModelArgs_);
+ node->setModelProperties(std::move(model));
+ } else {
+ node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
+ center_, rho0, -c, *composition_, earthRadius_);
+ }
+ */
+ layers_.push(std::move(node));
+ }
+
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
diff --git a/corsika/detail/media/LinearApproximationIntegrator.inl b/corsika/detail/media/LinearApproximationIntegrator.inl
index ced2dee38..8ee60d984 100644
--- a/corsika/detail/media/LinearApproximationIntegrator.inl
+++ b/corsika/detail/media/LinearApproximationIntegrator.inl
@@ -19,10 +19,14 @@ namespace corsika {
template <typename TDerived>
inline auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
- BaseTrajectory const& line, LengthType length) const {
+ BaseTrajectory const& line) const {
+ LengthType const length = line.getLength();
auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
line.getDirection(0));
+ CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1,
+ line.getPosition(0), line.getDirection(0),
+ (c0 + 0.5 * c1 * length) * length);
return (c0 + 0.5 * c1 * length) * length;
}
diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl
index 61911772a..80fc9bf6b 100644
--- a/corsika/detail/media/SlidingPlanarExponential.inl
+++ b/corsika/detail/media/SlidingPlanarExponential.inl
@@ -16,19 +16,17 @@ namespace corsika {
inline SlidingPlanarExponential<T>::SlidingPlanarExponential(
Point const& p0, MassDensityType rho0, LengthType lambda,
NuclearComposition const& nuclComp, LengthType referenceHeight)
- : BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
+ : BaseExponential<SlidingPlanarExponential<T>>(p0, referenceHeight, rho0, lambda)
, nuclComp_(nuclComp)
, referenceHeight_(referenceHeight) {}
template <typename T>
inline MassDensityType SlidingPlanarExponential<T>::getMassDensity(
Point const& point) const {
- auto const height =
+ auto const heightFull =
(point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
- .getNorm() -
- referenceHeight_;
- return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
- exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * height);
+ .getNorm();
+ return BaseExponential<SlidingPlanarExponential<T>>::getMassDensity(heightFull);
}
template <typename T>
@@ -39,11 +37,11 @@ namespace corsika {
template <typename T>
inline GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
- BaseTrajectory const& traj, LengthType l) const {
+ BaseTrajectory const& traj) const {
auto const axis = (traj.getPosition(0) -
BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
.normalized();
- return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, l,
+ return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj,
axis);
}
diff --git a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
index 7772f2660..537a3d85f 100644
--- a/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+++ b/corsika/detail/modules/energy_loss/BetheBlochPDG.inl
@@ -154,8 +154,7 @@ namespace corsika {
if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
- GrammageType const dX =
- p.getNode()->getModelProperties().getIntegratedGrammage(t, t.getLength());
+ GrammageType const dX = p.getNode()->getModelProperties().getIntegratedGrammage(t);
CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
p.getChargeNumber(), dX / 1_g * square(1_cm));
HEPEnergyType dE = getTotalEnergyLoss(p, dX);
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index aa50f5866..30f0b10a1 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -97,8 +97,7 @@ namespace corsika::proposal {
if (vT.getLength() == 0_m) return ProcessReturn::Ok;
// calculate passed grammage
- auto dX =
- vP.getNode()->getModelProperties().getIntegratedGrammage(vT, vT.getLength());
+ auto dX = vP.getNode()->getModelProperties().getIntegratedGrammage(vT);
// get or build corresponding track integral calculator and solve the
// integral
diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp
index 66f907943..27a267df6 100644
--- a/corsika/media/BaseExponential.hpp
+++ b/corsika/media/BaseExponential.hpp
@@ -23,9 +23,20 @@ namespace corsika {
*/
template <typename TDerived>
class BaseExponential {
+
+ public:
+ BaseExponential(Point const& point, LengthType const referenceHeight,
+ MassDensityType rho0, LengthType lambda);
+
+ Point const& getAnchorPoint() const { return point_; }
+ MassDensityType getRho0() const { return rho0_; }
+ InverseLengthType getInvLambda() const { return invLambda_; }
+
protected:
auto const& getImplementation() const;
+ MassDensityType getMassDensity(LengthType const height) const;
+
// clang-format off
/**
* For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
@@ -40,7 +51,7 @@ namespace corsika {
* \f]
*/
// clang-format on
- GrammageType getIntegratedGrammage(BaseTrajectory const& line, LengthType vL,
+ GrammageType getIntegratedGrammage(BaseTrajectory const& line,
DirectionVector const& axis) const;
// clang-format off
@@ -64,18 +75,12 @@ namespace corsika {
LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage,
DirectionVector const& axis) const;
- public:
- BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda);
-
- Point const& getAnchorPoint() const { return point_; }
- MassDensityType getRho0() const { return rho0_; }
- InverseLengthType getInvLambda() const { return invLambda_; }
-
private:
MassDensityType const rho0_;
LengthType const lambda_;
InverseLengthType const invLambda_;
Point const point_;
+ LengthType const referenceHeight_;
}; // class BaseExponential
diff --git a/corsika/media/FlatExponential.hpp b/corsika/media/FlatExponential.hpp
index 71767d705..6f64295e8 100644
--- a/corsika/media/FlatExponential.hpp
+++ b/corsika/media/FlatExponential.hpp
@@ -41,8 +41,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(BaseTrajectory const& line,
- LengthType to) const override;
+ GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& line,
GrammageType grammage) const override;
diff --git a/corsika/media/HomogeneousMedium.hpp b/corsika/media/HomogeneousMedium.hpp
index 95a9fb0ae..16331957d 100644
--- a/corsika/media/HomogeneousMedium.hpp
+++ b/corsika/media/HomogeneousMedium.hpp
@@ -30,8 +30,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(BaseTrajectory const&,
- LengthType to) const override;
+ GrammageType getIntegratedGrammage(BaseTrajectory const&) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const&,
GrammageType grammage) const override;
diff --git a/corsika/media/IMediumModel.hpp b/corsika/media/IMediumModel.hpp
index 87e702c5d..e3e030392 100644
--- a/corsika/media/IMediumModel.hpp
+++ b/corsika/media/IMediumModel.hpp
@@ -21,11 +21,25 @@ namespace corsika {
virtual MassDensityType getMassDensity(Point const&) const = 0;
- // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
- // approach; for now, only lines are supported
- virtual GrammageType getIntegratedGrammage(BaseTrajectory const&,
- LengthType) const = 0;
-
+ /**
+ * Integrate the matter density along trajectory.
+ *
+ * @return GrammageType as integrated matter density along the BaseTrajectory
+ *
+ * @todo think about the mixin inheritance of the trajectory vs the BaseTrajectory
+ * approach; for now, only lines are supported (?).
+ */
+ virtual GrammageType getIntegratedGrammage(BaseTrajectory const&) const = 0;
+
+ /**
+ * Calculates the length along the trajectory.
+ *
+ * The length along the trajectory is determined at which the integrated matter
+ * density is reached. If the specified matter density cannot be reached (is too
+ * large) the result becomes meaningless and could be "infinity" (discuss this).
+ *
+ * @return LengthType the length corresponding to grammage.
+ */
virtual LengthType getArclengthFromGrammage(BaseTrajectory const&,
GrammageType) const = 0;
diff --git a/corsika/media/InhomogeneousMedium.hpp b/corsika/media/InhomogeneousMedium.hpp
index a257bc508..66cf00f73 100644
--- a/corsika/media/InhomogeneousMedium.hpp
+++ b/corsika/media/InhomogeneousMedium.hpp
@@ -32,8 +32,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(BaseTrajectory const& line,
- LengthType to) const override;
+ GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& pLine,
GrammageType grammage) const override;
diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
index 5de160d13..41ad15801 100644
--- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+++ b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
@@ -78,6 +78,10 @@ namespace corsika {
void addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary);
void addLinearLayer(LengthType c, LengthType upperBoundary);
+ void addTabularLayer(GrammageType b,
+ std::function<MassDensityType(LengthType)> const& funcRho,
+ LengthType upperBoundary);
+
int getSize() const { return layers_.size(); }
void assemble(Environment<TMediumInterface>& env);
diff --git a/corsika/media/LinearApproximationIntegrator.hpp b/corsika/media/LinearApproximationIntegrator.hpp
index 0e2eb0a29..ef442ebdb 100644
--- a/corsika/media/LinearApproximationIntegrator.hpp
+++ b/corsika/media/LinearApproximationIntegrator.hpp
@@ -20,7 +20,7 @@ namespace corsika {
auto const& getImplementation() const;
public:
- auto getIntegrateGrammage(BaseTrajectory const& line, LengthType length) const;
+ auto getIntegrateGrammage(BaseTrajectory const& line) const;
auto getArclengthFromGrammage(BaseTrajectory const& line,
GrammageType grammage) const;
diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp
index a06c293a4..63f8238ba 100644
--- a/corsika/media/SlidingPlanarExponential.hpp
+++ b/corsika/media/SlidingPlanarExponential.hpp
@@ -46,8 +46,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
- GrammageType getIntegratedGrammage(BaseTrajectory const& line,
- LengthType l) const override;
+ GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& line,
GrammageType grammage) const override;
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index b39fbb357..88e143b7b 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -97,33 +97,35 @@ TEST_CASE("FlatExponential") {
auto const rho0 = 1_g / static_pow<3>(1_cm);
FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
protonComposition);
- auto const tEnd = 5_s;
+ SpeedType const speed = 20_m / second;
+ LengthType const length = 2_m;
+ TimeType const tEnd = length / speed;
SECTION("horizontal") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {20_cm / second, 0_m / second, 0_m / second}));
+ gCS, {speed, 0_m / second, 0_m / second}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
- CHECK((medium.getIntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
- CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
+ CHECK((medium.getIntegratedGrammage(trajectory) / (rho0 * length)) == Approx(1));
+ CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * length) / length) ==
+ Approx(1));
}
SECTION("vertical") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {0_m / second, 0_m / second, 5_m / second}));
+ gCS, {0_m / second, 0_m / second, speed}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
- LengthType const length = 2 * lambda;
GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
- CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1));
+ CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
SECTION("escape grammage") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {0_m / second, 0_m / second, -5_m / second}));
+ gCS, {SpeedType::zero(), SpeedType::zero(), -speed}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
GrammageType const escapeGrammage = rho0 * lambda;
@@ -134,15 +136,15 @@ TEST_CASE("FlatExponential") {
}
SECTION("inclined") {
- Line const line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {0_m / second, 5_m / second, 5_m / second}));
+ Line const line(gOrigin,
+ Vector<SpeedType::dimension_type>(
+ gCS, {0_m / second, speed / sqrt(2.), speed / sqrt(2.)}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
double const cosTheta = M_SQRT1_2;
- LengthType const length = 2 * lambda;
GrammageType const exact =
rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
- CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1));
+ CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
}
@@ -179,16 +181,16 @@ TEST_CASE("SlidingPlanarExponential") {
CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
flat.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
- CHECK(medium.getIntegratedGrammage(trajectory, 2_m).magnitude() ==
- flat.getIntegratedGrammage(trajectory, 2_m).magnitude());
+ CHECK(medium.getIntegratedGrammage(trajectory).magnitude() ==
+ flat.getIntegratedGrammage(trajectory).magnitude());
CHECK(medium.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
flat.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
}
}
-auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
+MassDensityType constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
-struct Exponential {
+struct ExponentialTest {
auto operator()(Point const& p) const {
return exp(p.getCoordinates()[0] / 1_m) * rho0;
}
@@ -213,41 +215,50 @@ TEST_CASE("InhomogeneousMedium") {
Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
+ SpeedType const speed = 20_m / second;
Line line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {20_m / second, 0_m / second, 0_m / second}));
+ gCS, {speed, SpeedType::zero(), SpeedType::zero()}));
- auto const tEnd = 5_s;
+ // the tested LinearApproximationIntegrator really does a single step only. It is very
+ // poor for exponentials with a bit larger step-width.
+ TimeType const tEnd = 0.001_s;
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
- Exponential const e;
- DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
+ ExponentialTest const expTest;
+ DensityFunction<ExponentialTest, LinearApproximationIntegrator> const rho(expTest);
SECTION("DensityFunction") {
- CHECK(e.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+ CHECK(expTest.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
Approx(1));
- CHECK(rho.evaluateAt(gOrigin) == e(gOrigin));
+ CHECK(rho.evaluateAt(gOrigin) == expTest(gOrigin));
}
auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); };
- auto constexpr l = 15_cm;
+ LengthType constexpr length = tEnd * speed;
NuclearComposition const composition{{Code::Proton}, {1.f}};
InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
+ CORSIKA_LOG_INFO("test={} l={} {} {}", rho.getIntegrateGrammage(trajectory), length,
+ exactGrammage(length), 1_m * rho0 * (exp(length / 1_m) - 1));
+
SECTION("Integration") {
- CHECK(rho.getIntegrateGrammage(trajectory, l) / exactGrammage(l) ==
+ CORSIKA_LOG_INFO("test={} {} {}", rho.getIntegrateGrammage(trajectory),
+ exactGrammage(length),
+ rho.getIntegrateGrammage(trajectory) / exactGrammage(length));
+ CHECK(rho.getIntegrateGrammage(trajectory) / exactGrammage(length) ==
Approx(1).epsilon(1e-2));
- CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(l)) /
- exactLength(exactGrammage(l)) ==
+ CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(length)) /
+ exactLength(exactGrammage(length)) ==
Approx(1).epsilon(1e-2));
CHECK(rho.getMaximumLength(trajectory, 1e-2) >
- l); // todo: write reasonable test when implementation is working
+ length); // todo: write reasonable test when implementation is working
- CHECK(rho.getIntegrateGrammage(trajectory, l) ==
- inhMedium.getIntegratedGrammage(trajectory, l));
+ CHECK(rho.getIntegrateGrammage(trajectory) ==
+ inhMedium.getIntegratedGrammage(trajectory));
CHECK(rho.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
inhMedium.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
CHECK(inhMedium.getNuclearComposition() == composition);
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index b3ccf4d75..6e9f347b8 100644
--- a/tests/media/testMagneticField.cpp
+++ b/tests/media/testMagneticField.cpp
@@ -25,7 +25,6 @@ using namespace corsika;
TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
@@ -85,8 +84,4 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
// and the associated trajectory
setup::Trajectory const track =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
-
- // and check the integrated grammage
- CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
- CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
}
diff --git a/tests/media/testMedium.cpp b/tests/media/testMedium.cpp
index a10348717..bb5a9f0e8 100644
--- a/tests/media/testMedium.cpp
+++ b/tests/media/testMedium.cpp
@@ -92,18 +92,21 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.getNuclearComposition();
+ SpeedType const speed = 1_m / second;
+
// create a line of length 1 m
- Line const line(gOrigin,
- VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second}));
+ Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second}));
// the end time of our line
auto const tEnd = 1_s;
+ LengthType const length = tEnd * speed;
+
// and the associated trajectory
setup::Trajectory const trajectory =
corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
- CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.getIntegratedGrammage(trajectory) / (density * length)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index acef32c4f..62422f7a8 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -76,19 +76,22 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.getNuclearComposition();
+ SpeedType const speed = 1_m / second;
+
// create a line of length 1 m
- Line const line(gOrigin,
- VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second}));
+ Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second}));
// the end time of our line
auto const tEnd = 1_s;
+ LengthType const length = tEnd * speed;
+
// and the associated trajectory
setup::Trajectory const track =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
- CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.getIntegratedGrammage(track) / (density * length)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
}
@@ -148,24 +151,26 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
REQUIRE(density == medium__.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium__.getNuclearComposition();
- // create a line of length 1 m
- Line const line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {1_m / second, 0_m / second, 0_m / second}));
+ SpeedType const velocity = 1_m / second;
// the end time of our line
- auto const tEnd = 1_s;
+ TimeType const tEnd = 1_s;
+
+ LengthType const length = tEnd * velocity;
+
+ // create a line of length 1 m
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {velocity, 0_m / second, 0_m / second}));
// and the associated trajectory
setup::Trajectory const track =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
- // // and the associated trajectory
- // Trajectory<Line> const trajectory(line, tEnd);
// and check the integrated grammage
- REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium.getIntegratedGrammage(track) / (density * length)) == Approx(1));
REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
- REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium_.getIntegratedGrammage(track) / (density * length)) == Approx(1));
REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
- REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium__.getIntegratedGrammage(track) / (density * length)) == Approx(1));
REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
}
--
GitLab