From f90cdda1df48e5575a9bb4e9d54c0c8d3bc0efd5 Mon Sep 17 00:00:00 2001
From: Alan Coleman <alanc@udel.edu>
Date: Tue, 31 Jan 2023 10:56:47 +0100
Subject: [PATCH] Fix issue 550, added unit tests to catch this
---
corsika/detail/media/FlatExponential.inl | 2 +-
tests/media/testEnvironment.cpp | 12 ++++++++++++
2 files changed, 13 insertions(+), 1 deletion(-)
diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl
index ea0fe5b0e..7604352b5 100644
--- a/corsika/detail/media/FlatExponential.inl
+++ b/corsika/detail/media/FlatExponential.inl
@@ -29,7 +29,7 @@ namespace corsika {
template <typename T>
inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
return BaseExponential<FlatExponential<T>>::getMassDensity(
- (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm());
+ (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_));
}
template <typename T>
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 81c1ac757..28368aa96 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -101,11 +101,19 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2_m;
TimeType const tEnd = length / speed;
+ CHECK((medium.getMassDensity(gOrigin)) == rho0);
CHECK(medium.getNuclearComposition().getFractions() == std::vector<double>{1.});
CHECK(medium.getNuclearComposition().getComponents() ==
std::vector<Code>{Code::Proton});
SECTION("horizontal") {
+ //Check that not moving along axis does not change density
+ CHECK(medium.getMassDensity(Point(gCS, 1_m, 0_m, 0_m)) == rho0);
+ CHECK(medium.getMassDensity(Point(gCS, 1_m, 0_m, 0_m)) == rho0);
+ CHECK(medium.getMassDensity(Point(gCS, -1_m, 0_m, 0_m)) == rho0);
+ CHECK(medium.getMassDensity(Point(gCS, 0_m, 1_m, 0_m)) == rho0);
+ CHECK(medium.getMassDensity(Point(gCS, 0_m,-1_m, 0_m)) == rho0);
+
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {speed, 0_m / second, 0_m / second}));
setup::Trajectory const trajectory =
@@ -117,6 +125,10 @@ TEST_CASE("FlatExponential") {
}
SECTION("vertical") {
+ //Moving along axis does change density
+ CHECK(medium.getMassDensity(Point(gCS, 0_m, 0_m, 1_m)) > rho0);
+ CHECK(medium.getMassDensity(Point(gCS, 0_m, 0_m, -1_m)) < rho0);
+
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, speed}));
setup::Trajectory const trajectory =
--
GitLab