Test

Framework/Cascade/Cascade_interpolation.h

+/*
+ * (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

Framework/Utilities/quartic.cpp

+/***************************************************************************
+ *   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;
+}
+

Framework/Utilities/quartic.h

+/***************************************************************************
+ *   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

Processes/TrackingLine/TrackingLine_interpolation.h

+/*
+ * (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