diff --git a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
new file mode 100644
index 0000000000000000000000000000000000000000..9ff574321434aa6c7c495ea3fc8db432e4a47879
--- /dev/null
+++ b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
@@ -0,0 +1,61 @@
+/*
+ * (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/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+namespace corsika {
+
+ inline Line LeapFrogTrajectory::getLine() const {
+ auto D = getPosition(1) - getPosition(0);
+ auto d = D.getNorm();
+ auto v = initialVelocity_;
+ if (d > 1_um) { // if trajectory is ultra-short, we do not
+ // re-calculate velocity, just use initial
+ // value. Otherwise, this is numerically unstable
+ v = D / d * getVelocity(0).getNorm();
+ }
+ return Line(getPosition(0), v);
+ }
+
+ inline Point LeapFrogTrajectory::getPosition(double const u) const {
+ Point position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
+ VelocityVector velocity =
+ initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ return position + velocity * timeStep_ * u / 2;
+ }
+
+ inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
+ return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+ }
+
+ inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
+ return getVelocity(u).normalized();
+ }
+
+ inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
+ return u * timeStep_ *
+ (double(getVelocity(u).getNorm() / initialVelocity_.getNorm()) + 1.0) / 2;
+ }
+
+ inline LengthType LeapFrogTrajectory::getLength(double const u) const {
+ return timeStep_ * initialVelocity_.getNorm() * u;
+ }
+
+ inline void LeapFrogTrajectory::setLength(LengthType const limit) {
+ if (initialVelocity_.getNorm() == 0_m / 1_s) setDuration(0_s);
+ setDuration(limit / initialVelocity_.getNorm());
+ }
+
+ inline void LeapFrogTrajectory::setDuration(TimeType const limit) { timeStep_ = limit; }
+
+} // namespace corsika
diff --git a/corsika/detail/framework/geometry/StraightTrajectory.inl b/corsika/detail/framework/geometry/StraightTrajectory.inl
new file mode 100644
index 0000000000000000000000000000000000000000..a8fef9b035a38745095976a17f6a0644b720a2a5
--- /dev/null
+++ b/corsika/detail/framework/geometry/StraightTrajectory.inl
@@ -0,0 +1,70 @@
+/*
+ * (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/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+namespace corsika {
+
+ inline VelocityVector StraightTrajectory::getVelocity(double const u) const {
+ return initialVelocity_ * (1 - u) + finalVelocity_ * u;
+ }
+
+ inline TimeType StraightTrajectory::getDuration(double const u) const {
+ return u * timeStep_;
+ }
+
+ inline LengthType StraightTrajectory::getLength(double const u) const {
+ if (timeLength_ == 0_s) return 0_m;
+ if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s)
+ return std::numeric_limits<LengthType::value_type>::infinity() * 1_m;
+ return getDistance(u) * timeStep_ / timeLength_;
+ }
+
+ inline void StraightTrajectory::setLength(LengthType const limit) {
+ setDuration(line_.getTimeFromArclength(limit));
+ }
+
+ inline void StraightTrajectory::setDuration(TimeType const limit) {
+ if (timeStep_ == 0_s) {
+ timeLength_ = 0_s;
+ setFinalVelocity(getVelocity(0));
+ timeStep_ = limit;
+ } else {
+ // for infinite steps there can't be a difference between
+ // curved and straight trajectory, this is fundamentally
+ // undefined: assume they are the same (which, i.e. is always correct for a
+ // straight line trajectory).
+ //
+ // Final note: only straight-line trajectories should have
+ // infinite steps! Everything else is ill-defined.
+ if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s ||
+ timeLength_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s) {
+ timeLength_ = limit;
+ timeStep_ = limit;
+ // ...and don't touch velocity
+ } else {
+ const double scale = limit / timeStep_;
+ timeLength_ *= scale;
+ setFinalVelocity(getVelocity(scale));
+ timeStep_ = limit;
+ }
+ }
+ }
+
+ inline LengthType StraightTrajectory::getDistance(double const u) const {
+ assert(u <= 1);
+ assert(u >= 0);
+ return line_.getArcLength(0 * second, u * timeLength_);
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/tracking/Intersect.inl b/corsika/detail/modules/tracking/Intersect.inl
new file mode 100644
index 0000000000000000000000000000000000000000..9c30ff87d1156bfdfd7c1a94d2ba32b4b771c73b
--- /dev/null
+++ b/corsika/detail/modules/tracking/Intersect.inl
@@ -0,0 +1,111 @@
+/*
+ * (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/core/PhysicalUnits.hpp>
+#include <corsika/framework/logging/Logging.hpp>
+#include <corsika/framework/geometry/Intersections.hpp>
+
+#include <limits>
+
+namespace corsika {
+
+ template <typename TDerived>
+ template <typename TParticle>
+ inline auto Intersect<TDerived>::nextIntersect(TParticle const& particle,
+ TimeType const step_limit) const {
+
+ typedef typename std::remove_reference<decltype(*particle.getNode())>::type node_type;
+
+ node_type& volumeNode =
+ *particle.getNode(); // current "logical" node, from previous tracking step
+
+ CORSIKA_LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
+ volumeNode.getVolume().contains(particle.getPosition()));
+
+ // start values:
+ TimeType minTime = step_limit;
+ node_type* minNode = &volumeNode;
+
+ // determine the first geometric collision with any other Volume boundary
+
+ // first check, where we leave the current volume
+ // this assumes our convention, that all Volume primitives must be convex
+ // thus, the last entry is always the exit point
+ const Intersections time_intersections_curr =
+ TDerived::intersect(particle, volumeNode);
+ CORSIKA_LOG_TRACE("curr node {}, parent node {}, hasIntersections={} ",
+ fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
+ time_intersections_curr.hasIntersections());
+ if (time_intersections_curr.hasIntersections()) {
+ CORSIKA_LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
+ time_intersections_curr.getEntry() / 1_s,
+ time_intersections_curr.getExit() / 1_s);
+ if (time_intersections_curr.getExit() <= minTime) {
+ minTime =
+ time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
+ minNode = volumeNode.getParent();
+ }
+ }
+
+ // where do we collide with any of the next-tree-level volumes
+ // entirely contained by currentLogicalVolumeNode
+ for (const auto& node : volumeNode.getChildNodes()) {
+
+ const Intersections time_intersections = TDerived::intersect(particle, *node);
+ if (!time_intersections.hasIntersections()) { continue; }
+ CORSIKA_LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
+ fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ time_intersections.getExit() / 1_s);
+
+ const auto t_entry = time_intersections.getEntry();
+ const auto t_exit = time_intersections.getExit();
+ CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+ t_entry <= minTime);
+ // note, theoretically t can even be smaller than 0 since we
+ // KNOW we can't yet be in this volume yet, so we HAVE TO
+ // enter it IF exit point is not also in the "past", AND if
+ // extry point is [[much much]] closer than exit point
+ // (because we might have already numerically "exited" it)!
+ if (t_exit > 0_s && t_entry <= minTime &&
+ -t_entry < t_exit) { // protection agains numerical problem, when we already
+ // _exited_ before
+ // enter chile volume here
+ minTime = t_entry;
+ minNode = node.get();
+ }
+ }
+
+ // these are volumes from the previous tree-level that are cut-out partly from the
+ // current volume
+ for (node_type* node : volumeNode.getExcludedNodes()) {
+
+ const Intersections time_intersections = TDerived::intersect(particle, *node);
+ if (!time_intersections.hasIntersections()) { continue; }
+ CORSIKA_LOG_DEBUG(
+ "intersection times with exclusion volume {} : enter {} s, exit {} s",
+ fmt::ptr(node), time_intersections.getEntry() / 1_s,
+ time_intersections.getExit() / 1_s);
+ const auto t_entry = time_intersections.getEntry();
+ const auto t_exit = time_intersections.getExit();
+ CORSIKA_LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
+ t_entry <= minTime);
+ // note, theoretically t can even be smaller than 0 since we
+ // KNOW we can't yet be in this volume yet, so we HAVE TO
+ // enter it IF exit point is not also in the "past"!
+ if (t_exit > 0_s && t_entry <= minTime) { // enter volumen child here
+ minTime = t_entry;
+ minNode = node;
+ }
+ }
+ CORSIKA_LOG_TRACE("t-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
+ return std::make_tuple(minTime, minNode);
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
new file mode 100644
index 0000000000000000000000000000000000000000..c57267525342b26205ddc3d8737252e8dbde6673
--- /dev/null
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogCurved.inl
@@ -0,0 +1,259 @@
+/*
+ * (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/LeapFrogTrajectory.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 {
+
+ template <typename TParticle>
+ inline auto make_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);
+ }
+
+ template <typename TParticle>
+ inline auto Tracking::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>
+ inline Intersections Tracking::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<TParticle, TMedium>(particle, sphere,
+ medium);
+ }
+
+ bool const numericallyInside = sphere.contains(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>
+ inline Intersections Tracking::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)");
+ }
+
+ template <typename TParticle>
+ inline auto Tracking::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
+ }
+
+ } // namespace tracking_leapfrog_curved
+
+} // namespace corsika
diff --git a/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
new file mode 100644
index 0000000000000000000000000000000000000000..132a87e5d80abc129bc3375b4cab921c701a565c
--- /dev/null
+++ b/corsika/detail/modules/tracking/TrackingLeapFrogStraight.inl
@@ -0,0 +1,186 @@
+/*
+ * (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/StraightTrajectory.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 {
+
+ template <typename Particle>
+ inline auto Tracking::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(
+ StraightTrajectory(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(
+ StraightTrajectory(
+ 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));
+ }
+
+ } // namespace tracking_leapfrog_straight
+
+} // namespace corsika
diff --git a/corsika/detail/stack/SimpleStack.inl b/corsika/detail/stack/SimpleStack.inl
new file mode 100644
index 0000000000000000000000000000000000000000..4b0e13ecb29592049deb275bf37b1ecddc7b4c48
--- /dev/null
+++ b/corsika/detail/stack/SimpleStack.inl
@@ -0,0 +1,93 @@
+/*
+ * (c) Copyright 2018 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/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <string>
+#include <tuple>
+#include <vector>
+
+namespace corsika {
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setParticleData(
+ std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
+ this->setPID(std::get<0>(v));
+ this->setEnergy(std::get<1>(v));
+ this->setMomentum(std::get<2>(v));
+ this->setPosition(std::get<3>(v));
+ this->setTime(std::get<4>(v));
+ }
+
+ template <typename StackIteratorInterface>
+ void ParticleInterface<StackIteratorInterface>::setParticleData(
+ ParticleInterface<StackIteratorInterface> const&,
+ std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v) {
+ this->setPID(std::get<0>(v));
+ this->setEnergy(std::get<1>(v));
+ this->setMomentum(std::get<2>(v));
+ this->setPosition(std::get<3>(v));
+ this->setTime(std::get<4>(v));
+ }
+
+ inline void SimpleStackImpl::clear() {
+ dataPID_.clear();
+ dataE_.clear();
+ momentum_.clear();
+ position_.clear();
+ time_.clear();
+ }
+
+ inline void SimpleStackImpl::copy(size_t i1, size_t i2) {
+ dataPID_[i2] = dataPID_[i1];
+ dataE_[i2] = dataE_[i1];
+ momentum_[i2] = momentum_[i1];
+ position_[i2] = position_[i1];
+ time_[i2] = time_[i1];
+ }
+
+ inline void SimpleStackImpl::swap(size_t i1, size_t i2) {
+ std::swap(dataPID_[i2], dataPID_[i1]);
+ std::swap(dataE_[i2], dataE_[i1]);
+ std::swap(momentum_[i2], momentum_[i1]);
+ std::swap(position_[i2], position_[i1]);
+ std::swap(time_[i2], time_[i1]);
+ }
+
+ inline void SimpleStackImpl::incrementSize() {
+ dataPID_.push_back(Code::Unknown);
+ dataE_.push_back(0 * electronvolt);
+
+ CoordinateSystemPtr const& dummyCS = get_root_CoordinateSystem();
+
+ momentum_.push_back(
+ MomentumVector(dummyCS, {0 * electronvolt, 0 * electronvolt, 0 * electronvolt}));
+
+ position_.push_back(Point(dummyCS, {0 * meter, 0 * meter, 0 * meter}));
+ time_.push_back(0 * second);
+ }
+
+ inline void SimpleStackImpl::decrementSize() {
+ if (dataE_.size() > 0) {
+ dataPID_.pop_back();
+ dataE_.pop_back();
+ momentum_.pop_back();
+ position_.pop_back();
+ time_.pop_back();
+ }
+ }
+
+} // namespace corsika
diff --git a/corsika/framework/geometry/LeapFrogTrajectory.hpp b/corsika/framework/geometry/LeapFrogTrajectory.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6b67bf6345a5c79400d72d9ec16e13afa380e13f
--- /dev/null
+++ b/corsika/framework/geometry/LeapFrogTrajectory.hpp
@@ -0,0 +1,83 @@
+/*
+ * (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/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+namespace corsika {
+
+ /**
+ * The LeapFrogTrajectory stores information on one leap-frog step.
+ *
+ * The leap-frog algorithm uses a half-step and is used in magnetic
+ * field tracking. The LeapFrogTrajectory will solve the leap-frog
+ * algorithm equation for a given constant $k$ that has to be
+ * specified during construction (essentially fixing the magnetic
+ * field). Thus, different steps (length) can be dynamically
+ * generated here. The velocity vector will correctly point into the
+ * direction as calculated by the algorithm for any steplength, or
+ * intermediate position.
+ *
+ **/
+
+ class LeapFrogTrajectory {
+
+ public:
+ LeapFrogTrajectory() = delete;
+ LeapFrogTrajectory(LeapFrogTrajectory const&) = default;
+ LeapFrogTrajectory(LeapFrogTrajectory&&) = default;
+ LeapFrogTrajectory& operator=(LeapFrogTrajectory const&) = delete;
+
+ LeapFrogTrajectory(Point const& pos, VelocityVector const& initialVelocity,
+ MagneticFieldVector const& Bfield,
+ decltype(square(meter) / (square(second) * volt)) const k,
+ TimeType const timeStep) // leap-from total length
+ : initialPosition_(pos)
+ , initialVelocity_(initialVelocity)
+ , initialDirection_(initialVelocity.normalized())
+ , magneticfield_(Bfield)
+ , k_(k)
+ , timeStep_(timeStep) {}
+
+ Line getLine() const;
+
+ Point getPosition(double const u) const;
+
+ VelocityVector getVelocity(double const u) const;
+
+ DirectionVector getDirection(double const u) const;
+
+ ///! duration along potentially bend trajectory
+ TimeType getDuration(double const u = 1) const;
+
+ ///! total length along potentially bend trajectory
+ LengthType getLength(double const u = 1) const;
+
+ ///! set new duration along potentially bend trajectory.
+ void setLength(LengthType const limit);
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ void setDuration(TimeType const limit);
+
+ private:
+ Point initialPosition_;
+ VelocityVector initialVelocity_;
+ DirectionVector initialDirection_;
+ MagneticFieldVector magneticfield_;
+ decltype(square(meter) / (square(second) * volt)) k_;
+ TimeType timeStep_;
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/LeapFrogTrajectory.inl>
diff --git a/corsika/framework/geometry/StraightTrajectory.hpp b/corsika/framework/geometry/StraightTrajectory.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..b05746269b017f7a80fe5c5807ca93c7afe336ab
--- /dev/null
+++ b/corsika/framework/geometry/StraightTrajectory.hpp
@@ -0,0 +1,124 @@
+/*
+ * (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/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+namespace corsika {
+
+ /**
+ *
+ * A Trajectory is a description of a momvement of an object in
+ * three-dimensional space that describes the trajectory (connection
+ * between two Points in space), as well as the direction of motion
+ * at any given point.
+ *
+ * A Trajectory has a start `0` and an end `1`, where
+ * e.g. getPosition(0) returns the start point and getDirection(1)
+ * the direction of motion at the end. Values outside 0...1 are not
+ * defined.
+ *
+ * A Trajectory has a length in [m], getLength, a duration in [s], getDuration.
+ *
+ * Note: so far it is assumed that the speed (d|vec{r}|/dt) between
+ * start and end does not change and is constant for the entire
+ * Trajectory.
+ *
+ **/
+
+ class StraightTrajectory {
+
+ public:
+ StraightTrajectory() = delete;
+ StraightTrajectory(StraightTrajectory const&) = default;
+ StraightTrajectory(StraightTrajectory&&) = default;
+ StraightTrajectory& operator=(StraightTrajectory const&) = delete;
+
+ /**
+ * \param theLine The geometric Line object that represents a straight-line
+ * connection
+ *
+ * \param timeLength The time duration to traverse the straight trajectory
+ * in units of TimeType
+ */
+ StraightTrajectory(Line const& theLine, TimeType timeLength)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeLength)
+ , initialVelocity_(theLine.getVelocity(TimeType::zero()))
+ , finalVelocity_(theLine.getVelocity(timeLength)) {}
+
+ /**
+ * \param theLine The geometric Line object that represents a straight-line
+ * connection
+ *
+ * \param timeLength The time duration to traverse the straight trajectory
+ * in units of TimeType
+ *
+ * \param timeStep Time duration to folow eventually curved
+ * trajectory in units of TimesType
+ *
+ * \param initialV Initial velocity vector at
+ * start of trajectory
+ *
+ * \param finalV Final velocity vector at start of trajectory
+ */
+ StraightTrajectory(Line const& theLine,
+ TimeType const timeLength, // length of theLine (straight)
+ TimeType const timeStep, // length of bend step (curved)
+ VelocityVector const& initialV, VelocityVector const& finalV)
+ : line_(theLine)
+ , timeLength_(timeLength)
+ , timeStep_(timeStep)
+ , initialVelocity_(initialV)
+ , finalVelocity_(finalV) {}
+
+ Line const& getLine() const { return line_; }
+
+ Point getPosition(double const u) const { return line_.getPosition(timeLength_ * u); }
+
+ VelocityVector getVelocity(double const u) const;
+
+ DirectionVector getDirection(double const u) const {
+ return getVelocity(u).normalized();
+ }
+
+ ///! duration along potentially bend trajectory
+ TimeType getDuration(double const u = 1) const;
+
+ ///! total length along potentially bend trajectory
+ LengthType getLength(double const u = 1) const;
+
+ ///! set new duration along potentially bend trajectory.
+ void setLength(LengthType const limit);
+
+ ///! set new duration along potentially bend trajectory.
+ // Scale other properties by "limit/timeLength_"
+ void setDuration(TimeType const limit);
+
+ protected:
+ ///! total length along straight trajectory
+ LengthType getDistance(double const u) const;
+
+ void setFinalVelocity(VelocityVector const& v) { finalVelocity_ = v; }
+
+ private:
+ Line line_;
+ TimeType timeLength_; ///! length of straight step (shortest connecting line)
+ TimeType timeStep_; ///! length of bend step (curved)
+ VelocityVector initialVelocity_;
+ VelocityVector finalVelocity_;
+ };
+
+} // namespace corsika
+
+#include <corsika/detail/framework/geometry/StraightTrajectory.inl>
diff --git a/corsika/stack/SimpleStack.hpp b/corsika/stack/SimpleStack.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a40d379b75f95ba6929c55bc2ab3f9afbb083899
--- /dev/null
+++ b/corsika/stack/SimpleStack.hpp
@@ -0,0 +1,176 @@
+/*
+ * (c) Copyright 2018 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/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/stack/Stack.hpp>
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+#include <string>
+#include <tuple>
+#include <vector>
+
+namespace corsika {
+
+ /**
+ * Example of a particle object on the stack.
+ */
+
+ template <typename StackIteratorInterface>
+ struct ParticleInterface : public ParticleBase<StackIteratorInterface> {
+
+ private:
+ typedef ParticleBase<StackIteratorInterface> super_type;
+
+ public:
+ std::string asString() const {
+ return fmt::format("particle: i={}, PID={}, E={}GeV", super_type::getIndex(),
+ get_name(this->getPID()), this->getEnergy() / 1_GeV);
+ }
+
+ void setParticleData(
+ std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
+
+ void setParticleData(
+ ParticleInterface<StackIteratorInterface> const&,
+ std::tuple<Code, HEPEnergyType, MomentumVector, Point, TimeType> const& v);
+
+ /// individual setters
+ void setPID(Code const id) {
+ super_type::getStackData().setPID(super_type::getIndex(), id);
+ }
+ void setEnergy(HEPEnergyType const& e) {
+ super_type::getStackData().setEnergy(super_type::getIndex(), e);
+ }
+ void setMomentum(MomentumVector const& v) {
+ super_type::getStackData().setMomentum(super_type::getIndex(), v);
+ }
+ void setPosition(Point const& v) {
+ super_type::getStackData().setPosition(super_type::getIndex(), v);
+ }
+ void setTime(TimeType const& v) {
+ super_type::getStackData().setTime(super_type::getIndex(), v);
+ }
+
+ /// individual getters
+ Code getPID() const {
+ return super_type::getStackData().getPID(super_type::getIndex());
+ }
+ HEPEnergyType getEnergy() const {
+ return super_type::getStackData().getEnergy(super_type::getIndex());
+ }
+ MomentumVector getMomentum() const {
+ return super_type::getStackData().getMomentum(super_type::getIndex());
+ }
+ Point getPosition() const {
+ return super_type::getStackData().getPosition(super_type::getIndex());
+ }
+ TimeType getTime() const {
+ return super_type::getStackData().getTime(super_type::getIndex());
+ }
+ /**
+ * @name derived quantities
+ *
+ * @{
+ */
+ DirectionVector getDirection() const {
+ return this->getMomentum() / this->getEnergy();
+ }
+
+ HEPMassType getMass() const { return get_mass(this->getPID()); }
+
+ int16_t getChargeNumber() const { return get_charge_number(this->getPID()); }
+ ///@}
+ };
+
+ /**
+ * Memory implementation of the most simple (stupid) particle stack object.
+ *
+ */
+
+ class SimpleStackImpl {
+
+ public:
+ typedef std::vector<Code> code_vector_type;
+ typedef std::vector<HEPEnergyType> energy_vector_type;
+ typedef std::vector<Point> point_vector_type;
+ typedef std::vector<TimeType> time_vector_type;
+ typedef std::vector<MomentumVector> momentum_vector_type;
+
+ SimpleStackImpl() = default;
+
+ SimpleStackImpl(SimpleStackImpl const& other) = default;
+
+ SimpleStackImpl(SimpleStackImpl&& other) = default;
+
+ SimpleStackImpl& operator=(SimpleStackImpl const& other) = default;
+
+ SimpleStackImpl& operator=(SimpleStackImpl&& other) = default;
+
+ void dump() const {}
+
+ inline void clear();
+
+ unsigned int getSize() const { return dataPID_.size(); }
+ unsigned int getCapacity() const { return dataPID_.size(); }
+
+ void setPID(size_t i, Code const id) { dataPID_[i] = id; }
+ void setEnergy(size_t i, HEPEnergyType const& e) { dataE_[i] = e; }
+ void setMomentum(size_t i, MomentumVector const& v) { momentum_[i] = v; }
+ void setPosition(size_t i, Point const& v) { position_[i] = v; }
+ void setTime(size_t i, TimeType const& v) { time_[i] = v; }
+
+ Code getPID(size_t i) const { return dataPID_[i]; }
+
+ HEPEnergyType getEnergy(size_t i) const { return dataE_[i]; }
+
+ MomentumVector getMomentum(size_t i) const { return momentum_[i]; }
+
+ Point getPosition(size_t i) const { return position_[i]; }
+ TimeType getTime(size_t i) const { return time_[i]; }
+
+ HEPEnergyType getDataE(size_t i) const { return dataE_[i]; }
+
+ void setDataE(size_t i, HEPEnergyType const& dataE) { dataE_[i] = dataE; }
+
+ Code getDataPid(size_t i) const { return dataPID_[i]; }
+
+ void setDataPid(size_t i, Code const& dataPid) { dataPID_[i] = dataPid; }
+ /**
+ * Function to copy particle at location i2 in stack to i1
+ */
+ inline void copy(size_t i1, size_t i2);
+
+ /**
+ * Function to copy particle at location i2 in stack to i1
+ */
+ inline void swap(size_t i1, size_t i2);
+
+ inline void incrementSize();
+ inline void decrementSize();
+
+ private:
+ /// the actual memory to store particle data
+ code_vector_type dataPID_;
+ energy_vector_type dataE_;
+ momentum_vector_type momentum_;
+ point_vector_type position_;
+ time_vector_type time_;
+
+ }; // end class SimpleStackImpl
+
+ typedef Stack<SimpleStackImpl, ParticleInterface> SimpleStack;
+
+} // namespace corsika
+
+#include <corsika/detail/stack/SimpleStack.inl>