diff --git a/Framework/Cascade/Cascade_interpolation.h b/Framework/Cascade/Cascade_interpolation.h
new file mode 100644
index 0000000000000000000000000000000000000000..dccd3f8e705ca264fa27e7b637bdc4ab332b74b8
--- /dev/null
+++ b/Framework/Cascade/Cascade_interpolation.h
@@ -0,0 +1,341 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_corsika_cascade_Cascade_h_
+#define _include_corsika_cascade_Cascade_h_
+
+#include <corsika/environment/Environment.h>
+#include <corsika/process/ProcessReturn.h>
+#include <corsika/random/ExponentialDistribution.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/random/UniformRealDistribution.h>
+#include <corsika/stack/SecondaryView.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+
+/* see Issue 161, we need to include SetupStack only because we need
+ to globally define StackView. This is clearly not nice and should
+ be changed, when possible. It might be that StackView needs to be
+ templated in Cascade, but this would be even worse... so we don't
+ do that until it is really needed.
+ */
+#include <corsika/setup/SetupStack.h>
+
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <limits>
+#include <type_traits>
+
+#include <boost/type_index.hpp>
+using boost::typeindex::type_id_with_cvr;
+
+/**
+ * The cascade namespace assembles all objects needed to simulate full particles cascades.
+ */
+
+namespace corsika::cascade {
+
+ /**
+ * \class Cascade
+ *
+ * The Cascade class is constructed from template arguments making
+ * it very versatile. Via the template arguments physics models are
+ * plugged into the cascade simulation.
+ *
+ * <b>TTracking</b> must be a class according to the
+ * TrackingInterface providing the functions:
+ * <code>auto GetTrack(Particle const& p)</auto>,
+ * with the return type <code>geometry::Trajectory<corsika::geometry::Line>
+ * </code>
+ *
+ * <b>TProcessList</b> must be a ProcessSequence. *
+ * <b>Stack</b> is the storage object for particle data, i.e. with
+ * Particle class type <code>Stack::ParticleType</code>
+ *
+ *
+ */
+
+ template <typename TTracking, typename TProcessList, typename TStack,
+ /*
+ TStackView is needed as template parameter because of issue 161 and the
+ inability of clang to understand "MakeView" so far.
+ */
+ typename TStackView = corsika::setup::StackView>
+ class Cascade {
+ using Particle = typename TStack::ParticleType;
+ using VolumeTreeNode =
+ std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
+ using MediumInterface = typename VolumeTreeNode::IModelProperties;
+
+ // we only want fully configured objects
+ Cascade() = delete;
+
+ public:
+ /**
+ * Cascade class cannot be default constructed, but needs a valid
+ * list of physics processes for configuration at construct time.
+ */
+ Cascade(corsika::environment::Environment<MediumInterface> const& env, TTracking& tr,
+ TProcessList& pl, TStack& stack)
+ : fEnvironment(env)
+ , fTracking(tr)
+ , fProcessSequence(pl)
+ , fStack(stack) {}
+
+ /**
+ * The Init function is called before the actual cascade simulations.
+ * All components of the Cascade simulation must be configured here.
+ */
+ void Init() {
+ fProcessSequence.Init();
+ fStack.Init();
+ }
+
+ /**
+ * set the nodes for all particles on the stack according to their numerical
+ * position
+ */
+ void SetNodes() {
+ std::for_each(fStack.begin(), fStack.end(), [&](auto& p) {
+ auto const* numericalNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ p.SetNode(numericalNode);
+ });
+ }
+
+ /**
+ * The Run function is the main simulation loop, which processes
+ * particles from the Stack until the Stack is empty.
+ */
+ void Run() {
+ SetNodes();
+
+ while (!fStack.IsEmpty()) {
+ while (!fStack.IsEmpty()) {
+ auto pNext = fStack.GetNextParticle();
+ std::cout << "========= next: " << pNext.GetPID() << std::endl;
+ Step(pNext);
+ std::cout << "========= stack ============" << std::endl;
+ fProcessSequence.DoStack(fStack);
+ }
+ // do cascade equations, which can put new particles on Stack,
+ // thus, the double loop
+ // DoCascadeEquations();
+ }
+ }
+
+ /**
+ * Force an interaction of the top particle of the stack at its current position.
+ * Note that SetNodes() or an equivalent procedure needs to be called first if you
+ * want to call forceInteraction() for the primary interaction.
+ */
+ void forceInteraction() {
+ std::cout << "forced interaction!" << std::endl;
+ auto vParticle = fStack.GetNextParticle();
+ TStackView secondaries(vParticle);
+ auto projectile = secondaries.GetProjectile();
+ interaction(vParticle, projectile);
+ fProcessSequence.DoSecondaries(secondaries);
+ vParticle.Delete(); // todo: this should be reviewed, see below
+ }
+
+ private:
+ /**
+ * The Step function is executed for each particle from the
+ * stack. It will calcualte geometric transport of the particles,
+ * and apply continuous and stochastic processes to it, which may
+ * lead to energy losses, scattering, absorption, decays and the
+ * production of secondary particles.
+ *
+ * New particles produced in one step are subject to further
+ * processing, e.g. thinning, etc.
+ */
+ void Step(Particle& vParticle) {
+ using namespace corsika;
+ using namespace corsika::units::si;
+
+ // determine combined total interaction length (inverse)
+ InverseGrammageType const total_inv_lambda =
+ fProcessSequence.GetTotalInverseInteractionLength(vParticle);
+
+ // sample random exponential step length in grammage
+ corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
+ GrammageType const next_interact = expDist(fRNG);
+
+ std::cout << "total_inv_lambda=" << total_inv_lambda
+ << ", next_interact=" << next_interact << std::endl;
+
+ auto const* currentLogicalNode = vParticle.GetNode();
+
+ // assert that particle stays outside void Universe if it has no
+ // model properties set
+ assert(currentLogicalNode != &*fEnvironment.GetUniverse() ||
+ fEnvironment.GetUniverse()->HasModelProperties());
+
+ // determine combined total inverse decay time
+ InverseTimeType const total_inv_lifetime =
+ fProcessSequence.GetTotalInverseLifetime(vParticle);
+
+ // sample random exponential decay time
+ corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
+ TimeType const next_decay = expDistDecay(fRNG);
+ std::cout << "total_inv_lifetime=" << total_inv_lifetime
+ << ", next_decay=" << next_decay << std::endl;
+
+ // convert next_decay from time to length [m]
+ LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
+ vParticle.GetEnergy() * units::constants::c;
+
+ // determine geometric tracking
+ auto [step, geomMaxLength, nextVol, magMaxLength, directionBefore, directionAfter] = fTracking.GetTrack(vParticle);
+ [[maybe_unused]] auto const& dummy_nextVol = nextVol;
+
+ // convert next_step from grammage to length
+ LengthType const distance_interact =
+ currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
+ next_interact);
+
+ // determine the maximum geometric step length
+ LengthType const distance_max = fProcessSequence.MaxStepLength(vParticle, step);
+ std::cout << "distance_max=" << distance_max << std::endl;
+
+ // take minimum of geometry, interaction, decay for next step
+ auto const min_distance = std::min(
+ {distance_interact, distance_decay, distance_max, geomMaxLength, magMaxLength});
+
+ std::cout << " move particle by : " << min_distance << std::endl;
+
+ // here the particle is actually moved along the trajectory to new position:
+ // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
+ vParticle.SetPosition(step.PositionFromArclength(min_distance));
+ // .... also update time, momentum, direction, ...
+ vParticle.SetMomentum((directionBefore * (1 - min_distance / magMaxLength) +
+ directionAfter * min_distance /magMaxLength) * vParticle.GetMomentum().GetNorm());
+ vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
+
+ step.LimitEndTo(min_distance);
+
+ // apply all continuous processes on particle + track
+ process::EProcessReturn status = fProcessSequence.DoContinuous(vParticle, step);
+
+ if (status == process::EProcessReturn::eParticleAbsorbed) {
+ std::cout << "Cascade: delete absorbed particle " << vParticle.GetPID() << " "
+ << vParticle.GetEnergy() / 1_GeV << "GeV" << std::endl;
+ vParticle.Delete();
+ return;
+ }
+
+ std::cout << "sth. happening before geometric limit ? "
+ << ((min_distance < geomMaxLength) ? "yes" : "no") << std::endl;
+
+ if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
+
+ // check whether decay or interaction limits this step the
+ // outcome of decay or interaction MAY be a) new particles in
+ // secondaries, b) the projectile particle deleted (or
+ // changed)
+
+ TStackView secondaries(vParticle);
+
+ if (min_distance != distance_max && min_distance != magMaxLength) {
+ /*
+ Create SecondaryView object on Stack. The data container
+ remains untouched and identical, and 'projectil' is identical
+ to 'vParticle' above this line. However,
+ projectil.AddSecondaries populate the SecondaryView, which can
+ then be used afterwards for further processing. Thus: it is
+ important to use projectle (and not vParticle) for Interaction,
+ and Decay!
+ */
+
+ [[maybe_unused]] auto projectile = secondaries.GetProjectile();
+
+ if (min_distance == distance_interact) {
+ interaction(vParticle, projectile);
+ } else {
+ assert(min_distance == distance_decay);
+ decay(vParticle, projectile);
+ // make sure particle actually did decay if it should have done so
+ if (secondaries.GetSize() == 1 &&
+ projectile.GetPID() == secondaries.GetNextParticle().GetPID())
+ throw std::runtime_error("Cascade::Step: Particle decays into itself!");
+ }
+
+ fProcessSequence.DoSecondaries(secondaries);
+ vParticle.Delete(); // todo: this should be reviewed. Where
+ // exactly are particles best deleted, and
+ // where they should NOT be
+ // deleted... maybe Delete function should
+ // be "protected" and not accessible to physics
+
+ } else { // step-length limitation within volume
+
+ std::cout << "step-length limitation" << std::endl;
+ fProcessSequence.DoSecondaries(secondaries);
+ }
+
+ [[maybe_unused]] auto const assertion = [&] {
+ auto const* numericalNodeAfterStep =
+ fEnvironment.GetUniverse()->GetContainingNode(vParticle.GetPosition());
+ return numericalNodeAfterStep == currentLogicalNode;
+ };
+
+ assert(assertion()); // numerical and logical nodes don't match
+ } else { // boundary crossing, step is limited by volume boundary
+ std::cout << "boundary crossing! next node = " << nextVol << std::endl;
+ vParticle.SetNode(nextVol);
+ // DoBoundary may delete the particle (or not)
+ fProcessSequence.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
+ }
+ }
+
+ auto decay(Particle& particle,
+ decltype(std::declval<TStackView>().GetProjectile()) projectile) {
+ std::cout << "decay" << std::endl;
+ units::si::InverseTimeType const actual_decay_time =
+ fProcessSequence.GetTotalInverseLifetime(particle);
+
+ random::UniformRealDistribution<units::si::InverseTimeType> uniDist(
+ actual_decay_time);
+ const auto sample_process = uniDist(fRNG);
+ units::si::InverseTimeType inv_decay_count = units::si::InverseTimeType::zero();
+ return fProcessSequence.SelectDecay(particle, projectile, sample_process,
+ inv_decay_count);
+ }
+
+ auto interaction(Particle& particle,
+ decltype(std::declval<TStackView>().GetProjectile()) projectile) {
+ std::cout << "collide" << std::endl;
+
+ units::si::InverseGrammageType const current_inv_length =
+ fProcessSequence.GetTotalInverseInteractionLength(particle);
+
+ random::UniformRealDistribution<units::si::InverseGrammageType> uniDist(
+ current_inv_length);
+ const auto sample_process = uniDist(fRNG);
+ auto inv_lambda_count = units::si::InverseGrammageType::zero();
+ return fProcessSequence.SelectInteraction(particle, projectile, sample_process,
+ inv_lambda_count);
+ }
+
+ private:
+ corsika::environment::Environment<MediumInterface> const& fEnvironment;
+ TTracking& fTracking;
+ TProcessList& fProcessSequence;
+ TStack& fStack;
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
+ }; // namespace corsika::cascade
+
+} // namespace corsika::cascade
+
+#endif
diff --git a/Framework/Utilities/quartic.cpp b/Framework/Utilities/quartic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f11a8e083e887a3f11ac05dba394f232fdcedd3f
--- /dev/null
+++ b/Framework/Utilities/quartic.cpp
@@ -0,0 +1,152 @@
+/***************************************************************************
+ * Copyright (C) 2016 by Саша Миленковић *
+ * sasa.milenkovic.xyz@gmail.com *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#include <cmath>
+#include "quartic.h"
+
+//---------------------------------------------------------------------------
+// solve cubic equation x^3 + a*x^2 + b*x + c
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+// 2 real roots: x[0], x[1], return 2
+// 1 real root : x[0], x[1] ± i*x[2], return 1
+unsigned int solveP3(double *x,double a,double b,double c) {
+ double a2 = a*a;
+ double q = (a2 - 3*b)/9;
+ double r = (a*(2*a2-9*b) + 27*c)/54;
+ double r2 = r*r;
+ double q3 = q*q*q;
+ double A,B;
+ if(r2<q3)
+ {
+ double t=r/sqrt(q3);
+ if( t<-1) t=-1;
+ if( t> 1) t= 1;
+ t=acos(t);
+ a/=3; q=-2*sqrt(q);
+ x[0]=q*cos(t/3)-a;
+ x[1]=q*cos((t+M_2PI)/3)-a;
+ x[2]=q*cos((t-M_2PI)/3)-a;
+ return 3;
+ }
+ else
+ {
+ A =-pow(fabs(r)+sqrt(r2-q3),1./3);
+ if( r<0 ) A=-A;
+ B = (0==A ? 0 : q/A);
+
+ a/=3;
+ x[0] =(A+B)-a;
+ x[1] =-0.5*(A+B)-a;
+ x[2] = 0.5*sqrt(3.)*(A-B);
+ if(fabs(x[2])<eps) { x[2]=x[1]; return 2; }
+
+ return 1;
+ }
+}
+
+//---------------------------------------------------------------------------
+// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+// Attention - this function returns dynamically allocated array. It has to be released afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d)
+{
+ double a3 = -b;
+ double b3 = a*c -4.*d;
+ double c3 = -a*a*d - c*c + 4.*b*d;
+
+ // cubic resolvent
+ // y^3 − b*y^2 + (ac−4d)*y − a^2*d−c^2+4*b*d = 0
+
+ double x3[3];
+ unsigned int iZeroes = solveP3(x3, a3, b3, c3);
+
+ double q1, q2, p1, p2, D, sqD, y;
+
+ y = x3[0];
+ // The essence - choosing Y with maximal absolute value.
+ if(iZeroes != 1)
+ {
+ if(fabs(x3[1]) > fabs(y)) y = x3[1];
+ if(fabs(x3[2]) > fabs(y)) y = x3[2];
+ }
+
+ // h1+h2 = y && h1*h2 = d <=> h^2 -y*h + d = 0 (h === q)
+
+ D = y*y - 4*d;
+ if(fabs(D) < eps) //in other words - D==0
+ {
+ q1 = q2 = y * 0.5;
+ // g1+g2 = a && g1+g2 = b-y <=> g^2 - a*g + b-y = 0 (p === g)
+ D = a*a - 4*(b-y);
+ if(fabs(D) < eps) //in other words - D==0
+ p1 = p2 = a * 0.5;
+
+ else
+ {
+ sqD = sqrt(D);
+ p1 = (a + sqD) * 0.5;
+ p2 = (a - sqD) * 0.5;
+ }
+ }
+ else
+ {
+ sqD = sqrt(D);
+ q1 = (y + sqD) * 0.5;
+ q2 = (y - sqD) * 0.5;
+ // g1+g2 = a && g1*h2 + g2*h1 = c ( && g === p ) Krammer
+ p1 = (a*q1-c)/(q1-q2);
+ p2 = (c-a*q2)/(q1-q2);
+ }
+
+ DComplex* retval = new DComplex[4];
+
+ // solving quadratic eq. - x^2 + p1*x + q1 = 0
+ D = p1*p1 - 4*q1;
+ if(D < 0.0)
+ {
+ retval[0].real( -p1 * 0.5 );
+ retval[0].imag( sqrt(-D) * 0.5 );
+ retval[1] = std::conj(retval[0]);
+ }
+ else
+ {
+ sqD = sqrt(D);
+ retval[0].real( (-p1 + sqD) * 0.5 );
+ retval[1].real( (-p1 - sqD) * 0.5 );
+ }
+
+ // solving quadratic eq. - x^2 + p2*x + q2 = 0
+ D = p2*p2 - 4*q2;
+ if(D < 0.0)
+ {
+ retval[2].real( -p2 * 0.5 );
+ retval[2].imag( sqrt(-D) * 0.5 );
+ retval[3] = std::conj(retval[2]);
+ }
+ else
+ {
+ sqD = sqrt(D);
+ retval[2].real( (-p2 + sqD) * 0.5 );
+ retval[3].real( (-p2 - sqD) * 0.5 );
+ }
+
+ return retval;
+}
+
diff --git a/Framework/Utilities/quartic.h b/Framework/Utilities/quartic.h
new file mode 100644
index 0000000000000000000000000000000000000000..968ebbbd79ca170d8a7ae2ec27737240203b57b5
--- /dev/null
+++ b/Framework/Utilities/quartic.h
@@ -0,0 +1,69 @@
+/***************************************************************************
+ * Copyright (C) 2016 by Саша Миленковић *
+ * sasa.milenkovic.xyz@gmail.com *
+ * *
+ * This program is free software; you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation; either version 2 of the License, or *
+ * (at your option) any later version. *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * ( http://www.gnu.org/licenses/gpl-3.0.en.html ) *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program; if not, write to the *
+ * Free Software Foundation, Inc., *
+ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
+ ***************************************************************************/
+
+#ifndef QUARTIC_H_INCLUDED
+#define QUARTIC_H_INCLUDED
+
+#include <complex>
+
+const double PI = 3.141592653589793238463L;
+const double M_2PI = 2*PI;
+const double eps=1e-12;
+
+typedef std::complex<double> DComplex;
+
+//---------------------------------------------------------------------------
+// useful for testing
+ inline DComplex polinom_2(DComplex x, double a, double b)
+ {
+ //Horner's scheme for x*x + a*x + b
+ return x * (x + a) + b;
+ }
+
+//---------------------------------------------------------------------------
+// useful for testing
+ inline DComplex polinom_3(DComplex x, double a, double b, double c)
+ {
+ //Horner's scheme for x*x*x + a*x*x + b*x + c;
+ return x * (x * (x + a) + b) + c;
+ }
+
+//---------------------------------------------------------------------------
+// useful for testing
+ inline DComplex polinom_4(DComplex x, double a, double b, double c, double d)
+ {
+ //Horner's scheme for x*x*x*x + a*x*x*x + b*x*x + c*x + d;
+ return x * (x * (x * (x + a) + b) + c) + d;
+ }
+
+//---------------------------------------------------------------------------
+// x - array of size 3
+// In case 3 real roots: => x[0], x[1], x[2], return 3
+// 2 real roots: x[0], x[1], return 2
+// 1 real root : x[0], x[1] ± i*x[2], return 1
+unsigned int solveP3(double* x, double a, double b, double c);
+
+//---------------------------------------------------------------------------
+// solve quartic equation x^4 + a*x^3 + b*x^2 + c*x + d
+// Attention - this function returns dynamically allocated array. It has to be released afterwards.
+DComplex* solve_quartic(double a, double b, double c, double d);
+
+
+#endif // QUARTIC_H_INCLUDED
diff --git a/Processes/TrackingLine/TrackingLine_interpolation.h b/Processes/TrackingLine/TrackingLine_interpolation.h
new file mode 100644
index 0000000000000000000000000000000000000000..228855cfe2c0ba93e69f0f1aebc1d890fa5f9442
--- /dev/null
+++ b/Processes/TrackingLine/TrackingLine_interpolation.h
@@ -0,0 +1,202 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_corsika_processes_TrackingLine_h_
+#define _include_corsika_processes_TrackingLine_h_
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/quartic.h>
+#include <cmath>
+#include <optional>
+#include <type_traits>
+#include <utility>
+
+namespace corsika::environment {
+ template <typename IEnvironmentModel>
+ class Environment;
+ template <typename IEnvironmentModel>
+ class VolumeTreeNode;
+} // namespace corsika::environment
+
+namespace corsika::process {
+
+ namespace tracking_line {
+
+ std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+ TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
+
+ corsika::units::si::TimeType TimeOfIntersection(geometry::Line const&,
+ geometry::Plane const&);
+
+ class TrackingLine {
+
+ public:
+ TrackingLine(){};
+
+ template <typename Particle> // was Stack previously, and argument was
+ // Stack::StackIterator
+ auto GetTrack(Particle const& p) {
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ geometry::Vector<SpeedType::dimension_type> velocity =
+ p.GetMomentum() / p.GetEnergy() * corsika::units::constants::c;
+
+ std::complex<double>* solutions = solve_quartic(1, 0, 1, -20);
+ std::vector<double> tmp;
+ for(int i = 0; i < 4; i++) {
+ if(solutions[i].imag() == 0 && solutions[i].real() >= 0) {
+ tmp.push_back(solutions[i].real());
+ }
+ }
+ double s = *std::min_element(tmp.begin(),tmp.end());
+ std::cout << "s = " << s << std::endl;
+ std::cout << "x1 = " << (solutions[0].real()>=0. ? " " : "") << solutions[0].real(); if(solutions[0].imag()!=0.0) std::cout << " + i * " << solutions[0].imag(); std::cout << std::endl;
+ std::cout << "x2 = " << (solutions[1].real()>=0. ? " " : "") << solutions[1].real(); if(solutions[1].imag()!=0.0) std::cout << " - i * " << -solutions[1].imag(); std::cout << std::endl;
+ std::cout << "x3 = " << (solutions[2].real()>=0. ? " " : "") << solutions[2].real(); if(solutions[2].imag()!=0.0) std::cout << " + i * " << solutions[2].imag(); std::cout << std::endl;
+ std::cout << "x4 = " << (solutions[3].real()>=0. ? " " : "") << solutions[3].real(); if(solutions[3].imag()!=0.0) std::cout << " - i * " << -solutions[3].imag(); std::cout << std::endl;
+
+ auto const currentPosition = p.GetPosition();
+ std::cout << "TrackingLine pid: " << p.GetPID()
+ << " , E = " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+ std::cout << "TrackingLine pos: "
+ << currentPosition.GetCoordinates()
+ // << " [" << p.GetNode()->GetModelProperties().GetName() << "]"
+ << std::endl;
+ std::cout << "TrackingLine E: " << p.GetEnergy() / 1_GeV << " GeV" << std::endl;
+
+ // determine velocity after adding magnetic field
+ auto const* currentLogicalVolumeNode = p.GetNode();
+ int chargeNumber;
+ if(corsika::particles::IsNucleus(p.GetPID())) {
+ chargeNumber = p.GetNuclearZ();
+ } else {
+ chargeNumber = corsika::particles::GetChargeNumber(p.GetPID());
+ }
+ geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
+ auto magMaxLength = 1_m/0;
+ auto directionAfter = directionBefore;
+ if(chargeNumber != 0) {
+ auto magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(currentPosition);
+ std::cout << "TrackingLine B: " << magneticfield.GetComponents() / 1_uT << " uT " << std::endl;
+ geometry::Vector<SpeedType::dimension_type> const velocityVerticalMag = velocity -
+ velocity.parallelProjectionOnto(magneticfield);
+ LengthType const gyroradius = p.GetEnergy() * velocityVerticalMag.GetNorm() * 1_V /
+ (corsika::units::constants::cSquared * abs(chargeNumber) *
+ magneticfield.GetNorm() * 1_eV);
+ //steplength should consider more things than just gyroradius
+ double maxAngle = 0.1;
+ LengthType const Steplength = 2 * sin(maxAngle * M_PI / 180) * gyroradius;
+ // First Movement
+ auto position = currentPosition + directionBefore * Steplength / 2;
+ // Change of direction by magnetic field at position
+ magneticfield = currentLogicalVolumeNode->GetModelProperties().GetMagneticField(position);
+ directionAfter = directionBefore + directionBefore.cross(magneticfield) * chargeNumber *
+ Steplength * corsika::units::constants::cSquared * 1_eV /
+ (p.GetEnergy() * velocity.GetNorm() * 1_V);
+ // Second Movement
+ position = position + directionAfter * Steplength / 2;
+ magMaxLength = (position - currentPosition).GetNorm();
+ geometry::Vector<dimensionless_d> const direction = (position - currentPosition) /
+ magMaxLength;
+ velocity = direction * velocity.GetNorm();
+ std::cout << "TrackingLine p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
+ << " GeV " << std::endl;
+
+ } else {
+ std::cout << "TrackingLine p: " << p.GetMomentum().GetComponents() / 1_GeV
+ << " GeV " << std::endl;
+ }
+
+ std::cout << "TrackingLine v: " << velocity.GetComponents() << std::endl;
+
+ geometry::Line line(currentPosition, velocity);
+
+ //auto const* currentLogicalVolumeNode = p.GetNode();
+ //~ auto const* currentNumericalVolumeNode =
+ //~ fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+ auto const numericallyInside =
+ currentLogicalVolumeNode->GetVolume().Contains(currentPosition);
+
+ std::cout << "numericallyInside = " << (numericallyInside ? "true" : "false");
+
+ auto const& children = currentLogicalVolumeNode->GetChildNodes();
+ auto const& excluded = currentLogicalVolumeNode->GetExcludedNodes();
+
+ std::vector<std::pair<TimeType, decltype(p.GetNode())>> intersections;
+
+ // for entering from outside
+ auto addIfIntersects = [&](auto const& vtn) {
+ auto const& volume = vtn.GetVolume();
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ volume); // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+
+ if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+ auto const [t1, t2] = *opt;
+ std::cout << "intersection times: " << t1 / 1_s << "; "
+ << t2 / 1_s
+ // << " " << vtn.GetModelProperties().GetName()
+ << std::endl;
+ if (t1.magnitude() > 0)
+ intersections.emplace_back(t1, &vtn);
+ else if (t2.magnitude() > 0)
+ std::cout << "inside other volume" << std::endl;
+ }
+ };
+
+ for (auto const& child : children) { addIfIntersects(*child); }
+ for (auto const* ex : excluded) { addIfIntersects(*ex); }
+
+ {
+ auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+ currentLogicalVolumeNode->GetVolume());
+ // for the moment we are a bit bold here and assume
+ // everything is a sphere, crashes with exception if not
+ [[maybe_unused]] auto const [t1, t2] = *TimeOfIntersection(line, sphere);
+ [[maybe_unused]] auto dummy_t1 = t1;
+ intersections.emplace_back(t2, currentLogicalVolumeNode->GetParent());
+ }
+
+ auto const minIter = std::min_element(
+ intersections.cbegin(), intersections.cend(),
+ [](auto const& a, auto const& b) { return a.first < b.first; });
+
+ TimeType min;
+
+ if (minIter == intersections.cend()) {
+ min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+ // to handle the numerics properly
+ throw std::runtime_error("no intersection with anything!");
+ } else {
+ min = minIter->first;
+ }
+
+ std::cout << " t-intersect: "
+ << min
+ // << " " << minIter->second->GetModelProperties().GetName()
+ << std::endl;
+
+ return std::make_tuple(geometry::Trajectory<geometry::Line>(line, min),
+ velocity.norm() * min, minIter->second, magMaxLength,
+ directionBefore, directionAfter);
+ }
+ };
+
+ } // namespace tracking_line
+
+} // namespace corsika::process
+
+#endif