From 5d07a85c6c9744ce1ad5afaa82f9723cb7890d05 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 20 Apr 2021 14:44:27 +0200
Subject: [PATCH] solver update
---
conanfile.txt | 2 +-
corsika/detail/framework/core/Cascade.inl | 61 +++-
corsika/detail/framework/geometry/Helix.inl | 2 +-
.../framework/process/ProcessSequence.inl | 14 +-
.../process/SwitchProcessSequence.inl | 36 +--
.../detail/framework/utility/CubicSolver.inl | 299 ++++++++++++++++++
.../detail/framework/utility/LinearSolver.inl | 36 +++
.../framework/utility/QuadraticSolver.inl | 81 +++++
.../framework/utility/QuarticSolver.inl | 265 +++++++++-------
corsika/detail/modules/ObservationPlane.inl | 9 +-
.../tracking/TrackingLeapFrogCurved.inl | 273 ++++++++--------
.../tracking/TrackingLeapFrogStraight.inl | 25 +-
.../modules/tracking/TrackingStraight.inl | 1 +
corsika/detail/modules/urqmd/UrQMD.inl | 4 +-
.../detail/stack/NuclearStackExtension.inl | 8 +
corsika/framework/core/Cascade.hpp | 5 +-
.../framework/geometry/LeapFrogTrajectory.hpp | 4 +-
corsika/framework/utility/CubicSolver.hpp | 88 ++++++
corsika/framework/utility/LinearSolver.hpp | 24 ++
corsika/framework/utility/QuadraticSolver.hpp | 27 ++
corsika/framework/utility/QuarticSolver.hpp | 69 ++--
corsika/modules/ObservationPlane.hpp | 17 +-
corsika/stack/NuclearStackExtension.hpp | 4 +
corsika/stack/VectorStack.hpp | 6 +
examples/vertical_EAS.cpp | 42 ++-
tests/framework/CMakeLists.txt | 1 +
tests/framework/testCOMBoost.cpp | 3 -
tests/framework/testCascade.cpp | 1 -
tests/framework/testClassTimer.cpp | 1 -
tests/framework/testCombinedStack.cpp | 2 -
tests/framework/testCorsikaFenv.cpp | 1 -
tests/framework/testFourVector.cpp | 1 -
tests/framework/testFunctionTimer.cpp | 1 -
tests/framework/testHelix.cpp | 1 -
tests/framework/testInteractionCounter.cpp | 1 -
tests/framework/testLogging.cpp | 2 +-
tests/framework/testNullModel.cpp | 1 -
tests/framework/testParticles.cpp | 1 -
tests/framework/testRandom.cpp | 2 -
tests/framework/testSaveBoostHistogram.cpp | 1 -
tests/framework/testSecondaryView.cpp | 1 -
tests/framework/testSolver.cpp | 278 ++++++++++++++++
tests/framework/testStackInterface.cpp | 1 -
tests/framework/testUnits.cpp | 1 -
tests/modules/testCONEX.cpp | 2 -
tests/modules/testExecTime.cpp | 1 -
tests/modules/testNullModel.cpp | 1 -
tests/modules/testOnShellCheck.cpp | 1 -
tests/modules/testQGSJetII.cpp | 11 +-
tests/modules/testSibyll.cpp | 2 -
tests/modules/testStackInspector.cpp | 1 -
tests/modules/testTracking.cpp | 14 +-
tests/modules/testUrQMD.cpp | 19 --
53 files changed, 1340 insertions(+), 415 deletions(-)
create mode 100644 corsika/detail/framework/utility/CubicSolver.inl
create mode 100644 corsika/detail/framework/utility/LinearSolver.inl
create mode 100644 corsika/detail/framework/utility/QuadraticSolver.inl
create mode 100644 corsika/framework/utility/CubicSolver.hpp
create mode 100644 corsika/framework/utility/LinearSolver.hpp
create mode 100644 corsika/framework/utility/QuadraticSolver.hpp
create mode 100644 tests/framework/testSolver.cpp
diff --git a/conanfile.txt b/conanfile.txt
index 85c1238b4..b4471ff3d 100644
--- a/conanfile.txt
+++ b/conanfile.txt
@@ -1,5 +1,5 @@
[requires]
-spdlog/1.8.0
+spdlog/1.8.5
catch2/2.13.3
eigen/3.3.8
boost/1.74.0
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 28306c590..f30ec982f 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -60,7 +60,7 @@ namespace corsika {
setNodes();
auto vParticle = stack_.getNextParticle();
TStackView secondaries(vParticle);
- interaction(secondaries);
+ interaction(secondaries, sequence_.getInverseInteractionLength(vParticle));
sequence_.doSecondaries(secondaries);
vParticle.erase(); // primary particle is done
}
@@ -243,15 +243,15 @@ namespace corsika {
[[maybe_unused]] auto projectile = secondaries.getProjectile();
if (distance_interact < distance_decay) {
- interaction(secondaries);
+ interaction(secondaries, total_inv_lambda);
} else {
- decay(secondaries);
- // make sure particle actually did decay if it should have done so
- // \todo this should go to a validation code and not be included here
- if (secondaries.getSize() == 1 &&
- projectile.getPID() == secondaries.getNextParticle().getPID())
- throw std::runtime_error(fmt::format("Particle {} decays into itself!",
- get_name(projectile.getPID())));
+ if (decay(secondaries, total_inv_lifetime) == ProcessReturn::Decayed) {
+ if (secondaries.getSize() == 1 &&
+ projectile.getPID() == secondaries.getNextParticle().getPID()) {
+ throw std::runtime_error(fmt::format("Particle {} decays into itself!",
+ get_name(projectile.getPID())));
+ }
+ }
}
sequence_.doSecondaries(secondaries);
@@ -261,16 +261,29 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::decay(
- TStackView& view) {
+ TStackView& view, InverseTimeType initial_inv_decay_time) {
CORSIKA_LOG_DEBUG("decay");
+
+#ifdef DEBUG
InverseTimeType const actual_decay_time = sequence_.getInverseLifetime(view.parent());
+ if (actual_decay_time > initial_inv_decay_time) {
+ CORSIKA_LOG_WARN(
+ "Decay time decreased during step! This leads to un-physical step length. "
+ "initial_decay_time={}, actual_decay_time={}",
+ 1 / initial_inv_decay_time, 1 / actual_decay_time);
+ }
+#endif
- UniformRealDistribution<InverseTimeType> uniDist(actual_decay_time);
+ // one option is that decay_time is now larger (less
+ // probability for decay) than it was before the step, thus,
+ // no decay might actually occur and is allowed
+
+ UniformRealDistribution<InverseTimeType> uniDist(initial_inv_decay_time);
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectDecay(view, sample_process);
if (returnCode != ProcessReturn::Decayed) {
- CORSIKA_LOG_WARN("Particle did not decay!");
+ CORSIKA_LOG_DEBUG("Particle did not decay!");
}
setEventType(view, history::EventType::Decay);
return returnCode;
@@ -279,18 +292,32 @@ namespace corsika {
template <typename TTracking, typename TProcessList, typename TStack,
typename TStackView>
inline ProcessReturn Cascade<TTracking, TProcessList, TStack, TStackView>::interaction(
- TStackView& view) {
+ TStackView& view, InverseGrammageType initial_inv_int_length) {
+
CORSIKA_LOG_DEBUG("collide");
- InverseGrammageType const current_inv_length =
- sequence_.getInverseInteractionLength(view.parent());
+#ifdef DEBUG
+ InverseGrammageType const actual_inv_length = sequence_.getInverseInteractionLength(
+ view.parent()); // 1/lambda_int after step, -dE/dX etc.
+
+ if (actual_inv_length > initial_inv_int_length) {
+ CORSIKA_LOG_WARN(
+ "Interaction length decreased during step! This leads to un-physical step "
+ "length. "
+ "initial_inv_int_length={}, actual_inv_length={}",
+ 1 / initial_inv_int_length, 1 / actual_inv_length);
+ }
+#endif
- UniformRealDistribution<InverseGrammageType> uniDist(current_inv_length);
+ // one option is that interaction_length is now larger (less
+ // probability for collision) than it was before the step, thus,
+ // no interaction might actually occur and is allowed
+ UniformRealDistribution<InverseGrammageType> uniDist(initial_inv_int_length);
const auto sample_process = uniDist(rng_);
auto const returnCode = sequence_.selectInteraction(view, sample_process);
if (returnCode != ProcessReturn::Interacted) {
- CORSIKA_LOG_WARN("Particle did not interace!");
+ CORSIKA_LOG_DEBUG("Particle did not interact!");
}
setEventType(view, history::EventType::Interaction);
return returnCode;
diff --git a/corsika/detail/framework/geometry/Helix.inl b/corsika/detail/framework/geometry/Helix.inl
index 03c0fd85d..5a029d8e2 100644
--- a/corsika/detail/framework/geometry/Helix.inl
+++ b/corsika/detail/framework/geometry/Helix.inl
@@ -19,7 +19,7 @@ namespace corsika {
inline Point Helix::getPosition(TimeType const t) const {
return r0_ + vPar_ * t +
- (vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
+ (vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
omegaC_;
}
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index 6c9bd7892..7b9ec0dfb 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -221,7 +221,7 @@ namespace corsika {
IndexProcess2>::doStack(TStack& stack) {
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// errors if process1_type is invalid
- if constexpr (std::is_base_of_v<StackProcess<process1_type>, process1_type> ||
+ if constexpr (is_stack_process_v<process1_type> ||
(process1_type::is_process_sequence &&
!process1_type::is_switch_process_sequence)) {
@@ -237,7 +237,7 @@ namespace corsika {
}
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
- if constexpr (std::is_base_of_v<StackProcess<process2_type>, process2_type> ||
+ if constexpr (is_stack_process_v<process2_type> ||
(process2_type::is_process_sequence &&
!process2_type::is_switch_process_sequence)) {
@@ -316,14 +316,14 @@ namespace corsika {
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// errors if process1_type is invalid
- if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> ||
+ if constexpr (is_interaction_process_v<process1_type> ||
process1_type::is_process_sequence) {
tot += A_.getInverseInteractionLength(particle);
}
}
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
- if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> ||
+ if constexpr (is_interaction_process_v<process2_type> ||
process2_type::is_process_sequence) {
tot += B_.getInverseInteractionLength(particle);
}
@@ -351,8 +351,7 @@ namespace corsika {
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != ProcessReturn::Ok) { return ret; }
} else if constexpr (is_process_v<process1_type> &&
- std::is_base_of_v<InteractionProcess<process1_type>,
- process1_type>) {
+ is_interaction_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -377,8 +376,7 @@ namespace corsika {
// if B_ is a process sequence --> check inside
return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
} else if constexpr (is_process_v<process2_type> &&
- std::is_base_of_v<InteractionProcess<process2_type>,
- process2_type>) {
+ is_interaction_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// soon as SecondaryView::parent() is migrated!
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index 2f87ade6d..634395f55 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -36,13 +36,11 @@ namespace corsika {
typename TParticle::node_type const& from,
typename TParticle::node_type const& to) {
if (select_(particle)) {
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
- process1_type> ||
+ if constexpr (is_boundary_process_v<process1_type> ||
process1_type::is_process_sequence) {
// interface checking on TSequence
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process1_type>,
- process1_type>) {
+ if constexpr (is_boundary_process_v<process1_type>) {
static_assert(
has_method_doBoundaryCrossing_v<TSequence, ProcessReturn, TParticle&>,
@@ -56,13 +54,11 @@ namespace corsika {
}
} else {
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
- process2_type> ||
+ if constexpr (is_boundary_process_v<process2_type> ||
process2_type::is_process_sequence) {
// interface checking on USequence
- if constexpr (std::is_base_of_v<BoundaryCrossingProcess<process2_type>,
- process2_type>) {
+ if constexpr (is_boundary_process_v<process2_type>) {
static_assert(
has_method_doBoundaryCrossing_v<USequence, ProcessReturn, TParticle>,
@@ -136,7 +132,7 @@ namespace corsika {
IndexProcess2>::doSecondaries(TSecondaries& vS) {
const auto& particle = vS.parent();
if (select_(particle)) {
- if constexpr (std::is_base_of_v<SecondariesProcess<process1_type>, process1_type> ||
+ if constexpr (is_secondaries_process_v<process1_type> ||
process1_type::is_process_sequence) {
// interface checking on TSequence
@@ -150,7 +146,7 @@ namespace corsika {
A_.doSecondaries(vS);
}
} else {
- if constexpr (std::is_base_of_v<SecondariesProcess<process2_type>, process2_type> ||
+ if constexpr (is_secondaries_process_v<process2_type> ||
process2_type::is_process_sequence) {
// interface checking on USequence
@@ -219,14 +215,14 @@ namespace corsika {
IndexProcess2>::getInverseInteractionLength(TParticle&& particle) {
if (select_(particle)) {
- if constexpr (std::is_base_of_v<InteractionProcess<process1_type>, process1_type> ||
+ if constexpr (is_interaction_process_v<process1_type> ||
process1_type::is_process_sequence) {
return A_.getInverseInteractionLength(particle);
}
} else {
- if constexpr (std::is_base_of_v<InteractionProcess<process2_type>, process2_type> ||
+ if constexpr (is_interaction_process_v<process2_type> ||
process2_type::is_process_sequence) {
return B_.getInverseInteractionLength(particle);
}
@@ -249,8 +245,7 @@ namespace corsika {
A_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
// if A_ did succeed, stop routine. Not checking other static branch B_.
if (ret != ProcessReturn::Ok) { return ret; }
- } else if constexpr (std::is_base_of_v<InteractionProcess<process1_type>,
- process1_type>) {
+ } else if constexpr (is_interaction_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += A_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -272,8 +267,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
return B_.selectInteraction(view, lambda_inv_select, lambda_inv_sum);
- } else if constexpr (std::is_base_of_v<InteractionProcess<process2_type>,
- process2_type>) {
+ } else if constexpr (is_interaction_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
lambda_inv_sum += B_.getInverseInteractionLength(view.parent());
// check if we should execute THIS process and then EXIT
@@ -301,14 +295,14 @@ namespace corsika {
IndexProcess2>::getInverseLifetime(TParticle&& particle) {
if (select_(particle)) {
- if constexpr (std::is_base_of_v<DecayProcess<process1_type>, process1_type> ||
+ if constexpr (is_decay_process_v<process1_type> ||
process1_type::is_process_sequence) {
return A_.getInverseLifetime(particle);
}
} else {
- if constexpr (std::is_base_of_v<DecayProcess<process2_type>, process2_type> ||
+ if constexpr (is_decay_process_v<process2_type> ||
process2_type::is_process_sequence) {
return B_.getInverseLifetime(particle);
}
@@ -331,8 +325,7 @@ namespace corsika {
ProcessReturn const ret = A_.selectDecay(view, decay_inv_select, decay_inv_sum);
// if A_ did succeed, stop routine here (not checking other static branch B_)
if (ret != ProcessReturn::Ok) { return ret; }
- } else if constexpr (std::is_base_of_v<DecayProcess<process1_type>,
- process1_type>) {
+ } else if constexpr (is_decay_process_v<process1_type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += A_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
@@ -355,8 +348,7 @@ namespace corsika {
if constexpr (process2_type::is_process_sequence) {
// if B_ is a process sequence --> check inside
return B_.selectDecay(view, decay_inv_select, decay_inv_sum);
- } else if constexpr (std::is_base_of_v<DecayProcess<process2_type>,
- process2_type>) {
+ } else if constexpr (is_decay_process_v<process2_type>) {
// if this is not a ContinuousProcess --> evaluate probability
decay_inv_sum += B_.getInverseLifetime(view.parent());
// check if we should execute THIS process and then EXIT
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
new file mode 100644
index 000000000..efeb48c34
--- /dev/null
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -0,0 +1,299 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <cmath>
+
+namespace corsika {
+
+ namespace andre {
+
+ //---------------------------------------------------------------------------
+ // solve cubic equation x^3 + a*x^2 + b*x + c = 0
+ // x - array of size 3
+ // In case 3 real roots: => x[0], x[1], x[2], return 3
+ // 2 real roots: x[0], x[1], return 2
+ // 1 real root : x[0], x[1] ± i*x[2], return 1
+ inline std::vector<double> solveP3(long double A, long double B, long double C,
+ long double D, double const epsilon) {
+
+ if (std::abs(A) < epsilon) { return solve_quadratic_real(B, C, epsilon); }
+
+ long double a = B / A;
+ long double b = C / A;
+ long double c = D / A;
+
+ long double a2 = a * a;
+ long double q = (a2 - 3 * b) / 9;
+ long double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
+ long double r2 = r * r;
+ long double q3 = q * q * q;
+
+ // disc = q**3 + r**2
+ // w:e = r*r exactly
+ long double w = r * r;
+ long double e = std::fma(r, r, -w);
+ // s:t = q*q exactly
+ long double s = -q * q;
+ long double t = std::fma(-q, q, -s);
+ // s:t * q + w:e = s*q + w + t*q +e = s*q+w + u:v + e = f + u:v + e
+ long double f = std::fma(s, q, w);
+ long double u = t * q;
+ long double v = std::fma(t, q, -u);
+ // error-free sum f+u. See Knuth, TAOCP vol. 2
+ long double uf = u + f;
+ long double au = uf - u;
+ long double ab = uf - au;
+ au = f - au;
+ ab = u - ab;
+ // sum all terms into final result
+ long double const disc = -(((e + uf) + au) + ab) + v;
+
+ CORSIKA_LOG_TRACE("solveP3 disc={} r2-r3={}", disc, r2 - q3);
+
+ if (disc >= 0) {
+ long double t = r / std::sqrt(q3);
+ if (t < -1) t = -1;
+ if (t > 1) t = 1;
+ t = std::acos(t);
+ a /= 3;
+ q = -2 * std::sqrt(q);
+ return {double(q * std::cos(t / 3) - a),
+ double(q * std::cos((t + 2 * M_PI) / 3) - a),
+ double(q * std::cos((t - 2 * M_PI) / 3) - a)};
+ } else {
+ long double A = -std::pow(std::fabs(r) + std::sqrt(-disc), 1. / 3);
+ if (r < 0) A = -A;
+ long double B = (0 == A ? 0 : q / A);
+
+ a /= 3;
+ long double test = 0.5 * std::sqrt(3.) * (A - B);
+ if (std::fabs(test) < epsilon) {
+ return {double((A + B) - 1), double(-0.5 * (A + B) - a)};
+ }
+
+ return {double((A + B) - a)};
+ }
+ }
+ } // namespace andre
+
+ inline std::vector<double> solve_cubic_depressed_disciminant_real(
+ long double p, long double q, long double const disc, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("p={}, q={}, disc={}", p, q, disc);
+
+ if (std::abs(disc) < epsilon) { // disc==0 multiple roots !
+ if (std::abs(p) < epsilon) { // tripple root
+ return {0};
+ }
+ // double root, single root
+ CORSIKA_LOG_TRACE("cubic double root");
+ return {double(-3 * q / (2 * p)), double(3 * q / p)};
+ }
+
+ if (std::abs(p) < epsilon) { // p==0 --> x^3 + q = 0
+ return {double(std::cbrt(-q))};
+ }
+
+ if (disc > 0) { // three real roots
+ CORSIKA_LOG_TRACE("cubic three roots");
+ long double const p_third = p / 3;
+ long double const sqrt_minus_p_third = std::sqrt(-p_third);
+
+ long double const arg = std::acos(q / (2 * p_third) / sqrt_minus_p_third) / 3;
+
+ long double constexpr pi = M_PI;
+ return {double(2 * sqrt_minus_p_third * std::cos(arg)),
+ double(2 * sqrt_minus_p_third * std::cos(arg - 2 * pi / 3)),
+ double(2 * sqrt_minus_p_third * std::cos(arg - 4 * pi / 3))};
+ }
+
+ // thus disc < 0 -> one real root
+ if (p < 0) {
+ CORSIKA_LOG_TRACE("cubic cosh");
+ long double const abs_q = std::abs(q);
+ long double const p_third = p / 3;
+ long double const sqrt_minus_p_third = std::sqrt(-p_third);
+ CORSIKA_LOG_TRACE("sqrt_minus_p_third={}, arcosh({})={}", sqrt_minus_p_third,
+ -abs_q / (2 * p_third) / sqrt_minus_p_third,
+ std::acosh(-abs_q / (2 * p_third) / sqrt_minus_p_third));
+ CORSIKA_LOG_TRACE(
+ "complex: {}",
+ -2 * abs_q / q * sqrt_minus_p_third *
+ std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3));
+ double const z =
+ -2 * abs_q / q * sqrt_minus_p_third *
+ std::cosh(std::acosh(-abs_q / (2 * p_third * sqrt_minus_p_third)) / 3);
+ return {z};
+ } else { // p > 0
+ CORSIKA_LOG_TRACE("cubic sinh");
+ long double const p_third = p / 3;
+ long double const sqrt_p_third = std::sqrt(p_third);
+ return {double(-2 * sqrt_p_third *
+ std::sinh(std::asinh(q / (2 * p_third * sqrt_p_third)) / 3))};
+ }
+ }
+
+ inline std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+ double const epsilon) {
+
+ // thanks!
+ // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+ // long double const disc = -(4 * p * p * p + 27 * q * q);
+ // disc = (p/3)**3 + (q/2)**2
+ long double p_third = p / 3;
+ long double q_half = q / 2;
+ // w:e = (q/2)*(q/2) exactly
+ long double w = q_half * q_half;
+ long double e = std::fma(q_half, q_half, -w);
+ // s:t = (p/3)*(p/3) exactly
+ long double s = p_third * p_third;
+ long double t = std::fma(p_third, p_third, -s);
+ // s:t * p + w:e = s*p + w + t*p +e = s*p+w + u:v + e = f + u:v + e
+ long double f = std::fma(s, p_third, w);
+ long double u = t * p_third;
+ long double v = std::fma(t, p_third, -u);
+ // error-free sum f+u. See Knuth, TAOCP vol. 2
+ long double a = u + f;
+ long double b = a - u;
+ long double c = a - b;
+ b = f - b;
+ c = u - c;
+ // sum all terms into final result
+ long double const disc = -(((e + a) + b) + c) + v;
+ return solve_cubic_depressed_disciminant_real(p, q, disc, epsilon);
+ }
+
+ /**
+ Analytical approach. Not very stable in all conditions.
+ */
+ inline std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+ long double c, long double d,
+ double const epsilon) {
+
+ CORSIKA_LOG_TRACE("cubic: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b, c,
+ d, epsilon, (std::abs(a - 1) < epsilon), (std::abs(b) < epsilon));
+
+ if (std::abs(a) < epsilon) { // this is just a quadratic
+ return solve_quadratic_real(b, c, d, epsilon);
+ }
+
+ if ((std::abs(a - 1) < epsilon) &&
+ (std::abs(b) < epsilon)) { // this is a depressed cubic
+ return solve_cubic_depressed_real(c, d, epsilon);
+ }
+
+ // p = (3ac - b^2) / (3a^2) = 3 * ( c/(3*a) - b^2/(9*a^2) )
+ long double b_over_a = b / a;
+ long double const p_third = std::fma(-b_over_a, b_over_a / 9, c / (a * 3));
+
+ // p = std::fma(a * 3, c, -b2) / (3 * a2);
+ // q = (2*b^3 - 9*abc + 27*a^2*d) / (27a^3) = 2 * ( b^3/(27*a^3) - bc/(6*a^2) +
+ // d/(2*a) )
+ long double const q_half_term1 = std::fma(b_over_a, b_over_a / 27, -c / (a * 6));
+ long double const q_half = std::fma(b_over_a, q_half_term1, d / (a * 2));
+
+ std::vector<double> zs = solve_cubic_depressed_real(3 * p_third, 2 * q_half, epsilon);
+ CORSIKA_LOG_TRACE("cubic: solve_depressed={}, b/3a={}", fmt::join(zs, ", "),
+ b / (a * 3));
+ for (auto& z : zs) {
+ z -= b / (a * 3);
+ double const b1 = z + b / a;
+ double const b0 = b1 * z + c / a;
+ std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+ CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+ std::pow(z, 3) * a + std::pow(z, 2) * b + std::pow(z, 1) * c + d);
+ }
+ CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(zs, ", "));
+ return zs;
+ }
+
+ inline long double cubic_function(long double x, long double a, long double b,
+ long double c, long double d) {
+ return std::pow(x, 3) * a + std::pow(x, 2) * b + x * c + d;
+ }
+
+ inline long double cubic_function_dfdx(long double x, long double a, long double b,
+ long double c) {
+ return std::pow(x, 2) * a * 3 + x * b * 2 + c;
+ }
+
+ inline long double cubic_function_d2fd2x(long double x, long double a, long double b) {
+ return x * a * 6 + b * 2;
+ }
+
+ /**
+ Iterative approach.
+ */
+
+ inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
+ long double d, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("cubic_2: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b,
+ c, d, epsilon, (std::abs(a - 1) < epsilon),
+ (std::abs(b) < epsilon));
+
+ if (std::abs(a) < epsilon) { // this is just a quadratic
+ return solve_quadratic_real(b, c, d, epsilon);
+ }
+
+ long double const dist = std::fma(b / a, b / a, -3 * c / a);
+ long double const xinfl = -b / (a * 3);
+
+ long double x1 = xinfl;
+ long double f_x1 = cubic_function(xinfl, a, b, c, d);
+
+ CORSIKA_LOG_TRACE("dist={} xinfl={} f_x1={} {} {}", dist, xinfl, f_x1,
+ cubic_function(xinfl - 2 / 3 * std::sqrt(dist), a, b, c, d),
+ cubic_function(xinfl + 2 / 3 * std::sqrt(dist), a, b, c, d));
+
+ if (std::abs(f_x1) > epsilon) {
+ if (std::abs(dist) < epsilon) {
+ x1 = xinfl - std::cbrt(f_x1);
+ } else if (dist > 0) {
+ if (f_x1 > 0)
+ x1 = xinfl - 2 / 3 * std::sqrt(dist);
+ else
+ x1 = xinfl + 2 / 3 * std::sqrt(dist);
+ }
+
+ int niter = 0;
+ const int maxiter = 100;
+ do {
+ long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
+ long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
+ // if (potential) saddle point... avoid
+ if (std::abs(f_prime_x1) < epsilon) {
+ x1 -= std::cbrt(f_x1);
+ } else {
+ x1 -= f_x1 * f_prime_x1 / (std::pow(f_prime_x1, 2) - 0.5 * f_x1 * f_prime2_x1);
+ }
+ f_x1 = cubic_function(x1, a, b, c, d);
+ CORSIKA_LOG_TRACE("niter={} x1={} f_x1={} f_prime={} f_prime2={} eps={}", niter,
+ x1, f_x1, f_prime_x1, f_prime2_x1, epsilon);
+ } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon));
+
+ CORSIKA_LOG_TRACE("niter={}", niter);
+ }
+
+ CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1);
+
+ double const b1 = x1 + b / a;
+ double const b0 = b1 * x1 + c / a;
+ std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
+ CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
+ cubic_function(x1, a, b, c, d));
+
+ quad_check.push_back(x1);
+ CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(quad_check, ", "));
+ return quad_check;
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/LinearSolver.inl b/corsika/detail/framework/utility/LinearSolver.inl
new file mode 100644
index 000000000..45ed7a517
--- /dev/null
+++ b/corsika/detail/framework/utility/LinearSolver.inl
@@ -0,0 +1,36 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+ inline std::vector<double> solve_linear_real(double a, double b, double const epsilon) {
+
+ if (std::abs(a) < epsilon) {
+ return {}; // no (b!=0), or infinite number (b==0) of solutions....
+ }
+
+ return {-b / a};
+ }
+
+ inline std::vector<std::complex<double>> solve_linear(double a, double b,
+ double const epsilon) {
+
+ if (std::abs(a) < epsilon) {
+ return {}; // no (b!=0), or infinite number (b==0) of solutions....
+ }
+
+ return {{-b / a, 0}};
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl
new file mode 100644
index 000000000..d3b3f8671
--- /dev/null
+++ b/corsika/detail/framework/utility/QuadraticSolver.inl
@@ -0,0 +1,81 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+namespace corsika {
+
+ inline std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+ long double c,
+ double const epsilon) {
+ if (std::abs(a) < epsilon) { return solve_linear(b, c, epsilon); }
+
+ if (std::abs(c) < epsilon) {
+ std::vector<std::complex<double>> lin_result = solve_linear(a, b, epsilon);
+ lin_result.push_back({0.});
+ return lin_result;
+ }
+
+ long double const radicant = std::pow(b, 2) - a * c * 4;
+
+ if (radicant < -epsilon) { // two complex solutions
+ double const rpart = -b / 2 * a;
+ double const ipart = std::sqrt(-radicant);
+ return {{rpart, ipart}, {rpart, -ipart}};
+ }
+
+ if (radicant < epsilon) { // just one real solution
+ return {{double(-b / 2 * a), 0}};
+ }
+
+ // two real solutions, use Citardauq formula and actively avoid cancellation in the
+ // denominator
+
+ const long double x1 =
+ 2 * c / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+ return {{double(x1), 0}, {double(c / (a * x1)), 0}};
+ }
+
+ inline std::vector<double> solve_quadratic_real(long double a, long double b,
+ long double c, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("quadratic: a={} b={} c={}", a, b, c);
+
+ if (std::abs(a) < epsilon) { return solve_linear_real(b, c, epsilon); }
+ if (std::abs(c) < epsilon) {
+ std::vector<double> lin_result = solve_linear_real(a, b, epsilon);
+ lin_result.push_back(0.);
+ return lin_result;
+ }
+
+ // long double const radicant = std::pow(b, 2) - a * c * 4;
+ // Thanks!
+ // https://math.stackexchange.com/questions/2434354/numerically-stable-scheme-for-the-3-real-roots-of-a-cubic
+ long double w = a * 4 * c;
+ long double e = std::fma(a * 4, c, -w);
+ long double f = std::fma(b, b, -w);
+ long double radicant = f + e;
+
+ CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4);
+
+ if (std::abs(radicant) < epsilon) { // just one real solution
+ return {double(-b / (2 * a))};
+ }
+
+ if (radicant < 0) { return {}; }
+
+ // two real solutions, use Citardauq formula and actively avoid cancellation in the
+ // denominator
+
+ long double const x1 =
+ c * 2 / (b > 0 ? -b - std::sqrt(radicant) : -b + std::sqrt(radicant));
+
+ CORSIKA_LOG_TRACE("quadratic x1={} x2={}", double(x1), double(c / (a * x1)));
+
+ return {double(x1), double(c / (a * x1))};
+ }
+} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index a82abf7ac..2bed4396c 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -8,128 +8,179 @@
#pragma once
+#include <corsika/framework/utility/CubicSolver.hpp>
#include <cmath>
-namespace corsika::quartic_solver {
-
- //---------------------------------------------------------------------------
- // solve cubic equation x^3 + a*x^2 + b*x + c = 0
- // x - array of size 3
- // In case 3 real roots: => x[0], x[1], x[2], return 3
- // 2 real roots: x[0], x[1], return 2
- // 1 real root : x[0], x[1] ± i*x[2], return 1
- inline unsigned int solveP3(double* x, double a, double b, double c) {
- double a2 = a * a;
- double q = (a2 - 3 * b) / 9;
- double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
- double r2 = r * r;
- double q3 = q * q * q;
- double A, B;
- if (r2 < q3) {
- double t = r / sqrt(q3);
- if (t < -1) t = -1;
- if (t > 1) t = 1;
- t = acos(t);
- a /= 3;
- q = -2 * sqrt(q);
- x[0] = q * cos(t / 3) - a;
- x[1] = q * cos((t + M_2PI) / 3) - a;
- x[2] = q * cos((t - M_2PI) / 3) - a;
- return 3;
- } else {
- A = -pow(fabs(r) + sqrt(r2 - q3), 1. / 3);
- if (r < 0) A = -A;
- B = (0 == A ? 0 : q / A);
-
- a /= 3;
- x[0] = (A + B) - a;
- x[1] = -0.5 * (A + B) - a;
- x[2] = 0.5 * sqrt(3.) * (A - B);
- if (fabs(x[2]) < eps) {
- x[2] = x[1];
- return 2;
+namespace corsika {
+
+ namespace andre {
+
+ inline std::vector<double> solve_quartic_real(long double a, long double b,
+ long double c, long double d,
+ long double e, double const epsilon) {
+
+ if (std::abs(a) < epsilon) { return solve_cubic_real(b, c, d, e, epsilon); }
+
+ b /= a;
+ c /= a;
+ d /= a;
+ e /= a;
+
+ long double a3 = -c;
+ long double b3 = b * d - 4. * e;
+ long double c3 = -b * b * e - d * d + 4. * c * e;
+
+ // cubic resolvent
+ // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+
+ std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
+ long double y = x3[0]; // there is always at least one solution
+ // The essence - choosing Y with maximal absolute value.
+ if (x3.size() == 3) {
+ if (fabs(x3[1]) > fabs(y)) y = x3[1];
+ if (fabs(x3[2]) > fabs(y)) y = x3[2];
+ }
+
+ long double q1, q2, p1, p2;
+ // h1+h2 = y && h1*h2 = e <=> h^2 -y*h + e = 0 (h === q)
+
+ long double Det = y * y - 4 * e;
+ CORSIKA_LOG_TRACE("Det={}", Det);
+ if (fabs(Det) < epsilon) // in other words - D==0
+ {
+ q1 = q2 = y * 0.5;
+ // g1+g2 = b && g1+g2 = c-y <=> g^2 - b*g + c-y = 0 (p === g)
+ Det = b * b - 4 * (c - y);
+ if (fabs(Det) < epsilon) { // in other words - D==0
+ p1 = p2 = b * 0.5;
+ } else {
+ if (Det < 0) return {};
+ long double sqDet = sqrt(Det);
+ p1 = (b + sqDet) * 0.5;
+ p2 = (b - sqDet) * 0.5;
+ }
+ } else {
+ if (Det < 0) return {};
+ long double sqDet = sqrt(Det);
+ q1 = (y + sqDet) * 0.5;
+ q2 = (y - sqDet) * 0.5;
+ // g1+g2 = b && g1*h2 + g2*h1 = c ( && g === p ) Krammer
+ p1 = (b * q1 - d) / (q1 - q2);
+ p2 = (d - b * q2) / (q1 - q2);
}
- return 1;
+ // solving quadratic eqs. x^2 + p1*x + q1 = 0
+ // x^2 + p2*x + q2 = 0
+ {
+ long double Det1 = p1 * p1 - 4 * q1;
+ long double Det2 = p2 * p2 - 4 * q2;
+ CORSIKA_LOG_TRACE("Det1={} Det2={}", Det1, Det2);
+ if (Det1 >= 0 && Det2 >= 0) {
+ long double sqDet1 = sqrt(Det1);
+ long double sqDet2 = sqrt(Det2);
+ // return {double(-p1 + sqDet1) * 0.5, double(-p1 - sqDet1) * 0.5, double(-p2 +
+ // sqDet2) * 0.5, double(-p2 - sqDet2) * 0.5};
+ CORSIKA_LOG_TRACE("check1 {} {} {} {}", double(-p1 + sqDet1) * 0.5,
+ double(-p1 - sqDet1) * 0.5, double(-p2 + sqDet2) * 0.5,
+ double(-p2 - sqDet2) * 0.5);
+ }
+ if (Det1 >= 0) {
+ long double sqDet1 = sqrt(Det1);
+ // return {double(-p1 + sqDet1) * 0.5, double(-p1 - sqDet1) * 0.5};
+ CORSIKA_LOG_TRACE("check2 {} {} ", double(-p1 + sqDet1) * 0.5,
+ double(-p1 - sqDet1) * 0.5);
+ }
+ if (Det2 >= 0) {
+ long double sqDet2 = sqrt(Det2);
+ // return {double(-p2 + sqDet2) * 0.5, double(-p2 - sqDet2) * 0.5};
+ CORSIKA_LOG_TRACE("check3 {} {}", double(-p1 + sqDet2) * 0.5,
+ double(-p1 - sqDet2) * 0.5);
+ }
+ }
+
+ std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
+ std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
+ if (quad2.size() > 0) {
+ for (auto val : quad2) quad1.push_back(val);
+ }
+ return quad1;
}
- }
+ } // namespace andre
- //---------------------------------------------------------------------------
- // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
- // Attention - this function returns dynamically allocated array. It has to be released
- // afterwards.
- inline DComplex* solve_quartic(double a, double b, double c, double d) {
- double a3 = -b;
- double b3 = a * c - 4. * d;
- double c3 = -a * a * d - c * c + 4. * b * d;
+ inline std::vector<double> solve_quartic_depressed_real(long double p, long double q,
+ long double r,
+ double const epsilon) {
- // cubic resolvent
- // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+ CORSIKA_LOG_TRACE("quartic-depressed: p={:f}, q={:f}, r={:f}, epsilon={}", p, q, r,
+ epsilon);
- double x3[3];
- unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+ long double const p2 = std::pow(p, 2);
+ long double const q2 = std::pow(q, 2);
- double q1, q2, p1, p2, D, sqD, y;
+ std::vector<double> const resolve_cubic =
+ solve_cubic_real(1, p, p2 / 4 - r, -q2 / 8, epsilon);
- y = x3[0];
- // The essence - choosing Y with maximal absolute value.
- if (iZeroes != 1) {
- if (fabs(x3[1]) > fabs(y)) y = x3[1];
- if (fabs(x3[2]) > fabs(y)) y = x3[2];
- }
+ CORSIKA_LOG_TRACE("resolve_cubic: N={}, m=[{}]", resolve_cubic.size(),
+ fmt::join(resolve_cubic, ", "));
- // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
-
- D = y * y - 4 * d;
- if (fabs(D) < eps) // in other words - D==0
- {
- q1 = q2 = y * 0.5;
- // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
- D = a * a - 4 * (b - y);
- if (fabs(D) < eps) // in other words - D==0
- p1 = p2 = a * 0.5;
-
- else {
- sqD = sqrt(D);
- p1 = (a + sqD) * 0.5;
- p2 = (a - sqD) * 0.5;
- }
- } else {
- sqD = sqrt(D);
- q1 = (y + sqD) * 0.5;
- q2 = (y - sqD) * 0.5;
- // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
- p1 = (a * q1 - c) / (q1 - q2);
- p2 = (c - a * q2) / (q1 - q2);
+ if (!resolve_cubic.size()) return {};
+
+ long double m = 0;
+ for (auto const& v : resolve_cubic) {
+ CORSIKA_LOG_TRACE("check pol3(v)={}", (std::pow(v, 3) + std::pow(v, 2) * p +
+ v * (p2 / 4 - r) - q2 / 8));
+ if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
}
+ CORSIKA_LOG_TRACE("check m={}", m);
+ if (m == 0) { return {{0}}; }
+
+ CORSIKA_LOG_TRACE("check m={}", m);
+
+ long double const quad_term1 = p / 2 + m;
+ long double const quad_term2 = std::sqrt(2 * m);
+ long double const quad_term3 = q / (2 * quad_term2);
- DComplex* retval = new DComplex[4];
-
- // solving quadratic eq. - x^2 + p1*x + q1 = 0
- D = p1 * p1 - 4 * q1;
- if (D < 0.0) {
- retval[0].real(-p1 * 0.5);
- retval[0].imag(sqrt(-D) * 0.5);
- retval[1] = std::conj(retval[0]);
- } else {
- sqD = sqrt(D);
- retval[0].real((-p1 + sqD) * 0.5);
- retval[1].real((-p1 - sqD) * 0.5);
+ std::vector<double> z_quad1 =
+ solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, epsilon);
+ std::vector<double> z_quad2 =
+ solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, epsilon);
+ for (auto const& z : z_quad2) z_quad1.push_back(z);
+ return z_quad1;
+ }
+
+ inline std::vector<double> solve_quartic_real(long double a, long double b,
+ long double c, long double d,
+ long double e, double const epsilon) {
+
+ CORSIKA_LOG_TRACE("quartic: a={:f}, b={:f}, c={:f}, d={:f}, e={:f}, epsilon={}", a, b,
+ c, d, e, epsilon);
+
+ if (std::abs(a) < epsilon) { // this is just a quadratic
+ return solve_cubic_real(b, c, d, e, epsilon);
}
- // solving quadratic eq. - x^2 + p2*x + q2 = 0
- D = p2 * p2 - 4 * q2;
- if (D < 0.0) {
- retval[2].real(-p2 * 0.5);
- retval[2].imag(sqrt(-D) * 0.5);
- retval[3] = std::conj(retval[2]);
- } else {
- sqD = sqrt(D);
- retval[2].real((-p2 + sqD) * 0.5);
- retval[3].real((-p2 - sqD) * 0.5);
+ if ((std::abs(a - 1) < epsilon) &&
+ (std::abs(b) < epsilon)) { // this is a depressed quartic
+ return solve_quartic_depressed_real(c, d, e, epsilon);
}
- return retval;
+ long double const b2 = std::pow(b, 2);
+ long double const b3 = std::pow(b, 3);
+ long double const b4 = std::pow(b, 4);
+ long double const a2 = std::pow(a, 2);
+ long double const a3 = std::pow(a, 3);
+ long double const a4 = std::pow(a, 4);
+
+ long double const p = (c * a * 8 - b2 * 3) / (a4 * 8);
+ long double const q = (b3 - b * c * a * 4 + d * a2 * 8) / (a4 * 8);
+ long double const r =
+ (-b4 * 3 + e * a3 * 256 - b * d * a2 * 64 + b2 * c * a * 16) / (a4 * 256);
+
+ std::vector<double> zs = solve_quartic_depressed_real(p, q, r, epsilon);
+ CORSIKA_LOG_TRACE("quartic: solve_depressed={}, b/4a={}", fmt::join(zs, ", "),
+ b / (4 * a));
+ for (auto& z : zs) { z -= b / (4 * a); }
+ CORSIKA_LOG_TRACE("quartic: solve_quartic_real returns={}", fmt::join(zs, ", "));
+ return zs;
}
-
-} // namespace corsika::quartic_solver
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/detail/modules/ObservationPlane.inl b/corsika/detail/modules/ObservationPlane.inl
index c4da50e94..e1fe6a504 100644
--- a/corsika/detail/modules/ObservationPlane.inl
+++ b/corsika/detail/modules/ObservationPlane.inl
@@ -45,6 +45,8 @@ namespace corsika {
<< ' ' << displacement.dot(xAxis_) / 1_m << ' '
<< displacement.dot(yAxis_) / 1_m << '\n';
+ CORSIKA_LOG_TRACE("Particle detected absorbed={}", deleteOnHit_);
+
if (deleteOnHit_) {
count_ground_++;
energy_ground_ += energy;
@@ -58,13 +60,14 @@ namespace corsika {
corsika::setup::Stack::particle_type const& particle,
corsika::setup::Trajectory const& trajectory) {
+ CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(),
+ particle.getPosition(), particle.getDirection(), plane_.asString());
+
Intersections const intersection =
setup::Tracking::intersect<corsika::setup::Stack::particle_type>(particle,
plane_);
TimeType const timeOfIntersection = intersection.getEntry();
- CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}, timeOfIntersection={}",
- particle.asString(), particle.getPosition(),
- particle.getDirection(), plane_.asString(), timeOfIntersection);
+ CORSIKA_LOG_TRACE("timeOfIntersection={}", timeOfIntersection);
if (timeOfIntersection < TimeType::zero()) {
return std::numeric_limits<double>::infinity() * 1_m;
}
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 5770fe439..f780175b5 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -34,15 +34,14 @@ namespace corsika {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
double(0));
} // charge of the particle
- int const chargeNumber = particle.getChargeNumber();
+ ElectricChargeType const charge = particle.getCharge();
auto const* currentLogicalVolumeNode = particle.getNode();
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
- VelocityVector velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
+ VelocityVector velocity = particle.getVelocity();
decltype(meter / (second * volt)) k =
- chargeNumber * constants::cSquared * 1_eV /
+ charge * constants::cSquared * 1_eV /
(velocity.getNorm() * particle.getEnergy() * 1_V);
DirectionVector direction = velocity.normalized();
auto position = particle.getPosition(); // First Movement
@@ -58,8 +57,7 @@ namespace corsika {
template <typename TParticle>
inline auto Tracking::getTrack(TParticle const& particle) {
- VelocityVector const initialVelocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
+ VelocityVector const initialVelocity = particle.getVelocity();
auto const position = particle.getPosition();
CORSIKA_LOG_DEBUG(
@@ -74,7 +72,7 @@ namespace corsika {
typedef
typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
- node_type& volumeNode = *particle.getNode();
+ node_type const& volumeNode = *particle.getNode();
// for the event of magnetic fields and curved trajectories, we need to limit
// maximum step-length since we need to follow curved
@@ -83,25 +81,25 @@ namespace corsika {
MagneticFieldVector const& magneticfield =
volumeNode.getModelProperties().getMagneticField(position);
MagneticFluxType const magnitudeB = magneticfield.getNorm();
- int const chargeNumber = particle.getChargeNumber();
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+ ElectricChargeType const charge = particle.getCharge();
+ bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
if (no_deflection) { return getLinearTrajectory(particle); }
- HEPMomentumType const pAlongB_delta =
+ HEPMomentumType const p_perp =
(particle.getMomentum() -
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
- if (pAlongB_delta == 0_GeV) {
+ if (p_perp < 1_eV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
- CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+ CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return getLinearTrajectory(particle);
}
- LengthType const gyroradius =
- (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+ constants::c / (abs(charge) * magnitudeB));
double const maxRadians = 0.01;
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
@@ -113,8 +111,15 @@ namespace corsika {
// intersection
auto [minTime, minNode] = nextIntersect(particle, steplimit_time);
- auto const k =
- chargeNumber * constants::cSquared * 1_eV / (particle.getEnergy() * 1_V);
+ auto const p_norm =
+ constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()); // kg *m /s
+ // k = q/|p|
+ decltype(1 / (tesla * second)) const k =
+ charge / p_norm *
+ initialVelocity.getNorm(); // since we use steps in time and not length
+ // units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1/(T s)
+
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
minTime), // trajectory
@@ -129,99 +134,104 @@ namespace corsika {
return Intersections();
}
- int const chargeNumber = particle.getChargeNumber();
+ ElectricChargeType const charge = particle.getCharge();
auto const& position = particle.getPosition();
auto const* currentLogicalVolumeNode = particle.getNode();
MagneticFieldVector const& magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
- VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- DirectionVector const directionBefore =
- velocity.normalized(); // determine steplength to next volume
+ VelocityVector const velocity = particle.getVelocity();
+ DirectionVector const directionBefore = velocity.normalized();
auto const projectedDirection = directionBefore.cross(magneticfield);
auto const projectedDirectionSqrNorm = projectedDirection.getSquaredNorm();
bool const isParallel = (projectedDirectionSqrNorm == 0 * square(1_T));
- if (chargeNumber == 0 || magneticfield.getNorm() == 0_T || isParallel) {
+ if ((charge == 0 * constants::e) || magneticfield.getNorm() == 0_T || isParallel) {
return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
}
bool const numericallyInside = sphere.contains(particle.getPosition());
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
- auto const absVelocity = velocity.getNorm();
- auto const energy = particle.getEnergy();
+ SpeedType const absVelocity = velocity.getNorm();
+ auto const p_norm =
+ constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()); // km * m /s
// this is: k = q/|p|
- auto const k =
- chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
+ auto const k = charge / p_norm;
- auto const direction_B_perp = directionBefore.cross(magneticfield);
- auto const denom = 4. / (direction_B_perp.getSquaredNorm() * k * k);
+ MagneticFieldVector const direction_x_B = directionBefore.cross(magneticfield);
+ auto const denom = 4. / (direction_x_B.getSquaredNorm() * k * k);
Vector<length_d> const deltaPos = position - sphere.getCenter();
- double const b = (direction_B_perp.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
+ double const b = (direction_x_B.dot(deltaPos) * k + 1) * denom / (1_m * 1_m);
double const c = directionBefore.dot(deltaPos) * 2 * denom / (1_m * 1_m * 1_m);
LengthType const deltaPosLength = deltaPos.getNorm();
double const d = (deltaPosLength + sphere.getRadius()) *
(deltaPosLength - sphere.getRadius()) * denom /
(1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
- std::complex<double> const* solutions = quartic_solver::solve_quartic(0, b, c, d);
+ std::vector<double> solutions = andre::solve_quartic_real(1, 0, b, c, d);
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
- for (int i = 0; i < 4; i++) {
- if (solutions[i].imag() == 0) {
- LengthType const dist = solutions[i].real() * 1_m;
- CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
- if (numericallyInside) {
- // there must be an entry (negative) and exit (positive) solution
- if (dist < -0.0001_m) { // security margin to assure transfer to next
- // logical volume
- if (first_entry == 0) {
- d_enter = dist;
- } else {
- d_enter = std::max(d_enter, dist); // closest negative to zero (-1e-4) m
- }
- first_entry++;
-
- } else { // thus, dist > -0.0001_m
-
- if (first_exit == 0) {
- d_exit = dist;
- } else {
- d_exit = std::min(d_exit, dist); // closest positive to zero (-1e-4) m
- }
- first_exit++;
+ for (auto solution : solutions) {
+ LengthType const dist = solution * 1_m;
+ CORSIKA_LOG_TRACE("Solution (real) for current Volume: {} ", dist);
+ if (numericallyInside) {
+ // there must be an entry (negative) and exit (positive) solution
+ if (dist < 0.0001_m) { // security margin to assure
+ // transfer to next logical volume
+ // (even if dist suggest marginal
+ // entry already, which we
+ // classify as numerical artifact)
+ if (first_entry == 0) {
+ d_enter = dist;
+ } else {
+ d_enter = std::max(d_enter, dist); // closest negative to zero >1e-4 m
}
- first = int(first_exit > 0) + int(first_entry > 0);
+ first_entry++;
- } else { // thus, numericallyInside == false
+ } else { // thus, dist > +0.0001_m
- // both physical solutions (entry, exit) must be positive, and as small as
- // possible
- if (dist < -0.0001_m) { // need small numerical margin, to assure transport
- // into next logical volume
- continue;
+ if (first_exit == 0) {
+ d_exit = dist;
+ } else {
+ d_exit = std::min(d_exit, dist); // closest positive to zero >1e-4 m
}
- if (first == 0) {
+ first_exit++;
+ }
+ first = int(first_exit > 0) + int(first_entry > 0);
+
+ } else { // thus, numericallyInside == false
+
+ // both physical solutions (entry, exit) must be positive, and as small as
+ // possible
+ if (dist < -0.0001_m) { // need small numerical margin, to
+ // assure transport. We consider
+ // begin marginally already inside
+ // next volume (besides
+ // numericallyInside=false) as numerical glitch.
+ // into next logical volume
+ continue;
+ }
+ if (first == 0) {
+ d_enter = dist;
+ } else {
+ if (dist < d_enter) {
+ d_exit = d_enter;
d_enter = dist;
} else {
- if (dist < d_enter) {
- d_exit = d_enter;
- d_enter = dist;
- } else {
- d_exit = dist;
- }
+ d_exit = dist;
}
- first++;
}
- } // loop over solutions
- }
- delete[] solutions;
+ first++;
+ }
+ } // loop over solutions
- if (first != 2) { // entry and exit points found
- CORSIKA_LOG_DEBUG("no intersection! count={}", first);
+ if (first == 0) { // entry and exit points found
+ CORSIKA_LOG_DEBUG(
+ "no intersections found: count={}, first_entry={}, first_exit={}", first,
+ first_entry, first_exit);
return Intersections();
}
return Intersections(d_enter / absVelocity, d_exit / absVelocity);
@@ -231,79 +241,80 @@ namespace corsika {
inline Intersections Tracking::intersect(TParticle const& particle,
Plane const& plane) {
- int chargeNumber;
- if (is_nucleus(particle.getPID())) {
- chargeNumber = particle.getNuclearZ();
- } else {
- chargeNumber = get_charge_number(particle.getPID());
- }
- auto const* currentLogicalVolumeNode = particle.getNode();
- VelocityVector const velocity =
- particle.getMomentum() / particle.getEnergy() * constants::c;
- auto const absVelocity = velocity.getNorm();
- DirectionVector const direction =
- velocity.normalized(); // determine steplength to next volume
- Point const position = particle.getPosition();
-
- auto const magneticfield =
- currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
+ CORSIKA_LOG_TRACE("intersection particle with plane");
+
+ ElectricChargeType const charge = particle.getCharge();
- if (chargeNumber != 0 && abs(plane.getNormal().dot(velocity.cross(magneticfield))) >
- 1e-6_T * 1_m / 1_s) {
+ if (charge != 0 * constants::e) {
auto const* currentLogicalVolumeNode = particle.getNode();
+ VelocityVector const velocity = particle.getVelocity();
+ auto const absVelocity = velocity.getNorm();
+ DirectionVector const direction = velocity.normalized();
+ Point const position = particle.getPosition();
+
auto const magneticfield =
currentLogicalVolumeNode->getModelProperties().getMagneticField(position);
- auto const k =
- chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
- auto const direction_B_perp = direction.cross(magneticfield);
- auto const denom = plane.getNormal().dot(direction_B_perp) * k;
- auto const sqrtArg =
- direction.dot(plane.getNormal()) * direction.dot(plane.getNormal()) -
- (plane.getNormal().dot(position - plane.getCenter()) * denom * 2);
+ // solve: denom x^2 + p x + q =0 for x = delta-l
+
+ auto const direction_x_B = direction.cross(magneticfield);
+ double const denom = charge *
+ plane.getNormal().dot(direction_x_B) // unit: C*T = kg/s
+ / 1_kg * 1_s;
+
+ CORSIKA_LOG_TRACE("denom={}", denom);
+
+ auto const p_norm =
+ constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()); // unit: kg * m/s
+
+ double const p = (2 * p_norm * direction.dot(plane.getNormal())) // unit: kg*m/s
+ / (1_m * 1_kg) * 1_s;
+ double const q =
+ (2 * p_norm *
+ plane.getNormal().dot(position - plane.getCenter())) // unit: kg*m/s *m
+ / (1_m * 1_m * 1_kg) * 1_s;
+
+ std::vector<double> deltaLs = solve_quadratic_real(denom, p, q);
- if (sqrtArg < 0) {
+ if (deltaLs.size() == 0) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
}
- double const sqrSqrtArg = sqrt(sqrtArg);
- auto const norm_projected =
- direction.dot(plane.getNormal()) / direction.getNorm();
- LengthType const MaxStepLength1 = (sqrSqrtArg - norm_projected) / denom;
- LengthType const MaxStepLength2 = (-sqrSqrtArg - norm_projected) / denom;
- CORSIKA_LOG_TRACE("MaxStepLength1={}, MaxStepLength2={}", MaxStepLength1,
- MaxStepLength2);
+ // select smallest but positive solution
+ bool first = true;
+ LengthType maxStepLength = 0_m;
+ for (auto& deltaL : deltaLs) {
+ if (deltaL < 0) continue;
+ if (first) {
+ first = false;
+ maxStepLength = deltaL * meter;
+ } else if (maxStepLength > deltaL * meter) {
+ maxStepLength = deltaL * meter;
+ }
+ }
- // check: both intersections in past
- if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
+ // check: both intersections in past, or no valid intersection
+ if (first) {
return Intersections(std::numeric_limits<double>::infinity() * 1_s);
-
- // check: next intersection is MaxStepLength2
- } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {
- CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength2);
- return Intersections(
- MaxStepLength2 *
- (direction + direction_B_perp * MaxStepLength2 * k / 2).getNorm() /
- absVelocity);
-
- // check: next intersections is MaxStepLength1
- } else if (MaxStepLength2 <= 0_m || MaxStepLength1 < MaxStepLength2) {
- CORSIKA_LOG_TRACE(" steplength to obs plane 2: {} ", MaxStepLength1);
- return Intersections(
- MaxStepLength1 *
- (direction + direction_B_perp * MaxStepLength1 * k / 2).getNorm() /
- absVelocity);
}
- CORSIKA_LOG_WARN(
- "Particle wasn't tracked with curved trajectory -> straight (is this an "
- "ERROR?)");
+ CORSIKA_LOG_TRACE("maxStepLength={}", maxStepLength);
+
+ // with final length correction
+ auto const corr =
+ (direction + direction_x_B * maxStepLength * charge / (p_norm * 2))
+ .getNorm() /
+ absVelocity; // unit: s/m
+
+ return Intersections(maxStepLength * corr); // unit: s
- } // end if curved-tracking
+ } // no charge
- CORSIKA_LOG_TRACE("straight tracking with chargeNumber={}, B={}", chargeNumber,
- magneticfield);
+ CORSIKA_LOG_TRACE("(plane) straight tracking with charge={}, B={}", charge,
+ particle.getNode()->getModelProperties().getMagneticField(
+ particle.getPosition()));
return tracking_line::Tracking::intersect(particle, plane);
}
@@ -330,7 +341,7 @@ namespace corsika {
straightTrajectory.getLine().getVelocity(),
MagneticFieldVector(particle.getPosition().getCoordinateSystem(), 0_T, 0_T,
0_T),
- square(0_m) / (square(1_s) * 1_V),
+ 0 * square(meter) / (square(second) * volt),
straightTrajectory.getDuration()), // trajectory
minNode); // next volume node
}
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 61ad3793d..77cef175d 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -50,14 +50,13 @@ namespace corsika {
}
// charge of the particle, and magnetic field
- const int chargeNumber = particle.getChargeNumber();
+ auto const charge = particle.getCharge();
auto magneticfield =
volumeNode->getModelProperties().getMagneticField(initialPosition);
auto const magnitudeB = magneticfield.getNorm();
- CORSIKA_LOG_DEBUG("field={} uT, chargeNumber={}, magnitudeB={} uT",
- magneticfield.getComponents() / 1_uT, chargeNumber,
- magnitudeB / 1_T);
- bool const no_deflection = chargeNumber == 0 || magnitudeB == 0_T;
+ CORSIKA_LOG_DEBUG("field={} uT, charge={}, magnitudeB={} uT",
+ magneticfield.getComponents() / 1_uT, charge, magnitudeB / 1_T);
+ bool const no_deflection = (charge == 0 * constants::e) || magnitudeB == 0_T;
// check, where the first halve-step direction has geometric intersections
auto const [initialTrack, initialTrackNextVolume] =
@@ -75,20 +74,20 @@ namespace corsika {
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
- HEPMomentumType const pAlongB_delta =
+ HEPMomentumType const p_perp =
(particle.getMomentum() -
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
- if (pAlongB_delta == 0_GeV) {
+ if (p_perp == 0_GeV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
- CORSIKA_LOG_TRACE("pAlongB_delta is 0_GeV --> parallel");
+ CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
- LengthType const gyroradius =
- (pAlongB_delta * 1_V / (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
+ constants::c / (abs(charge) * magnitudeB));
// we need to limit maximum step-length since we absolutely
// need to follow strongly curved trajectories segment-wise,
@@ -111,9 +110,9 @@ namespace corsika {
firstHalveSteplength, steplimit, initialTrackLength);
// perform the first halve-step
Point const position_mid = initialPosition + direction * firstHalveSteplength;
- auto const k =
- chargeNumber * (constants::c * 1_eV / 1_V) / particle.getMomentum().getNorm();
- auto const new_direction =
+ auto const k = charge / (constants::c * convert_HEP_to_SI<MassType::dimension_type>(
+ particle.getMomentum().getNorm()));
+ DirectionVector const new_direction =
direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
auto const new_direction_norm = new_direction.getNorm(); // by design this is >1
CORSIKA_LOG_DEBUG(
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index 9a4659250..36e097560 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -20,6 +20,7 @@
#include <type_traits>
#include <utility>
+#include <cmath>
namespace corsika::tracking_line {
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index ee29642e0..f251c60ba 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -268,10 +268,8 @@ namespace corsika::urqmd {
2 * ::urqmd::options_.CTParam[30 - 1];
::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV);
- if (projectileCode == Code::K0Long) {
+ if (projectileCode == Code::K0Long || projectileCode == Code::K0Short) {
projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar;
- } else if (projectileCode == Code::K0Short) {
- throw std::runtime_error("K0Short should not interact");
}
auto const [ityp, iso3] = convertToUrQMD(projectileCode);
diff --git a/corsika/detail/stack/NuclearStackExtension.inl b/corsika/detail/stack/NuclearStackExtension.inl
index 728647ca0..5dab28ab2 100644
--- a/corsika/detail/stack/NuclearStackExtension.inl
+++ b/corsika/detail/stack/NuclearStackExtension.inl
@@ -113,6 +113,14 @@ namespace corsika::nuclear_stack {
return super_type::getMass();
}
+ template <template <typename> class InnerParticleInterface,
+ typename StackIteratorInterface>
+ inline ElectricChargeType NuclearParticleInterface<
+ InnerParticleInterface, StackIteratorInterface>::getCharge() const {
+ if (super_type::getPID() == Code::Nucleus) return getNuclearZ() * constants::e;
+ return super_type::getCharge();
+ }
+
template <template <typename> class InnerParticleInterface,
typename StackIteratorInterface>
inline int16_t NuclearParticleInterface<
diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp
index 7ee12186a..b9ee5361b 100644
--- a/corsika/framework/core/Cascade.hpp
+++ b/corsika/framework/core/Cascade.hpp
@@ -130,8 +130,9 @@ namespace corsika {
*/
void step(Particle& vParticle);
- ProcessReturn decay(TStackView& view);
- ProcessReturn interaction(TStackView& view);
+ ProcessReturn decay(TStackView& view, InverseTimeType initial_inv_decay_time);
+ ProcessReturn interaction(TStackView& view,
+ InverseGrammageType initial_inv_int_length);
void setEventType(TStackView& view, history::EventType);
// data members
diff --git a/corsika/framework/geometry/LeapFrogTrajectory.hpp b/corsika/framework/geometry/LeapFrogTrajectory.hpp
index 6b67bf634..11e3534a1 100644
--- a/corsika/framework/geometry/LeapFrogTrajectory.hpp
+++ b/corsika/framework/geometry/LeapFrogTrajectory.hpp
@@ -39,7 +39,7 @@ namespace corsika {
LeapFrogTrajectory(Point const& pos, VelocityVector const& initialVelocity,
MagneticFieldVector const& Bfield,
- decltype(square(meter) / (square(second) * volt)) const k,
+ decltype(1 / (tesla * second)) const k,
TimeType const timeStep) // leap-from total length
: initialPosition_(pos)
, initialVelocity_(initialVelocity)
@@ -74,7 +74,7 @@ namespace corsika {
VelocityVector initialVelocity_;
DirectionVector initialDirection_;
MagneticFieldVector magneticfield_;
- decltype(square(meter) / (square(second) * volt)) k_;
+ decltype(1 / (tesla * second)) k_;
TimeType timeStep_;
};
diff --git a/corsika/framework/utility/CubicSolver.hpp b/corsika/framework/utility/CubicSolver.hpp
new file mode 100644
index 000000000..d62757277
--- /dev/null
+++ b/corsika/framework/utility/CubicSolver.hpp
@@ -0,0 +1,88 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <numeric>
+#include <complex>
+
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+#include <boost/multiprecision/cpp_bin_float.hpp>
+
+/**
+ @file CubicSolver.hpp
+ */
+
+namespace corsika {
+
+ /**
+ @ingroup Utility
+ @{
+ */
+
+ namespace andre {
+ /**
+ Solve a x^3 + b x^2 + c x + d = 0
+ */
+ std::vector<double> solveP3(long double a, long double b, long double c,
+ long double d, double const epsilon = 1e-12);
+ } // namespace andre
+
+ /**
+ Solve depressed cubic: x^3 + p x + q = 0
+ */
+ inline std::vector<double> solve_cubic_depressed_disciminant_real(
+ long double p, long double q, long double const disc, double const epsilon = 1e-12);
+
+ std::vector<double> solve_cubic_depressed_real(long double p, long double q,
+ double const epsilon = 1e-12);
+
+ /**
+ Solve a x^3 + b x^2 + c x + d = 0
+
+ Analytical approach. Not very stable in all conditions.
+ */
+ std::vector<double> solve_cubic_real_analytic(long double a, long double b,
+ long double c, long double d,
+ double const epsilon = 1e-12);
+
+ /**
+ @return a x^3 + b x^2 + c x + d
+ */
+ long double cubic_function(long double x, long double a, long double b, long double c,
+ long double d);
+
+ /**
+ @return 3 a x^2 + 2 b x + c
+ */
+ long double cubic_function_dfdx(long double x, long double a, long double b,
+ long double c);
+
+ /**
+ @return 6 a x + 2 b
+ */
+ long double cubic_function_d2fd2x(long double x, long double a, long double b);
+
+ /**
+ Iterative approach to solve: a x^3 + b x^2 + c x + d = 0
+
+ This is often fastest and most precise.
+ */
+ std::vector<double> solve_cubic_real(long double a, long double b, long double c,
+ long double d, double const epsilon = 1e-12);
+
+ //! @}
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/CubicSolver.inl>
diff --git a/corsika/framework/utility/LinearSolver.hpp b/corsika/framework/utility/LinearSolver.hpp
new file mode 100644
index 000000000..95322958d
--- /dev/null
+++ b/corsika/framework/utility/LinearSolver.hpp
@@ -0,0 +1,24 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+namespace corsika {
+
+ std::vector<double> solve_linear_real(double a, double b, double const epsilon = 1e-12);
+
+ std::vector<std::complex<double>> solve_linear(double a, double b,
+ double const epsilon = 1e-12);
+
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/LinearSolver.inl>
diff --git a/corsika/framework/utility/QuadraticSolver.hpp b/corsika/framework/utility/QuadraticSolver.hpp
new file mode 100644
index 000000000..968c01bf1
--- /dev/null
+++ b/corsika/framework/utility/QuadraticSolver.hpp
@@ -0,0 +1,27 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include <complex>
+
+#include <corsika/framework/utility/LinearSolver.hpp>
+
+namespace corsika {
+
+ std::vector<std::complex<double>> solve_quadratic(long double a, long double b,
+ long double c,
+ double const epsilon = 1e-12);
+
+ std::vector<double> solve_quadratic_real(long double a, long double b, long double c,
+ double const epsilon = 1e-12);
+} // namespace corsika
+
+#include <corsika/detail/framework/utility/QuadraticSolver.inl>
diff --git a/corsika/framework/utility/QuarticSolver.hpp b/corsika/framework/utility/QuarticSolver.hpp
index ad9484939..d6e20f02f 100644
--- a/corsika/framework/utility/QuarticSolver.hpp
+++ b/corsika/framework/utility/QuarticSolver.hpp
@@ -8,59 +8,46 @@
#pragma once
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+
#include <complex>
+#include <algorithm>
+#include <cmath>
+#include <numeric>
/**
* @file QuarticSolver.hpp
- * @ingroup Utilities
- * @{
- * \todo convert to class
*/
-namespace corsika::quartic_solver {
-
- const double PI = 3.141592653589793238463L;
- const double M_2PI = 2 * PI;
- const double eps = 1e-12;
-
- typedef std::complex<double> DComplex;
-
- //---------------------------------------------------------------------------
- // useful for testing
- inline DComplex polinom_2(DComplex x, double a, double b) {
- // Horner's scheme for x*x + a*x + b
- return x * (x + a) + b;
- }
+namespace corsika {
- //---------------------------------------------------------------------------
- // useful for testing
- inline DComplex polinom_3(DComplex x, double a, double b, double c) {
- // Horner's scheme for x*x*x + a*x*x + b*x + c;
- return x * (x * (x + a) + b) + c;
- }
+ /**
+ * @ingroup Utilities
+ * @{
+ */
- //---------------------------------------------------------------------------
- // useful for testing
- inline DComplex polinom_4(DComplex x, double a, double b, double c, double d) {
- // Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
- return x * (x * (x * (x + a) + b) + c) + d;
- }
+ namespace andre {
- //---------------------------------------------------------------------------
- // x - array of size 3
- // In case 3 real roots: => x[0], x[1], x[2], return 3
- // 2 real roots: x[0], x[1], return 2
- // 1 real root : x[0], x[1] ± i*x[2], return 1
- unsigned int solveP3(double* x, double a, double b, double c);
+ /**
+ solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e
+ Attention - this function returns dynamically allocated array. It has to be
+ released afterwards.
+ */
+ std::vector<double> solve_quartic_real(long double a, long double b, long double c,
+ long double d, long double e,
+ double const epsilon = 1e-12);
+ } // namespace andre
- //---------------------------------------------------------------------------
- // solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
- // Attention - this function returns dynamically allocated array. It has to be released
- // afterwards.
- DComplex* solve_quartic(double a, double b, double c, double d);
+ /**
+ solve quartic equation a*x^4 + b*x^3 + c*x^2 + d*x + e
+ */
+ std::vector<double> solve_quartic_real(long double a, long double b, long double c,
+ long double d, long double e,
+ double const epsilon = 1e-12);
//! @}
-} // namespace corsika::quartic_solver
+} // namespace corsika
#include <corsika/detail/framework/utility/QuarticSolver.inl>
diff --git a/corsika/modules/ObservationPlane.hpp b/corsika/modules/ObservationPlane.hpp
index 8276c0d12..f7ccba999 100644
--- a/corsika/modules/ObservationPlane.hpp
+++ b/corsika/modules/ObservationPlane.hpp
@@ -20,9 +20,19 @@
namespace corsika {
/**
- * The ObservationPlane writes PDG codes, energies, and distances of particles to the
- * central point of the plane into its output file. The particles are considered
- * "absorbed" afterwards.
+ @ingroup Modules
+ @{
+
+ The ObservationPlane writes PDG codes, energies, and distances of particles to the
+ central point of the plane into its output file. The particles are considered
+ "absorbed" afterwards.
+
+ **Note/Limitation:** as discussed in
+ https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/397
+ you cannot put two ObservationPlanes exactly on top of each
+ other. Even if one of them is "permeable". You have to put a
+ small gap in between the two plane in such a scenario, or develop
+ another more specialized output class.
*/
class ObservationPlane : public ContinuousProcess<ObservationPlane> {
@@ -50,6 +60,7 @@ namespace corsika {
DirectionVector const xAxis_;
DirectionVector const yAxis_;
};
+ //! @}
} // namespace corsika
#include <corsika/detail/modules/ObservationPlane.inl>
diff --git a/corsika/stack/NuclearStackExtension.hpp b/corsika/stack/NuclearStackExtension.hpp
index 34bd1377d..19e9a44c3 100644
--- a/corsika/stack/NuclearStackExtension.hpp
+++ b/corsika/stack/NuclearStackExtension.hpp
@@ -92,6 +92,10 @@ namespace corsika::nuclear_stack {
* Overwrite normal getParticleMass function with nuclear version
*/
HEPMassType getMass() const;
+ /**
+ * Overwrite normal getParticleCharge function with nuclear version
+ */
+ ElectricChargeType getCharge() const;
/**
* Overwirte normal getChargeNumber function with nuclear version
**/
diff --git a/corsika/stack/VectorStack.hpp b/corsika/stack/VectorStack.hpp
index adaae996b..e5a765781 100644
--- a/corsika/stack/VectorStack.hpp
+++ b/corsika/stack/VectorStack.hpp
@@ -87,8 +87,14 @@ namespace corsika {
return this->getMomentum() / this->getEnergy();
}
+ VelocityVector getVelocity() const {
+ return this->getMomentum() / this->getEnergy() * constants::c;
+ }
+
HEPMassType getMass() const { return get_mass(this->getPID()); }
+ ElectricChargeType getCharge() const { return get_charge(this->getPID()); }
+
int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
///@}
};
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index e5f6455d3..1b60722d8 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -145,7 +145,7 @@ int main(int argc, char** argv) {
"environment setup: universe={}, layer1={}, layer2={}, layer3={}, layer4={}, "
"layer5={}",
fmt::ptr(env.getUniverse()->getContainingNode(
- Point(rootCS, {constants::EarthRadius::Mean + 130_km, 0_m, 0_m}))),
+ Point(rootCS, {constants::EarthRadius::Mean + 130000_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
Point(rootCS, {constants::EarthRadius::Mean + 110_km, 0_m, 0_m}))),
fmt::ptr(env.getUniverse()->getContainingNode(
@@ -289,6 +289,46 @@ int main(int argc, char** argv) {
// setup particle stack, and add primary particle
setup::Stack stack;
stack.clear();
+ unsigned short const A = std::stoi(std::string(argv[1]));
+ Code beamCode;
+ HEPEnergyType mass;
+ unsigned short Z = 0;
+ if (A > 0) {
+ beamCode = Code::Nucleus;
+ Z = std::stoi(std::string(argv[2]));
+ mass = get_nucleus_mass(A, Z);
+ } else {
+ int pdg = std::stoi(std::string(argv[2]));
+ beamCode = convert_from_PDG(PDGCode(pdg));
+ mass = get_mass(beamCode);
+ }
+ HEPEnergyType const E0 = 1_GeV * std::stof(std::string(argv[3]));
+ double theta = 0.;
+ auto const thetaRad = theta / 180. * M_PI;
+
+ auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
+ return sqrt((Elab - m) * (Elab + m));
+ };
+ HEPMomentumType P0 = elab2plab(E0, mass);
+ auto momentumComponents = [](double thetaRad, HEPMomentumType ptot) {
+ return std::make_tuple(ptot * sin(thetaRad), 0_eV, -ptot * cos(thetaRad));
+ };
+
+ auto const [px, py, pz] = momentumComponents(thetaRad, P0);
+ auto plab = MomentumVector(rootCS, {px, py, pz});
+ cout << "input particle: " << beamCode << endl;
+ cout << "input angles: theta=" << theta << endl;
+ cout << "input momentum: " << plab.getComponents() / 1_GeV
+ << ", norm = " << plab.getNorm() << endl;
+
+ auto const observationHeight = 0_km + builder.getEarthRadius();
+ auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+ auto const t = (injectionHeight - observationHeight) / cos(thetaRad);
+ Point const showerCore{rootCS, 0_m, 0_m, observationHeight};
+ Point const injectionPos =
+ showerCore + DirectionVector{rootCS, {-sin(thetaRad), 0, cos(thetaRad)}} * t;
+
+ std::cout << "point of injection: " << injectionPos.getCoordinates() << std::endl;
if (A > 1) {
stack.addParticle(std::make_tuple(beamCode, E0, plab, injectionPos, 0_ns, A, Z));
diff --git a/tests/framework/CMakeLists.txt b/tests/framework/CMakeLists.txt
index 35f61737b..829709779 100644
--- a/tests/framework/CMakeLists.txt
+++ b/tests/framework/CMakeLists.txt
@@ -20,6 +20,7 @@ set (test_framework_sources
testNullModel.cpp
testHelix.cpp
testInteractionCounter.cpp
+ testSolver.cpp
)
CORSIKA_ADD_TEST (testFramework SOURCES ${test_framework_sources})
diff --git a/tests/framework/testCOMBoost.cpp b/tests/framework/testCOMBoost.cpp
index caf6f0896..d092a4b61 100644
--- a/tests/framework/testCOMBoost.cpp
+++ b/tests/framework/testCOMBoost.cpp
@@ -40,7 +40,6 @@ auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
TEST_CASE("rotation") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// define projectile kinematics in lab frame
HEPMassType const projectileMass = 1_GeV;
@@ -172,7 +171,6 @@ TEST_CASE("rotation") {
TEST_CASE("boosts") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
@@ -329,7 +327,6 @@ TEST_CASE("boosts") {
TEST_CASE("rest frame") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
HEPMassType const projectileMass = 1_GeV;
HEPMomentumType const P0 = 1_TeV;
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index b860e9955..3e930911b 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -137,7 +137,6 @@ public:
TEST_CASE("Cascade", "[Cascade]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
HEPEnergyType E0 = 100_GeV;
diff --git a/tests/framework/testClassTimer.cpp b/tests/framework/testClassTimer.cpp
index d51a56b1a..b88d2fc97 100644
--- a/tests/framework/testClassTimer.cpp
+++ b/tests/framework/testClassTimer.cpp
@@ -106,7 +106,6 @@ public:
TEST_CASE("ClassTimer", "[Timer]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Measure runtime of a function without arguments") {
diff --git a/tests/framework/testCombinedStack.cpp b/tests/framework/testCombinedStack.cpp
index 1f76fafc2..de6defc6a 100644
--- a/tests/framework/testCombinedStack.cpp
+++ b/tests/framework/testCombinedStack.cpp
@@ -89,7 +89,6 @@ using StackTest = CombinedStack<TestStackData, TestStackData2, CombinedTestInter
TEST_CASE("Combined Stack", "[stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
@@ -395,7 +394,6 @@ using Particle2 = typename StackTest2::particle_type;
TEST_CASE("Combined Stack - secondary view") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("create secondaries via secondaryview") {
diff --git a/tests/framework/testCorsikaFenv.cpp b/tests/framework/testCorsikaFenv.cpp
index 2b6cdb580..030d33dfa 100644
--- a/tests/framework/testCorsikaFenv.cpp
+++ b/tests/framework/testCorsikaFenv.cpp
@@ -22,7 +22,6 @@ static void handle_fpe(int /*signo*/) { gRESULT = 0; }
TEST_CASE("CorsikaFenv", "[fenv]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Enable all exceptions") { feenableexcept(FE_ALL_EXCEPT); }
diff --git a/tests/framework/testFourVector.cpp b/tests/framework/testFourVector.cpp
index 09b458295..489d6261b 100644
--- a/tests/framework/testFourVector.cpp
+++ b/tests/framework/testFourVector.cpp
@@ -20,7 +20,6 @@ using namespace corsika;
TEST_CASE("four vectors") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
// this is just needed as a baseline
diff --git a/tests/framework/testFunctionTimer.cpp b/tests/framework/testFunctionTimer.cpp
index 46d97f710..1328ea7ba 100644
--- a/tests/framework/testFunctionTimer.cpp
+++ b/tests/framework/testFunctionTimer.cpp
@@ -32,7 +32,6 @@ public:
TEST_CASE("FunctionTimer", "[Timer]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Measure runtime of a free function") {
diff --git a/tests/framework/testHelix.cpp b/tests/framework/testHelix.cpp
index cbe997d20..042e4212f 100644
--- a/tests/framework/testHelix.cpp
+++ b/tests/framework/testHelix.cpp
@@ -22,7 +22,6 @@ double constexpr absMargin = 1.0e-8;
TEST_CASE("Helix class") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
const CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
Point r0(rootCS, {0_m, 0_m, 0_m});
diff --git a/tests/framework/testInteractionCounter.cpp b/tests/framework/testInteractionCounter.cpp
index 556d4733b..4b4b7ff25 100644
--- a/tests/framework/testInteractionCounter.cpp
+++ b/tests/framework/testInteractionCounter.cpp
@@ -43,7 +43,6 @@ struct DummyProcess {
TEST_CASE("InteractionCounter", "[process]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::info);
DummyProcess d;
diff --git a/tests/framework/testLogging.cpp b/tests/framework/testLogging.cpp
index b737d7dfa..92ad036e5 100644
--- a/tests/framework/testLogging.cpp
+++ b/tests/framework/testLogging.cpp
@@ -14,7 +14,7 @@ using namespace corsika;
TEST_CASE("Logging", "[Logging]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+ // corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("top level functions using default corsika logger") {
logging::info("(1) This is an info message!");
diff --git a/tests/framework/testNullModel.cpp b/tests/framework/testNullModel.cpp
index 5b1f351e1..872fc34db 100644
--- a/tests/framework/testNullModel.cpp
+++ b/tests/framework/testNullModel.cpp
@@ -21,7 +21,6 @@ using namespace corsika;
TEST_CASE("NullModel", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("interface") {
[[maybe_unused]] NullModel nm;
diff --git a/tests/framework/testParticles.cpp b/tests/framework/testParticles.cpp
index b74d836bc..d8c71e6e7 100644
--- a/tests/framework/testParticles.cpp
+++ b/tests/framework/testParticles.cpp
@@ -17,7 +17,6 @@ using namespace corsika;
TEST_CASE("ParticleProperties", "[Particles]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Types") {
CHECK(Electron::code == Code::Electron);
diff --git a/tests/framework/testRandom.cpp b/tests/framework/testRandom.cpp
index 5774575e6..2d0c0720a 100644
--- a/tests/framework/testRandom.cpp
+++ b/tests/framework/testRandom.cpp
@@ -55,7 +55,6 @@ SCENARIO("random-number streams can be registered and retrieved") {
TEST_CASE("UniformRealDistribution") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::mt19937 rng;
@@ -101,7 +100,6 @@ TEST_CASE("UniformRealDistribution") {
TEST_CASE("ExponentialDistribution") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
std::mt19937 rng;
diff --git a/tests/framework/testSaveBoostHistogram.cpp b/tests/framework/testSaveBoostHistogram.cpp
index d5c16438a..89182311a 100644
--- a/tests/framework/testSaveBoostHistogram.cpp
+++ b/tests/framework/testSaveBoostHistogram.cpp
@@ -16,7 +16,6 @@ namespace bh = boost::histogram;
TEST_CASE("SaveHistogram") {
- corsika::corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
corsika::logging::set_level(corsika::logging::level::info);
std::mt19937 rng;
diff --git a/tests/framework/testSecondaryView.cpp b/tests/framework/testSecondaryView.cpp
index 87cfd91c0..fa113f88b 100644
--- a/tests/framework/testSecondaryView.cpp
+++ b/tests/framework/testSecondaryView.cpp
@@ -44,7 +44,6 @@ using Particle = typename StackTest::particle_type;
TEST_CASE("SecondaryStack", "[stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
diff --git a/tests/framework/testSolver.cpp b/tests/framework/testSolver.cpp
new file mode 100644
index 000000000..52f1fa13f
--- /dev/null
+++ b/tests/framework/testSolver.cpp
@@ -0,0 +1,278 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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.
+ */
+
+#include <catch2/catch.hpp>
+
+#include <cmath>
+#include <vector>
+#include <complex>
+#include <algorithm>
+
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/QuadraticSolver.hpp>
+#include <corsika/framework/utility/LinearSolver.hpp>
+#include <corsika/framework/utility/CubicSolver.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
+
+using namespace corsika;
+using namespace std;
+
+double pol4(long double x, long double a, long double b, long double c, long double d,
+ long double e) {
+ return std::pow(x, 4) * a + std::pow(x, 3) * b + std::pow(x, 2) * c + x * d + e;
+}
+
+void remove_duplicates(vector<double>& v, double const eps) {
+ std::sort(v.begin(), v.end());
+ v.erase(std::unique(v.begin(), v.end(),
+ [eps](auto v1, auto v2) {
+ return (std::abs(v2) > eps ? std::abs(2 * (v1 - v2) / (v1 + v2)) <
+ eps // relative
+ : std::abs(v1 - v2) < eps); // absolute
+ }),
+ v.end());
+}
+
+TEST_CASE("Solver") {
+
+ logging::set_level(logging::level::trace);
+
+ double epsilon_check = 1e-3; // for catch2 asserts
+
+ SECTION("linear") {
+
+ std::vector<std::pair<double, double>> vals = {{13.33, 8338.3},
+ {-333.99, -8.15},
+ {-58633.3, 2343.54},
+ {87632.87, -982.37},
+ {7e8, -1e-7}};
+
+ for (auto v : vals) {
+
+ {
+ double a = v.first;
+ double b = v.second;
+
+ vector<double> s1 = solve_linear_real(a, b);
+
+ CORSIKA_LOG_INFO("linear: a={}, b={}, N={}, s1[0]={}", a, b, s1.size(), s1[0]);
+
+ double expected = -b / a;
+
+ CHECK(s1.size() == 1);
+ CHECK((s1[0] == Approx(expected).epsilon(epsilon_check)));
+
+ vector<complex<double>> s2 = solve_linear(a, b);
+ CHECK(s2.size() == 1);
+ CHECK((s2[0].real() == Approx(expected).epsilon(epsilon_check)));
+ }
+ }
+
+ CHECK(solve_linear_real(0, 55.).size() == 0);
+ }
+
+ SECTION("quadratic") {
+
+ // tests of type: (x-z1)*(x-z2) = 0 --> x1=z1, x2=z2 and a=1, b=-z2-z1,
+ // c=z1*z2
+
+ std::vector<std::vector<double>> zeros = {
+ {13.33, 8338.3}, {-333.99, -8.15}, {-58633.3, 2343.54}, {87632.87, -982.37},
+ {1.3e-6, 5e7}, {-4.2e-7, 65e6}, {0.1, -13e-5}, {7e8, -1e-7}};
+
+ for (unsigned int idegree = 0; idegree < 2; ++idegree) {
+ CORSIKA_LOG_INFO("------------- quadratic idegree={}", idegree);
+ for (auto z : zeros) {
+
+ {
+ long double const z1 = z[0];
+ long double const z2 = idegree <= 0 ? z1 : z[1];
+
+ long double a = 1;
+ long double b = -z1 - z2;
+ long double c = z1 * z2;
+
+ std::vector<double> s1 = solve_quadratic_real(a, b, c);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("quadratic: z1={}, z2={}, N={}, s1=[{}]", z1, z2, s1.size(),
+ fmt::join(s1, ", "));
+
+ CHECK(s1.size() == idegree + 1);
+ for (auto value : s1) {
+ if (std::abs(value) < epsilon_check) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check))));
+ }
+ }
+ /*
+ std::vector<complex<double>> s2 = solve_quadratic(a, b, c, epsilon);
+ CHECK(s2.size() == idegree + 1);
+ for (auto value : s2) {
+ if (std::abs(value) < epsilon) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check))));
+ }
+ }
+ }*/
+ }
+ }
+ }
+ }
+
+ SECTION("cubic") {
+
+ // tests of type:
+ // (x-z1) * (x-z2) * (x-z3) = 0 --> x1=z1, x2=z2, x3=z3 and
+
+ // (x^2 - x*(z1+z2) + z1*z2) * (x-z3) =
+ // x^3 - x^2*(z1+z2) + x*z1*z2 - x^2 z3 + x*(z1+z2)*z3 - z1*z2*z3 =
+ // x^3 + x^2 (-z3-z1-z2) + x (z1*z2 + (z1+z2)*z3) - z1*z2*z3
+
+ epsilon_check = 1e-2; // for catch2 asserts
+
+ std::vector<std::vector<double>> zeros = {
+ {13.33, 838.3, 44.}, {-333.99, -8.15, -33e8}, {-58633.3, 2343.54, -1e-5},
+ {87632.87, -982.37, 0.}, {1.3e-4, 5e7, 6.6e9}, {-4.2e-7, 65e6, 433.3334},
+ {23e-1, -13e-2, 10.}, {7e8, -1e-7, 1e8}};
+
+ for (unsigned int idegree = 0; idegree < 3; ++idegree) {
+ CORSIKA_LOG_INFO("------------- cubic idegree={}", idegree);
+ for (auto z : zeros) {
+
+ {
+ long double const z1 = z[0];
+ long double const z2 = idegree <= 0 ? z1 : z[1];
+ long double const z3 = idegree <= 1 ? z1 : z[2];
+
+ long double const a = 1;
+ long double const b = -z1 - z2 - z3;
+ long double const c = z1 * z2 + (z1 + z2) * z3;
+ long double const d = -z1 * z2 * z3;
+ //
+ CORSIKA_LOG_INFO(
+ "cubic: z1={}, z2={}, z3={}, "
+ "a={}, b={}, c={}, d={}",
+ z1, z2, z3, a, b, c, d);
+
+ vector<double> s1 = solve_cubic_real(a, b, c, d);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
+
+ CHECK(s1.size() == idegree + 1);
+ for (double value : s1) {
+ CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} eps_check={}", value, z1, z2,
+ z3, epsilon_check);
+ if (std::abs(value) < epsilon_check) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check)) ||
+ (value == Approx(z3).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check)) ||
+ (value == Approx(z3).epsilon(epsilon_check))));
+ }
+ }
+ }
+ }
+ }
+ } // cubic
+
+ SECTION("quartic") {
+
+ epsilon_check = 1e-4; // for catch2 asserts
+
+ // **clang-format-off**
+ // tests of type:
+ // (x-z1) * (x-z2) * (x-z3) * (x-z4) = 0 --> x1=z1, x2=z2, x3=z3, x4=z4 and
+
+ // (x^2 - x (z1+z2) + z1*z2) * (x-z3) *
+ // (x-z4) = (x^3 - x^2 (z1+z2) + x z1*z2 - x^2 z3 + x*(z1+z2)*z3 -
+ // z1*z2*z3) * (x-z4) = (x^3 + x^2 (-z1-z2-z3) + x (z1*z2 + (z1+z2)*z3) -
+ // z1*z2*z3) * (x-z4)
+ // =
+ //
+ // x^4 + x^3 (-z1-z2-z3) + x^2 (z1*z2 + (z1+z2)*z3) - x z1*z2*z3
+ // - x^3 z4 - x^2 (-z1-z2-z3)*z4 - x (z1*z2 + (z1+z2)*z3)*z4
+ // + z1*z2*z3*z4
+ //
+ // = x^4
+ // + x^3 (-z1-z2-z3-z4)
+ // + x^2 (z1*z2 + (z1+z2)*z3 + (z1+z2+z3)*z4)
+ // - x (z1*z2*z3 + (z1*z2 + (z1+z2)*z3)*z4))
+ // + z1*z2*z3*z4
+ // **clang-format-on**
+
+ std::vector<std::vector<double>> zeros = {
+ {13.33, 838.3, 44., 2.3}, {-3333.99, -8.15, -33e4, 8.8e3},
+ {-58633.3, 2343.54, -1e-1, 0.}, {87632.87, -982.37, 10., 1e-2},
+ {1.3e2, 5e5, 6.6e5, 1e3}, {-4.9, 65e2, 433.3334, 6e5},
+ {23e-1, -13e-2, 10., 3.4e6}, {7e6, -1e-1, 1e5, 2.55e4}};
+
+ for (unsigned int idegree = 2; idegree < 4;
+ ++idegree) { // idegree==1 is not very stable !!!!!!!
+ CORSIKA_LOG_INFO("------------- quartic idegree={}", idegree);
+ for (auto z : zeros) {
+
+ {
+ long double const z1 = z[0];
+ long double const z2 = idegree <= 0 ? z1 : z[1];
+ long double const z3 = idegree <= 1 ? z1 : z[2];
+ long double const z4 = idegree <= 2 ? z1 : z[3];
+
+ long double const a = 1;
+ long double const b = -z1 - z2 - z3 - z4;
+ long double const c = z1 * z2 + (z1 + z2) * z3 + (z1 + z2 + z3) * z4;
+ long double const d = -z1 * z2 * z3 - (z1 * z2 + (z1 + z2) * z3) * z4;
+ long double const e = z1 * z2 * z3 * z4;
+
+ //
+ CORSIKA_LOG_INFO(
+ "quartic: z1={}, z2={}, z3={}, z4={}, "
+ "a={}, b={}, c={}, d={}, e={}",
+ z1, z2, z3, z4, a, b, c, d, e);
+
+ // make sure the tested cases are all ZERO (printout)
+ CORSIKA_LOG_INFO("quartic trace: {} {} {} {}", pol4(z1, a, b, c, d, e),
+ pol4(z2, a, b, c, d, e), pol4(z3, a, b, c, d, e),
+ pol4(z4, a, b, c, d, e));
+
+ vector<double> s1 = andre::solve_quartic_real(a, b, c, d, e);
+ // vector<double> s1 = solve_quartic_real(a, b, c, d, e, epsilon);
+ remove_duplicates(s1, epsilon_check);
+
+ CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
+
+ CHECK(s1.size() == idegree + 1);
+ for (double value : s1) {
+ CORSIKA_LOG_INFO("value={}, z1={} z2={} z3={} z4={} eps_check={}", value, z1,
+ z2, z3, z4, epsilon_check);
+ if (std::abs(value) < epsilon_check) {
+ CHECK(((value == Approx(z1).margin(epsilon_check)) ||
+ (value == Approx(z2).margin(epsilon_check)) ||
+ (value == Approx(z3).margin(epsilon_check)) ||
+ (value == Approx(z4).margin(epsilon_check))));
+ } else {
+ CHECK(((value == Approx(z1).epsilon(epsilon_check)) ||
+ (value == Approx(z2).epsilon(epsilon_check)) ||
+ (value == Approx(z3).epsilon(epsilon_check)) ||
+ (value == Approx(z4).epsilon(epsilon_check))));
+ }
+ }
+ }
+ }
+ }
+ } // quartic
+}
diff --git a/tests/framework/testStackInterface.cpp b/tests/framework/testStackInterface.cpp
index b37dc495d..4171647af 100644
--- a/tests/framework/testStackInterface.cpp
+++ b/tests/framework/testStackInterface.cpp
@@ -27,7 +27,6 @@ typedef Stack<TestStackData, TestParticleInterface> StackTest;
TEST_CASE("Stack", "[Stack]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
// helper function for sum over stack data
auto sum = [](const StackTest& stack) {
diff --git a/tests/framework/testUnits.cpp b/tests/framework/testUnits.cpp
index 8d15544ac..5c90aa9cd 100644
--- a/tests/framework/testUnits.cpp
+++ b/tests/framework/testUnits.cpp
@@ -19,7 +19,6 @@ using namespace corsika;
TEST_CASE("PhysicalUnits", "[Units]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Consistency") {
CHECK(1_m / 1_m == Approx(1));
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 34718d69b..2128e0a47 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -39,7 +39,6 @@ using MExtraEnvirnoment = MediumPropertyModel<UniformMagneticField<T>>;
TEST_CASE("CONEXSourceCut") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
RNGManager::getInstance().registerRandomStream("cascade");
RNGManager::getInstance().registerRandomStream("sibyll");
@@ -123,7 +122,6 @@ TEST_CASE("CONEXSourceCut") {
TEST_CASE("ConexOutput", "[output validation]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto file = GENERATE(as<std::string>{}, "conex_fit", "conex_output");
diff --git a/tests/modules/testExecTime.cpp b/tests/modules/testExecTime.cpp
index 68a85d81b..4ebc03342 100644
--- a/tests/modules/testExecTime.cpp
+++ b/tests/modules/testExecTime.cpp
@@ -27,7 +27,6 @@ using namespace corsika::process::example_processors;
TEST_CASE("Timing process", "[proccesses][analytic_processors ExecTime]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
int tmp = 0;
diff --git a/tests/modules/testNullModel.cpp b/tests/modules/testNullModel.cpp
index a24433b08..6485b7df9 100644
--- a/tests/modules/testNullModel.cpp
+++ b/tests/modules/testNullModel.cpp
@@ -14,7 +14,6 @@ using namespace corsika;
TEST_CASE("NullModel", "[processes]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
SECTION("interface") {
diff --git a/tests/modules/testOnShellCheck.cpp b/tests/modules/testOnShellCheck.cpp
index eccf73f28..e28c5ffcd 100644
--- a/tests/modules/testOnShellCheck.cpp
+++ b/tests/modules/testOnShellCheck.cpp
@@ -25,7 +25,6 @@ using namespace corsika;
TEST_CASE("OnShellCheck", "[processes]") {
logging::set_level(logging::level::debug);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
feenableexcept(FE_INVALID);
using EnvType = setup::Environment;
diff --git a/tests/modules/testQGSJetII.cpp b/tests/modules/testQGSJetII.cpp
index 4a604fd8b..3efea7742 100644
--- a/tests/modules/testQGSJetII.cpp
+++ b/tests/modules/testQGSJetII.cpp
@@ -49,7 +49,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("CORSIKA_DATA", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("check CORSIKA_DATA") {
@@ -68,7 +67,6 @@ TEST_CASE("CORSIKA_DATA", "[processes]") {
TEST_CASE("QgsjetII", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("Corsika -> QgsjetII") {
CHECK(corsika::qgsjetII::convertToQgsjetII(PiMinus::code) ==
@@ -133,7 +131,6 @@ TEST_CASE("QgsjetII", "[processes]") {
TEST_CASE("QgsjetIIInterface", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
[[maybe_unused]] auto const& env_dummy = env;
@@ -184,7 +181,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
corsika::qgsjetII::Interaction model;
model.doInteraction(view); // this also should produce some fragments
- CHECK(view.getSize() == Approx(188).margin(2)); // this is not physics validation
+ CHECK(view.getSize() == Approx(228).margin(10)); // this is not physics validation
int countFragments = 0;
for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation
@@ -240,7 +237,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
setup::StackView& view = *(secViewPtr.get());
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
- CHECK(view.getSize() == Approx(7).margin(2)); // this is not physics validation
+ CHECK(view.getSize() == Approx(12).margin(4)); // this is not physics validation
}
{ // Lambda is internally converted into neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
@@ -249,11 +246,11 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
setup::StackView& view = *(secViewPtr.get());
corsika::qgsjetII::Interaction model;
model.doInteraction(view);
- CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+ CHECK(view.getSize() == Approx(10).margin(3)); // this is not physics validation
}
{ // AntiLambda is internally converted into anti neutron
auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::Lambda0Bar, 0, 0, 100_GeV,
+ Code::Lambda0Bar, 0, 0, 1000_GeV,
(setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
setup::StackView& view = *(secViewPtr.get());
corsika::qgsjetII::Interaction model;
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index dc237b4a3..dc5249a72 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -32,7 +32,6 @@ using namespace corsika::sibyll;
TEST_CASE("Sibyll", "[processes]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
SECTION("Sibyll -> Corsika") {
@@ -99,7 +98,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("SibyllInterface", "[processes]") {
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
logging::set_level(logging::level::trace);
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
diff --git a/tests/modules/testStackInspector.cpp b/tests/modules/testStackInspector.cpp
index 00187ed9a..0b2a3ab7b 100644
--- a/tests/modules/testStackInspector.cpp
+++ b/tests/modules/testStackInspector.cpp
@@ -24,7 +24,6 @@ using namespace corsika;
TEST_CASE("StackInspector", "[processes]") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
auto const& rootCS = get_root_CoordinateSystem();
Point const origin(rootCS, {0_m, 0_m, 0_m});
diff --git a/tests/modules/testTracking.cpp b/tests/modules/testTracking.cpp
index 87cd51e93..3f1d6ecaa 100644
--- a/tests/modules/testTracking.cpp
+++ b/tests/modules/testTracking.cpp
@@ -43,8 +43,6 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
tracking_leapfrog_straight::Tracking, tracking_line::Tracking) {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
-
logging::set_level(logging::level::trace);
const HEPEnergyType P0 = 10_GeV;
@@ -78,8 +76,9 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
int deflect = 0;
if (chargeNumber != 0 and Bfield != 0_T) {
deflect = -sgn(chargeNumber) * sgn(Bfield / 1_T); // direction of deflection
- LengthType const gyroradius =
- P0 * 1_V / (constants::c * abs(chargeNumber) * abs(Bfield) * 1_eV);
+ LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(P0) *
+ constants::c / (abs(get_charge(PID)) * abs(Bfield)));
+ CORSIKA_LOG_DEBUG("Rg={} deflect={}", gyroradius, deflect);
radius = gyroradius;
}
@@ -123,7 +122,7 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
particle.setMomentum(traj.getDirection(1) * particle.getMomentum().getNorm());
if (outer) {
// now we know we are in target volume, depending on "outer"
- CHECK(traj.getLength(1) == 0_m);
+ CHECK(traj.getLength(1) / 1_m == Approx(0).margin(1e-3));
CHECK(nextVol == targetPtr);
}
// move forward, until we leave target volume
@@ -133,6 +132,11 @@ TEMPLATE_TEST_CASE("TrackingLeapfrog_Curved", "tracking",
particle.setNode(nextVol);
particle.setPosition(traj2.getPosition(1));
particle.setMomentum(traj2.getDirection(1) * particle.getMomentum().getNorm());
+ CORSIKA_LOG_DEBUG("pos={}, p={}, |p|={} |v|={}, delta-l={}, delta-t={}",
+ particle.getPosition(), particle.getMomentum(),
+ particle.getMomentum().getNorm(),
+ particle.getVelocity().getNorm(), traj2.getLength(1),
+ traj2.getLength(1) / particle.getVelocity().getNorm());
}
CHECK(nextVol == worldPtr);
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 0f9f4a862..fd3ba64fe 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -59,7 +59,6 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
TEST_CASE("UrQMD") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
SECTION("conversion") {
CHECK_THROWS(corsika::urqmd::convertFromUrQMD(106, 0));
@@ -227,22 +226,4 @@ TEST_CASE("UrQMD") {
CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
Approx(0).margin(1e-2));
}
-
- SECTION("K0Short projectile") {
- auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
- [[maybe_unused]] auto const& env_dummy = env; // against warnings
- [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
-
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- Code::K0Short, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
-
- // must be assigned to variable, cannot be used as rvalue?!
- auto projectile = secViewPtr->getProjectile();
- auto const projectileMomentum = projectile.getMomentum();
-
- CHECK_THROWS(urqmd.doInteraction(*secViewPtr));
- }
}
--
GitLab