From 88c8ddc62f026fd19b587e8e0a2c1c4a9a5ea3f1 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Sat, 5 Jun 2021 12:30:08 +0200
Subject: [PATCH] improved leap-frog tracking
---
.../framework/geometry/LeapFrogTrajectory.inl | 19 ++-
.../detail/framework/utility/CubicSolver.inl | 153 ++++++++++--------
.../framework/utility/QuadraticSolver.inl | 2 +-
.../framework/utility/QuarticSolver.inl | 8 +-
.../tracking/TrackingLeapFrogCurved.inl | 82 ++++++----
.../tracking/TrackingLeapFrogStraight.inl | 9 +-
.../modules/tracking/TrackingStraight.inl | 16 +-
corsika/framework/geometry/Intersections.hpp | 14 +-
.../tracking/TrackingLeapFrogCurved.hpp | 21 ++-
corsika/modules/tracking/TrackingStraight.hpp | 6 +
examples/CMakeLists.txt | 4 +
tests/framework/testProcessSequence.cpp | 2 -
tests/framework/testSolver.cpp | 6 +-
13 files changed, 201 insertions(+), 141 deletions(-)
diff --git a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
index fd0e29d44..eed4221ce 100644
--- a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
+++ b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
@@ -35,19 +35,17 @@ namespace corsika {
}
inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
- return (initialDirection_ +
- initialDirection_.cross(magneticfield_) * timeStep_ * u * k_) *
- initialVelocity_.getNorm();
+ return getDirection(u) * initialVelocity_.getNorm();
}
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
- return getVelocity(u).normalized();
+ return (initialDirection_ +
+ initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
+ .normalized();
}
inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
- return u * timeStep_ *
- (1. + fabs(0.5 * initialDirection_.cross(magneticfield_).getNorm() * u *
- timeStep_ * k_));
+ return u * timeStep_;
}
inline LengthType LeapFrogTrajectory::getLength(double const u) const {
@@ -60,9 +58,10 @@ namespace corsika {
}
inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
- double const correction =
- (1. + fabs(0.5 * initialDirection_.cross(magneticfield_).getNorm() * limit * k_));
- timeStep_ = limit / correction;
+ // double const correction =
+ // (1. + fabs(0.5 * initialDirection_.cross(magneticfield_).getNorm() * limit *
+ // k_));
+ timeStep_ = limit;
}
} // namespace corsika
diff --git a/corsika/detail/framework/utility/CubicSolver.inl b/corsika/detail/framework/utility/CubicSolver.inl
index a48a0da0c..2a6dee559 100644
--- a/corsika/detail/framework/utility/CubicSolver.inl
+++ b/corsika/detail/framework/utility/CubicSolver.inl
@@ -18,6 +18,8 @@ namespace corsika {
//---------------------------------------------------------------------------
// solve cubic equation A x^3 + B*x^2 + C*x + D = 0
+ // x^3 + a*x^2 + b*x + c = 0
+ // mainly along WolframAlpha formulas
inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
long double C, long double D,
double const epsilon) {
@@ -29,8 +31,8 @@ namespace corsika {
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 q = (3 * b - a2) / 9;
+ long double r = (a * (9 * b - 2 * a2) - 27 * c) / 54;
long double q3 = q * q * q;
// disc = q**3 + r**2
@@ -38,8 +40,8 @@ namespace corsika {
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);
+ 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;
@@ -51,30 +53,34 @@ namespace corsika {
au = f - au;
ab = u - ab;
// sum all terms into final result
- long double const disc = -(((e + uf) + au) + ab) + v;
+ long double const disc = (((e + uf) + au) + ab) + v;
- if (disc >= 0) {
- long double t = r / std::sqrt(q3);
+ CORSIKA_LOG_TRACE("disc={} {}", disc, q3 + r * r);
+
+ if (std::abs(disc) < epsilon) {
+
+ a /= 3;
+ long double const cbrtR = std::cbrt(r);
+ return {double(2 * cbrtR - a), double(-cbrtR - a)}; // 2nd solution is doublet
+
+ } else if (disc > 0) {
+
+ long double const S = std::cbrt(r + std::sqrt(disc));
+ long double const T = std::cbrt(r - std::sqrt(disc));
+ a /= 3;
+ return {double((S + T) - a)}; // plus two imaginary solution
+
+ } else { // 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);
+ 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 term1 = -cbrt(std::fabs(r) + std::sqrt(-disc));
- if (r < 0) term1 = -term1;
- long double term2 = (0 == term1 ? 0 : q / term1);
-
- a /= 3;
- long double test = 0.5 * std::sqrt(3.) * (term1 - term2);
- if (std::fabs(test) < epsilon) {
- return {double((term1 + term2) - 1), double(-0.5 * (term1 + term2) - a)};
- }
-
- return {double((term1 + term2) - a)};
+ double(q * std::cos((t + 4 * M_PI) / 3) - a)};
}
}
} // namespace andre
@@ -228,7 +234,8 @@ namespace corsika {
}
/**
- * Iterative approach.
+ * Iterative approach. https://en.wikipedia.org/wiki/Halley%27s_method
+ * Halley's method
*/
inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
@@ -238,72 +245,80 @@ namespace corsika {
a, b, c, d, epsilon, (std::abs(a - 1) < epsilon),
(std::abs(b) < epsilon));
-#ifdef DEBUG
- {
- auto test = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
-
- for (long double test_v : test) {
- CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
- cubic_function(test_v, a, b, c, d));
- }
- }
-#endif
-
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);
+ auto pre_opt = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+ long double x1 = 0; // start value
- long double x1 = xinfl;
- long double f_x1 = cubic_function(xinfl, 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);
+ if (pre_opt.size()) {
+ x1 = pre_opt[0]; //*std::max_element(pre_opt.begin(), pre_opt.end());
+#ifdef DEBUG
+ for (long double test_v : pre_opt) {
+ CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
+ cubic_function(test_v, a, b, c, d));
}
-
- 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 / (static_pow<2>(f_prime_x1) - 0.5 * f_x1 * f_prime2_x1);
+#endif
+ } else {
+ long double const dist = std::fma(b / a, b / a, -3 * c / a);
+ long double const xinfl = -b / (a * 3);
+
+ x1 = xinfl;
+ long double f_test = cubic_function(xinfl, a, b, c, d);
+
+ if (std::abs(f_test) > epsilon) {
+ if (std::abs(dist) < epsilon) {
+ x1 = xinfl - std::cbrt(f_test);
+ } else if (dist > 0) {
+ if (f_test > 0)
+ x1 = xinfl - 2 / 3 * std::sqrt(dist);
+ else
+ x1 = xinfl + 2 / 3 * std::sqrt(dist);
}
- 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);
- if (niter >= maxiter) {
- // CORSIKA_LOG_TRACE("failure, no solution");
- // return std::vector<double>{};
}
}
+ long double f_x1 = cubic_function(x1, a, b, c, d);
+ long double dx1 = 0;
+
+ 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) {
+ dx1 = std::cbrt(f_x1);
+ } else {
+ dx1 = f_x1 * f_prime_x1 * 2 / (f_prime_x1 * f_prime_x1 * 2 - f_x1 * f_prime2_x1);
+ }
+ x1 -= dx1;
+ f_x1 = cubic_function(x1, a, b, c, d);
+ CORSIKA_LOG_TRACE(
+ "niter={} x1={:.20f} f_x1={:.20f} f_prime={:.20f} f_prime2={:.20f} dx1={}",
+ niter, x1, f_x1, f_prime_x1, f_prime2_x1,
+ f_x1 * f_prime_x1 / (f_prime_x1 * f_prime_x1 - f_x1 * f_prime2_x1 * 0.5));
+ } while ((++niter < maxiter) && (std::abs(f_x1) > epsilon * 1000) &&
+ (std::abs(dx1) > epsilon));
+
+ CORSIKA_LOG_TRACE("niter={}", niter);
+ if (niter >= maxiter) {
+ CORSIKA_LOG_DEBUG("niter reached max iterations {}", niter);
+ return andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
+ }
+
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);
+ std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, 1e-3);
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
} // namespace corsika
diff --git a/corsika/detail/framework/utility/QuadraticSolver.inl b/corsika/detail/framework/utility/QuadraticSolver.inl
index 9c9759e97..a335b00d6 100644
--- a/corsika/detail/framework/utility/QuadraticSolver.inl
+++ b/corsika/detail/framework/utility/QuadraticSolver.inl
@@ -62,7 +62,7 @@ namespace corsika {
long double f = std::fma(b, b, -w);
long double radicant = f + e;
- CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4);
+ CORSIKA_LOG_TRACE("radicant={} {} ", radicant, b * b - a * c * 4);
if (std::abs(radicant) < epsilon) { // just one real solution
return {double(-b / (2 * a))};
diff --git a/corsika/detail/framework/utility/QuarticSolver.inl b/corsika/detail/framework/utility/QuarticSolver.inl
index ebdfd932e..84a9c5165 100644
--- a/corsika/detail/framework/utility/QuarticSolver.inl
+++ b/corsika/detail/framework/utility/QuarticSolver.inl
@@ -77,8 +77,8 @@ namespace corsika {
// solving quadratic eqs. x^2 + p1*x + q1 = 0
// x^2 + p2*x + q2 = 0
- std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
- std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
+ std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
+ std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
if (quad2.size() > 0) {
for (auto val : quad2) quad1.push_back(val);
}
@@ -140,9 +140,9 @@ namespace corsika {
long double const quad_term3 = q / (2 * quad_term2);
std::vector<double> z_quad1 =
- solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, epsilon);
+ solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, 1e-5);
std::vector<double> z_quad2 =
- solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, epsilon);
+ solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, 1e-5);
for (auto const& z : z_quad2) z_quad1.push_back(z);
return z_quad1;
}
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
index 4542f2e56..a07e74963 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -29,10 +29,10 @@ namespace corsika {
namespace tracking_leapfrog_curved {
template <typename TParticle>
- inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) {
+ inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
if (particle.getMomentum().getNorm() == 0_GeV) {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
- double(0));
+ double(0) * 1_m);
} // charge of the particle
ElectricChargeType const charge = particle.getCharge();
auto const* currentLogicalVolumeNode = particle.getNode();
@@ -40,9 +40,14 @@ namespace corsika {
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
VelocityVector velocity = particle.getVelocity();
- decltype(meter / (second * volt)) k =
- charge * constants::cSquared * 1_eV /
- (velocity.getNorm() * 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 * meter)) const k =
+ charge / p_norm; // * initialVelocity.getNorm();
+
DirectionVector direction = velocity.normalized();
auto position = particle.getPosition(); // First Movement
// assuming magnetic field does not change during movement
@@ -133,22 +138,18 @@ namespace corsika {
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)
- // correction factor to move from leap-frog length L to half-step deltaL/2
- double const correction =
- 1. + fabs(0.5 * particle.getDirection().cross(magneticfield).getNorm() * k *
- minTime);
-
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
- minTime / correction), // --> trajectory
- minNode); // --> next volume node
+ minTime), // --> trajectory
+ minNode); // --> next volume node
}
template <typename TParticle>
inline Intersections Tracking::intersect(TParticle const& particle,
Sphere const& sphere) {
- if (sphere.getRadius() == 1_km * std::numeric_limits<double>::infinity()) {
+ LengthType const radius = sphere.getRadius();
+ if (radius == 1_km * std::numeric_limits<double>::infinity()) {
return Intersections();
}
@@ -174,6 +175,30 @@ namespace corsika {
bool const numericallyInside = sphere.contains(particle.getPosition());
CORSIKA_LOG_TRACE("numericallyInside={}", numericallyInside);
+ Vector<length_d> const deltaPos = position - sphere.getCenter();
+
+ { // check extreme cases we don't want to solve analytically explicit
+ HEPMomentumType const p_perp =
+ (particle.getMomentum() -
+ particle.getMomentum().getParallelProjectionOnto(magneticfield))
+ .getNorm();
+
+ LengthType const gyroradius =
+ (convert_HEP_to_SI<MassType::dimension_type>(p_perp) * constants::c /
+ (abs(charge) * magneticfield.getNorm()));
+
+ LengthType const trackDist = abs(deltaPos.getNorm() - radius);
+ if (trackDist > gyroradius) {
+ // there is never a solution
+ return Intersections();
+ }
+
+ if (gyroradius > 100 * trackDist) {
+ // the bending is negligible, use straight intersections instead
+ return tracking_line::Tracking::intersect(particle, sphere);
+ }
+ }
+
SpeedType const absVelocity = velocity.getNorm();
auto const p_norm =
constants::c * convert_HEP_to_SI<MassType::dimension_type>(
@@ -183,12 +208,10 @@ namespace corsika {
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_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 /
+ double const d = (deltaPosLength + radius) * (deltaPosLength - radius) * denom /
(1_m * 1_m * 1_m * 1_m);
CORSIKA_LOG_TRACE("denom={}, b={}, c={}, d={}", denom, b, c, d);
// solutions of deltaL are obtained from quartic equation. Note, deltaL/2 is the
@@ -196,18 +219,14 @@ namespace corsika {
// because of the non-conservation of norm/velocity.
// The leap-frog length L is deltaL/2 * (1+|u_{n+1}|)
std::vector<double> solutions = solve_quartic_real(1, 0, b, c, d);
- if (!solutions.size()) {
- return tracking_line::Tracking::intersect<TParticle>(particle, sphere);
- }
+ if (!solutions.size()) { return Intersections(); }
LengthType d_enter, d_exit;
int first = 0, first_entry = 0, first_exit = 0;
for (auto solution : solutions) {
LengthType const dist = solution * 1_m;
CORSIKA_LOG_TRACE(
- "Solution (real) for current Volume: deltaL/2*2={} (deltaL/2*2/v={}, "
- "corrected={}) ",
- dist, dist / absVelocity,
- dist / absVelocity * (1. + fabs(0.5 * direction_x_B.getNorm() * k * dist)));
+ "Solution (real) for current Volume: deltaL/2*2={} (deltaL/2*2/v={}) ", dist,
+ dist / absVelocity);
if (numericallyInside) {
// there must be an entry (negative) and exit (positive) solution
if (dist < 0.0001_m) { // security margin to assure
@@ -265,12 +284,9 @@ namespace corsika {
first_entry, first_exit);
return Intersections();
}
- // return in units of time, correct to leap-frog length L
- double const corr_enter = 1. + fabs(0.5 * direction_x_B.getNorm() * k * d_enter);
- double const corr_exit = 1. + fabs(0.5 * direction_x_B.getNorm() * k * d_exit);
+ // return in units of time
- return Intersections(d_enter * corr_enter / absVelocity,
- d_exit * corr_exit / absVelocity);
+ return Intersections(d_enter / absVelocity, d_exit / absVelocity);
}
template <typename TParticle>
@@ -345,11 +361,13 @@ namespace corsika {
CORSIKA_LOG_TRACE("maxStepLength={} s", maxStepLength / 1_s);
// with final length correction, |direction| becomes >1 during step
- double const corr =
- 1. + fabs(0.5 * direction_x_B.getNorm() * maxStepLength * charge / p_norm);
+ // double const correction =
+ // 1. + fabs(0.5 * direction_x_B.getNorm() * maxStepLength * charge /
+ // p_norm);
+
+ return Intersections(maxStepLength / absVelocity); // unit: s
- return Intersections(maxStepLength * corr / absVelocity); // unit: s
- } // no charge
+ } // no charge
CORSIKA_LOG_TRACE("(plane) straight tracking with charge={}, B={}", charge,
particle.getNode()->getModelProperties().getMagneticField(
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
index 4e7efed77..675636a17 100644
--- a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -63,7 +63,7 @@ namespace corsika {
// check, where the first halve-step direction has geometric intersections
auto const [initialTrack, initialTrackNextVolume] =
tracking_line::Tracking::getTrack(particle);
- { [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
+ //{ [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
auto const initialTrackLength = initialTrack.getLength(1);
CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
@@ -82,8 +82,8 @@ namespace corsika {
.getNorm();
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.
+ // particle travels along, parallel to magnetic field. Rg is
+ // "0", return straight track here.
CORSIKA_LOG_TRACE("p_perp is 0_GeV --> parallel");
return std::make_tuple(initialTrack, initialTrackNextVolume);
}
@@ -105,8 +105,9 @@ namespace corsika {
DirectionVector const direction = initialVelocity.normalized();
// avoid any intersections within first halve steplength
+ // aim for intersection in second halve step
LengthType const firstHalveSteplength =
- std::min(steplimit, initialTrackLength * firstFraction_);
+ std::min(steplimit / 2, initialTrackLength * firstFraction_);
CORSIKA_LOG_DEBUG("first halve step length {}, steplimit={}, initialTrackLength={}",
firstHalveSteplength, steplimit, initialTrackLength);
diff --git a/corsika/detail/modules/tracking/TrackingStraight.inl b/corsika/detail/modules/tracking/TrackingStraight.inl
index cad8ac39d..3bdbb73df 100644
--- a/corsika/detail/modules/tracking/TrackingStraight.inl
+++ b/corsika/detail/modules/tracking/TrackingStraight.inl
@@ -24,6 +24,15 @@
namespace corsika::tracking_line {
+ template <typename TParticle>
+ inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
+ if (particle.getMomentum().getNorm() == 0_GeV) {
+ return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV);
+ } // charge of the particle
+ DirectionVector const dir = particle.getDirection();
+ return std::make_tuple(particle.getPosition() + dir * steplength, dir.normalized());
+ }
+
template <typename TParticle>
inline auto Tracking::getTrack(TParticle const& particle) {
VelocityVector const initialVelocity =
@@ -92,9 +101,10 @@ namespace corsika::tracking_line {
CORSIKA_LOG_TRACE("n_dot_v={}, delta={}, momentum={}", n_dot_v, delta,
particle.getMomentum());
- return Intersections(n_dot_v.magnitude() == 0
- ? std::numeric_limits<TimeType::value_type>::infinity() * 1_s
- : n.dot(delta) / n_dot_v);
+ if (n_dot_v.magnitude() == 0)
+ return Intersections();
+ else
+ return Intersections(n.dot(delta) / n_dot_v);
}
} // namespace corsika::tracking_line
diff --git a/corsika/framework/geometry/Intersections.hpp b/corsika/framework/geometry/Intersections.hpp
index c9951213e..ce0131ee3 100644
--- a/corsika/framework/geometry/Intersections.hpp
+++ b/corsika/framework/geometry/Intersections.hpp
@@ -43,10 +43,20 @@ namespace corsika {
bool hasIntersections() const { return has_intersections_; }
///! where did the trajectory currently enter the volume
- TimeType getEntry() const { return intersections_.first; }
+ TimeType getEntry() const {
+ if (has_intersections_)
+ return intersections_.first;
+ else
+ return std::numeric_limits<TimeType::value_type>::infinity() * second;
+ }
///! where did the trajectory currently exit the volume
- TimeType getExit() const { return intersections_.second; }
+ TimeType getExit() const {
+ if (has_intersections_)
+ return intersections_.second;
+ else
+ return std::numeric_limits<TimeType::value_type>::infinity() * second;
+ }
private:
bool has_intersections_;
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
index c4ed5aafd..66bc8136f 100644
--- a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -29,24 +29,16 @@ namespace corsika {
namespace tracking_leapfrog_curved {
/**
- * \file TrackingLeapFrogCurved.hpp
- *
- * Performs one leap-frog step consistent of two halve-steps with steplength/2
- * The step is caluculated analytically precisely to reach to the next volume
- * boundary.
+ * \file TrackingLeapFrogCurved.hpp The leap-frog tracking.
*/
- template <typename TParticle>
- auto make_LeapFrogStep(TParticle const& particle, LengthType steplength);
/**
- *
* The class tracking_leapfrog_curved::Tracking is based on the
* Bachelor thesis of Andre Schmidt (KIT). It implements a
* two-step leap-frog algorithm, but with analytically exact geometric
* intersections between leap-frog steps and geometric volumes
* (spheres, planes).
- *
- **/
+ */
class Tracking : public Intersect<Tracking> {
@@ -59,6 +51,14 @@ namespace corsika {
template <typename TParticle>
auto getTrack(TParticle const& particle);
+ /**
+ * Performs one leap-frog step consistent of two halve-steps with steplength/2
+ * Due to the nature of the algorithm the second halve step is slightly longer than
+ * the first halve step.
+ */
+ template <typename TParticle>
+ static auto makeStep(TParticle const& particle, LengthType const steplength);
+
/**
* find intersection of Sphere with Track
*
@@ -77,7 +77,6 @@ namespace corsika {
*
* The intersection time(s) of a particle, assuming a curved leap-frog
* step, are calculated for any volume type.
- *
*/
template <typename TParticle, typename TBaseNodeType>
static Intersections intersect(TParticle const& particle,
diff --git a/corsika/modules/tracking/TrackingStraight.hpp b/corsika/modules/tracking/TrackingStraight.hpp
index a5800cd7e..3ebc18b37 100644
--- a/corsika/modules/tracking/TrackingStraight.hpp
+++ b/corsika/modules/tracking/TrackingStraight.hpp
@@ -36,6 +36,12 @@ namespace corsika::tracking_line {
using Intersect<Tracking>::nextIntersect;
public:
+ /**
+ * Performs one straight step of length steplength.
+ */
+ template <typename TParticle>
+ static auto makeStep(TParticle const& particle, LengthType const steplength);
+
//! Determine track of particle
template <typename TParticle>
auto getTrack(TParticle const& particle);
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 7aa0ad850..9cb9b3fd3 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -59,3 +59,7 @@ CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
add_executable (corsika corsika.cpp)
target_link_libraries (corsika CORSIKA8::CORSIKA8)
+
+add_executable (tracking tracking.cpp)
+target_link_libraries (tracking CORSIKA8::CORSIKA8)
+
diff --git a/tests/framework/testProcessSequence.cpp b/tests/framework/testProcessSequence.cpp
index 1335974cd..cf9dc523f 100644
--- a/tests/framework/testProcessSequence.cpp
+++ b/tests/framework/testProcessSequence.cpp
@@ -369,7 +369,6 @@ private:
TEST_CASE("ProcessSequence General", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
SECTION("BaseProcess") {
@@ -889,7 +888,6 @@ TEST_CASE("SwitchProcessSequence", "ProcessSequence") {
TEST_CASE("ProcessSequence Indexing", "ProcessSequence") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$]: %v");
SECTION("Indexing") {
diff --git a/tests/framework/testSolver.cpp b/tests/framework/testSolver.cpp
index 089335ed5..aa8885bb0 100644
--- a/tests/framework/testSolver.cpp
+++ b/tests/framework/testSolver.cpp
@@ -168,7 +168,7 @@ TEST_CASE("Solver") {
z1, z2, z3, a, b, c, d);
vector<double> s1 = solve_cubic_real(a, b, c, d);
- remove_duplicates(s1, epsilon_check);
+ remove_duplicates(s1, epsilon_check * 10);
CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
@@ -193,7 +193,7 @@ TEST_CASE("Solver") {
SECTION("quartic") {
- epsilon_check = 1e-4; // for catch2 asserts
+ epsilon_check = 1e-2; // for catch2 asserts
// **clang-format-off**
// tests of type:
@@ -252,7 +252,7 @@ TEST_CASE("Solver") {
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);
+ remove_duplicates(s1, epsilon_check * 10);
CORSIKA_LOG_INFO("N={}, s1=[{}]", s1.size(), fmt::join(s1, ", "));
--
GitLab