From 8dc6d688f9b86e2d94452d8e8c19e30ad11e1491 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Thu, 31 Dec 2020 12:00:08 +0100
Subject: [PATCH] missing files
---
corsika/detail/modules/TrackingLine.inl | 59 ++++
corsika/modules/Tracking.hpp | 18 ++
.../tracking/TrackingLeapFrogCurved.hpp | 298 ++++++++++++++++++
.../tracking/TrackingLeapFrogStraight.hpp | 225 +++++++++++++
4 files changed, 600 insertions(+)
create mode 100644 corsika/detail/modules/TrackingLine.inl
create mode 100644 corsika/modules/Tracking.hpp
create mode 100644 corsika/modules/tracking/TrackingLeapFrogCurved.hpp
create mode 100644 corsika/modules/tracking/TrackingLeapFrogStraight.hpp
diff --git a/corsika/detail/modules/TrackingLine.inl b/corsika/detail/modules/TrackingLine.inl
new file mode 100644
index 000000000..52200c42d
--- /dev/null
+++ b/corsika/detail/modules/TrackingLine.inl
@@ -0,0 +1,59 @@
+/*
+ * (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/geometry/Point.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/media/Environment.hpp>
+
+#include <limits>
+#include <stdexcept>
+#include <utility>
+
+namespace corsika::tracking_line {
+
+ std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(
+ corsika::Line const& line, corsika::Sphere const& sphere) {
+ auto const delta = line.getStartPoint() - sphere.getCenter();
+ auto const v = line.getVelocity();
+ auto const vSqNorm =
+ v.getSquaredNorm(); // todo: get rid of this by having V0 normalized always
+ auto const R = sphere.getRadius();
+
+ auto const vDotDelta = v.dot(delta);
+ auto const discriminant =
+ vDotDelta * vDotDelta - vSqNorm * (delta.getSquaredNorm() - R * R);
+
+ if (discriminant.magnitude() > 0) {
+ auto const sqDisc = sqrt(discriminant);
+ auto const invDenom = 1 / vSqNorm;
+ return std::make_pair((-vDotDelta - sqDisc) * invDenom,
+ (-vDotDelta + sqDisc) * invDenom);
+ } else {
+ return {};
+ }
+ }
+
+ TimeType getTimeOfIntersection(Line const& vLine, Plane const& vPlane) {
+
+ auto const delta = vPlane.getCenter() - vLine.getStartPoint();
+ auto const v = vLine.getVelocity();
+ auto const n = vPlane.getNormal();
+ auto const c = n.dot(v);
+
+ if (c.magnitude() == 0) {
+ return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
+ } else {
+ return n.dot(delta) / c;
+ }
+ }
+
+} // namespace corsika::tracking_line
diff --git a/corsika/modules/Tracking.hpp b/corsika/modules/Tracking.hpp
new file mode 100644
index 000000000..fb12412fb
--- /dev/null
+++ b/corsika/modules/Tracking.hpp
@@ -0,0 +1,18 @@
+/*
+ * (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/modules/tracking/TrackingStraight.hpp>
+#include <corsika/modules/tracking/TrackingLeapFrogStraight.hpp> // simple leap-frog implementation with two straight lines
+#include <corsika/modules/tracking/TrackingLeapFrogCurved.hpp> // // more complete, curved, leap-frog
+
+/**
+ * \todo add TrackingCurved, with simple circular trajectories
+ \todo add Boris algorithm
+ */
diff --git a/corsika/modules/tracking/TrackingLeapFrogCurved.hpp b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
new file mode 100644
index 000000000..bb75daf62
--- /dev/null
+++ b/corsika/modules/tracking/TrackingLeapFrogCurved.hpp
@@ -0,0 +1,298 @@
+/*
+ * (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/modules/tracking/TrackingStraight.hpp> // for neutral particles
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
+#include <corsika/framework/logging/Logging.hpp>
+#include <corsika/modules/tracking/Intersect.hpp>
+
+#include <type_traits>
+#include <utility>
+
+namespace corsika {
+
+ namespace tracking_leapfrog_curved {
+
+ /**
+ * \function LeapFrogStep
+ *
+ * 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.
+ **/
+ template <typename TParticle>
+ auto LeapFrogStep(TParticle const& particle, LengthType steplength) {
+ if (particle.getMomentum().norm() == 0_GeV) {
+ return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
+ double(0));
+ } // charge of the particle
+ int const chargeNumber = particle.getChargeNumber();
+ auto const* currentLogicalVolumeNode = particle.getNode();
+ MagneticFieldVector const& magneticfield =
+ currentLogicalVolumeNode->getModelProperties().getMagneticField(
+ particle.getPosition());
+ VelocityVector velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ decltype(meter / (second * volt)) k =
+ chargeNumber * constants::cSquared * 1_eV /
+ (velocity.getNorm() * particle.getEnergy() * 1_V);
+ DirectionVector direction = velocity.normalized();
+ auto position = particle.getPosition(); // First Movement
+ // assuming magnetic field does not change during movement
+ position =
+ position + direction * steplength / 2; // Change of direction by magnetic field
+ direction =
+ direction + direction.cross(magneticfield) * steplength * k; // Second Movement
+ position = position + direction * steplength / 2;
+ auto steplength_true = steplength * (1.0 + (double)direction.getNorm()) / 2;
+ return std::make_tuple(position, direction.normalized(), steplength_true);
+ }
+
+ /**
+ *
+ * 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> {
+
+ using Intersect<Tracking>::nextIntersect;
+
+ public:
+ Tracking()
+ : straightTracking_{tracking_line::Tracking()} {}
+
+ template <typename TParticle>
+ auto getTrack(TParticle const& particle) {
+ VelocityVector const initialVelocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+
+ auto const position = particle.getPosition();
+ CORSIKA_LOG_DEBUG(
+ "Tracking pid: {}"
+ " , E = {} GeV",
+ particle.getPID(), particle.getEnergy() / 1_GeV);
+ CORSIKA_LOG_DEBUG("Tracking pos: {}", position.getCoordinates());
+ CORSIKA_LOG_DEBUG("Tracking E: {} GeV", particle.getEnergy() / 1_GeV);
+ CORSIKA_LOG_DEBUG("Tracking p: {} GeV",
+ particle.getMomentum().getComponents() / 1_GeV);
+ CORSIKA_LOG_DEBUG("Tracking v: {} ", initialVelocity.getComponents());
+
+ typedef
+ typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
+ node_type& 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
+ // trajectories segment-wise -- at least if we don't employ concepts as "Helix
+ // Trajectories" or similar
+ 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;
+
+ if (no_deflection) { return getLinearTrajectory(particle); }
+
+ HEPMomentumType const pAlongB_delta =
+ (particle.getMomentum() -
+ particle.getMomentum().getParallelProjectionOnto(magneticfield))
+ .getNorm();
+
+ if (pAlongB_delta == 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");
+ return getLinearTrajectory(particle);
+ }
+
+ LengthType const gyroradius =
+ (pAlongB_delta * 1_V /
+ (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+
+ const double maxRadians = 0.01;
+ const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
+ const TimeType steplimit_time = steplimit / initialVelocity.getNorm();
+ CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
+ steplimit_time);
+
+ // traverse the environment volume tree and find next
+ // intersection
+ auto [minTime, minNode] = nextIntersect(particle, steplimit_time);
+
+ const auto k =
+ chargeNumber * constants::cSquared * 1_eV / (particle.getEnergy() * 1_V);
+ return std::make_tuple(
+ LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
+ minTime), // trajectory
+ minNode); // next volume node
+ }
+
+ template <typename TParticle, typename TMedium>
+ static Intersections intersect(const TParticle& particle, const Sphere& sphere,
+ const TMedium& medium) {
+
+ if (sphere.getRadius() == 1_km * std::numeric_limits<double>::infinity()) {
+ return Intersections();
+ }
+
+ int const chargeNumber = particle.getChargeNumber();
+ auto const& position = particle.getPosition();
+ MagneticFieldVector const& magneticfield = medium.getMagneticField(position);
+
+ VelocityVector const velocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+ DirectionVector const directionBefore =
+ velocity.normalized(); // determine steplength to next volume
+
+ 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) {
+ return tracking_line::Tracking::intersect(particle, sphere, medium);
+ }
+
+ bool const numericallyInside = sphere.isInside(particle.getPosition());
+
+ const auto absVelocity = velocity.getNorm();
+ auto energy = particle.getEnergy();
+ auto k = chargeNumber * constants::cSquared * 1_eV / (absVelocity * energy * 1_V);
+
+ auto const denom =
+ (directionBefore.cross(magneticfield)).getSquaredNorm() * k * k;
+ const double a =
+ ((directionBefore.cross(magneticfield)).dot(position - sphere.getCenter()) *
+ k +
+ 1) *
+ 4 / (1_m * 1_m * denom);
+ const double b = directionBefore.dot(position - sphere.getCenter()) * 8 /
+ (denom * 1_m * 1_m * 1_m);
+ const double c = ((position - sphere.getCenter()).getSquaredNorm() -
+ (sphere.getRadius() * sphere.getRadius())) *
+ 4 / (denom * 1_m * 1_m * 1_m * 1_m);
+ CORSIKA_LOG_TRACE("denom={}, a={}, b={}, c={}", denom, a, b, c);
+ std::complex<double>* solutions = quartic_solver::solve_quartic(0, a, b, c);
+ 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++;
+ }
+ 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
+ // into next logical volume
+ continue;
+ }
+ if (first == 0) {
+ d_enter = dist;
+ } else {
+ if (dist < d_enter) {
+ d_exit = d_enter;
+ d_enter = dist;
+ } else {
+ d_exit = dist;
+ }
+ }
+ first++;
+ }
+ } // loop over solutions
+ }
+ delete[] solutions;
+
+ if (first != 2) { // entry and exit points found
+ CORSIKA_LOG_DEBUG("no intersection! count={}", first);
+ return Intersections();
+ }
+ return Intersections(d_enter / absVelocity, d_exit / absVelocity);
+ }
+
+ template <typename TParticle, typename TBaseNodeType>
+ static Intersections intersect(const TParticle& particle,
+ const TBaseNodeType& volumeNode) {
+ const Sphere* sphere = dynamic_cast<const Sphere*>(&volumeNode.getVolume());
+ if (sphere) {
+ return intersect(particle, *sphere, volumeNode.getModelProperties());
+ }
+ throw std::runtime_error(
+ "The Volume type provided is not supported in intersect(particle, node)");
+ }
+
+ protected:
+ /**
+ * Use internally stored class tracking_line::Tracking to
+ * perform a straight line tracking, if no magnetic bendig was
+ * detected.
+ *
+ */
+ template <typename TParticle>
+ auto getLinearTrajectory(TParticle& particle) {
+
+ // perform simple linear tracking
+ auto [straightTrajectory, minNode] = straightTracking_.getTrack(particle);
+
+ // return as leap-frog trajectory
+ return std::make_tuple(
+ LeapFrogTrajectory(
+ straightTrajectory.getLine().getStartPoint(),
+ straightTrajectory.getLine().getVelocity(),
+ MagneticFieldVector(particle.getPosition().getCoordinateSystem(), 0_T,
+ 0_T, 0_T),
+ square(0_m) / (square(1_s) * 1_V),
+ straightTrajectory.getDuration()), // trajectory
+ minNode); // next volume node
+ }
+
+ protected:
+ tracking_line::Tracking
+ straightTracking_; ///! we want this for neutral and B=0T tracks
+
+ }; // namespace tracking_leapfrog_curved
+
+ } // namespace tracking_leapfrog_curved
+
+} // namespace corsika
diff --git a/corsika/modules/tracking/TrackingLeapFrogStraight.hpp b/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
new file mode 100644
index 000000000..0bc38b93d
--- /dev/null
+++ b/corsika/modules/tracking/TrackingLeapFrogStraight.hpp
@@ -0,0 +1,225 @@
+/*
+ * (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/geometry/Line.hpp>
+#include <corsika/framework/geometry/Plane.hpp>
+#include <corsika/framework/geometry/Sphere.hpp>
+#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <corsika/modules/tracking/TrackingStraight.hpp>
+
+#include <type_traits>
+#include <utility>
+
+namespace corsika {
+
+ namespace tracking_leapfrog_straight {
+
+ /**
+ *
+ * The class tracking_leapfrog_straight::Tracking inherits from
+ * tracking_line::Tracking and adds a (two-step) Leap-Frog
+ * algorithms with two halve-steps and magnetic deflection.
+ *
+ * The two halve steps are implemented as two straight explicit
+ * `tracking_line::Tracking`s and all geometry intersections are,
+ * thus, based on those two straight line elements.
+ *
+ * As a precaution for numerical instability, the steplength is
+ * limited to correspond to a straight line distance to the next
+ * volume intersection. In typical situations this leads to about
+ * (at least) one full leap-frog step to the next volume boundary.
+ *
+ **/
+
+ class Tracking : public tracking_line::Tracking {
+
+ public:
+ /**
+ * \param firstFraction fraction of first leap-frog halve step
+ * relative to full linear step to next volume boundary. This
+ * should not be less than 0.5, otherwise you risk that
+ * particles will never travel from one volume to the next
+ * one. A cross should be possible (even likely). If
+ * firstFraction is too big (~1) the resulting calculated error
+ * will be largest.
+ *
+ */
+ Tracking(double const firstFraction = 0.55)
+ : firstFraction_(firstFraction) {}
+
+ template <typename Particle>
+ auto getTrack(Particle& particle) {
+ VelocityVector initialVelocity =
+ particle.getMomentum() / particle.getEnergy() * constants::c;
+
+ const Point initialPosition = particle.getPosition();
+ CORSIKA_LOG_DEBUG(
+ "TrackingB pid: {}"
+ " , E = {} GeV",
+ particle.getPID(), particle.getEnergy() / 1_GeV);
+ CORSIKA_LOG_DEBUG("TrackingB pos: {}", initialPosition.getCoordinates());
+ CORSIKA_LOG_DEBUG("TrackingB E: {} GeV", particle.getEnergy() / 1_GeV);
+ CORSIKA_LOG_DEBUG("TrackingB p: {} GeV",
+ particle.getMomentum().getComponents() / 1_GeV);
+ CORSIKA_LOG_DEBUG("TrackingB v: {} ", initialVelocity.getComponents());
+
+ typedef decltype(particle.getNode()) node_type;
+ node_type const volumeNode = particle.getNode();
+
+ // check if particle is moving at all
+ auto const absVelocity = initialVelocity.getNorm();
+ if (absVelocity * 1_s == 0_m) {
+ return std::make_tuple(
+ LineTrajectory(Line(initialPosition, initialVelocity), 0_s), volumeNode);
+ }
+
+ // charge of the particle, and magnetic field
+ const int chargeNumber = particle.getChargeNumber();
+ auto magneticfield =
+ volumeNode->getModelProperties().getMagneticField(initialPosition);
+ const auto 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;
+
+ // check, where the first halve-step direction has geometric intersections
+ const auto [initialTrack, initialTrackNextVolume] =
+ tracking_line::Tracking::getTrack(particle);
+ { [[maybe_unused]] auto& initialTrackNextVolume_dum = initialTrackNextVolume; }
+ const auto initialTrackLength = initialTrack.getLength(1);
+
+ CORSIKA_LOG_DEBUG("initialTrack(0)={}, initialTrack(1)={}, initialTrackLength={}",
+ initialTrack.getPosition(0).getCoordinates(),
+ initialTrack.getPosition(1).getCoordinates(),
+ initialTrackLength);
+
+ // if particle is non-deflectable, we are done:
+ if (no_deflection) {
+ CORSIKA_LOG_DEBUG("no deflection. tracking finished");
+ return std::make_tuple(initialTrack, initialTrackNextVolume);
+ }
+
+ HEPMomentumType const pAlongB_delta =
+ (particle.getMomentum() -
+ particle.getMomentum().getParallelProjectionOnto(magneticfield))
+ .getNorm();
+
+ if (pAlongB_delta == 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");
+ return std::make_tuple(initialTrack, initialTrackNextVolume);
+ }
+
+ LengthType const gyroradius =
+ (pAlongB_delta * 1_V /
+ (constants::c * abs(chargeNumber) * magnitudeB * 1_eV));
+
+ // we need to limit maximum step-length since we absolutely
+ // need to follow strongly curved trajectories segment-wise,
+ // at least if we don't employ concepts as "Helix
+ // Trajectories" or similar
+ const double maxRadians = 0.01;
+ const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
+ CORSIKA_LOG_DEBUG("gyroradius {}, Steplimit: {}", gyroradius, steplimit);
+
+ // calculate first halve step for "steplimit"
+ auto const initialMomentum = particle.getMomentum();
+ auto const absMomentum = initialMomentum.getNorm();
+ DirectionVector const direction = initialVelocity.normalized();
+
+ // avoid any intersections within first halve steplength
+ LengthType const firstHalveSteplength =
+ std::min(steplimit, initialTrackLength * firstFraction_);
+
+ CORSIKA_LOG_DEBUG(
+ "first halve step length {}, steplimit={}, initialTrackLength={}",
+ firstHalveSteplength, steplimit, initialTrackLength);
+ // perform the first halve-step
+ const Point position_mid = initialPosition + direction * firstHalveSteplength;
+ const auto k =
+ chargeNumber * constants::c * 1_eV / (particle.getMomentum().getNorm() * 1_V);
+ const auto new_direction =
+ direction + direction.cross(magneticfield) * firstHalveSteplength * 2 * k;
+ const auto new_direction_norm = new_direction.getNorm(); // by design this is >1
+ CORSIKA_LOG_DEBUG(
+ "position_mid={}, new_direction={}, (new_direction_norm)={}, deflection={}",
+ position_mid.getCoordinates(), new_direction.getComponents(),
+ new_direction_norm,
+ acos(std::min(1.0, direction.dot(new_direction) / new_direction_norm)) * 180 /
+ M_PI);
+
+ // check, where the second halve-step direction has geometric intersections
+ particle.setPosition(position_mid);
+ particle.setMomentum(new_direction * absMomentum);
+ const auto [finalTrack, finalTrackNextVolume] =
+ tracking_line::Tracking::getTrack(particle);
+ particle.setPosition(initialPosition); // this is not nice...
+ particle.setMomentum(initialMomentum); // this is not nice...
+
+ LengthType const finalTrackLength = finalTrack.getLength(1);
+ LengthType const secondLeapFrogLength = firstHalveSteplength * new_direction_norm;
+
+ // check if volume transition is obvious, OR
+ // for numerical reasons, particles slighly bend "away" from a
+ // volume boundary have a very hard time to cross the border,
+ // thus, if secondLeapFrogLength is just slighly shorter (1e-4m) than
+ // finalTrackLength we better just [extend the
+ // secondLeapFrogLength slightly and] force the volume
+ // crossing:
+ bool const switch_volume = finalTrackLength - 0.0001_m <= secondLeapFrogLength;
+ LengthType const secondHalveStepLength =
+ std::min(secondLeapFrogLength, finalTrackLength);
+
+ CORSIKA_LOG_DEBUG(
+ "finalTrack(0)={}, finalTrack(1)={}, finalTrackLength={}, "
+ "secondLeapFrogLength={}, secondHalveStepLength={}, "
+ "secondLeapFrogLength-finalTrackLength={}, "
+ "secondHalveStepLength-finalTrackLength={}, "
+ "nextVol={}, transition={}",
+ finalTrack.getPosition(0).getCoordinates(),
+ finalTrack.getPosition(1).getCoordinates(), finalTrackLength,
+ secondLeapFrogLength, secondHalveStepLength,
+ secondLeapFrogLength - finalTrackLength,
+ secondHalveStepLength - finalTrackLength, fmt::ptr(finalTrackNextVolume),
+ switch_volume);
+
+ // perform the second halve-step
+ auto const new_direction_normalized = new_direction.normalized();
+ const Point finalPosition =
+ position_mid + new_direction_normalized * secondHalveStepLength;
+
+ const LengthType totalStep = firstHalveSteplength + secondHalveStepLength;
+ const auto delta_pos = finalPosition - initialPosition;
+ const auto distance = delta_pos.getNorm();
+
+ return std::make_tuple(
+ LineTrajectory(Line(initialPosition,
+ (distance == 0_m ? initialVelocity
+ : delta_pos.normalized() * absVelocity)),
+ distance / absVelocity, // straight distance
+ totalStep / absVelocity, // bend distance
+ initialVelocity,
+ new_direction_normalized * absVelocity), // trajectory
+ (switch_volume ? finalTrackNextVolume : volumeNode));
+ }
+
+ protected:
+ double firstFraction_;
+ };
+
+ } // namespace tracking_leapfrog_straight
+
+} // namespace corsika
--
GitLab