From 68bbace4d35de56ad703ae4762cc4fd51d64ead1 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Sat, 23 Feb 2019 17:19:31 +0100
Subject: [PATCH] fixed bug concerning escape grammage in FlatExponential
---
Environment/FlatExponential.h | 13 ++++++++++---
Environment/testEnvironment.cc | 15 +++++++++++++--
2 files changed, 23 insertions(+), 5 deletions(-)
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
index 7079a22bf..19e42c725 100644
--- a/Environment/FlatExponential.h
+++ b/Environment/FlatExponential.h
@@ -1,6 +1,5 @@
-
/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
@@ -21,6 +20,7 @@
#include <corsika/units/PhysicalUnits.h>
#include <cassert>
+#include <limits>
/**
*
@@ -76,7 +76,14 @@ namespace corsika::environment {
if (vDotA == 0) {
return pGrammage / GetMassDensity(line.GetR0());
} else {
- return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+ 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;
+ }
}
}
};
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 5beab6587..206849e7e 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -75,7 +75,19 @@ TEST_CASE("FlatExponential") {
REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
- SECTION("vertical") {
+ SECTION("escape grammage") {
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {0_m / second, 0_m / second, -5_m / second}));
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ GrammageType const escapeGrammage = rho0 * lambda;
+
+ REQUIRE(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
+ REQUIRE(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
+ std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
+ }
+
+ SECTION("inclined") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 5_m / second, 5_m / second}));
Trajectory<Line> const trajectory(line, tEnd);
@@ -83,7 +95,6 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2 * lambda;
GrammageType const exact =
rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
-
REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
--
GitLab