Added Intersection

Documentation/Examples/vertical_EAS.cc

@@ -97,7 +97,7 @@ int main(int argc, char** argv) {
   node1->SetModelProperties<UniformMagneticField<HomogeneousMedium<IEnvModel>>>(
       B, 1_g / (1_cm * 1_cm * 1e9_cm), composition);
 
-  auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
+  /*auto node2 = std::make_unique<VolumeTreeNode<IEnvModel>>(
       std::make_unique<geometry::Sphere>(center, seaLevel + 100.0_km));
   node2->SetModelProperties<UniformMagneticField<SlidingPlanarExponential<IEnvModel>>>(
       B, center, 540.1778_g / (1_cm * 1_cm) / 772170.16_cm, -772170.16_cm, composition,
@@ -124,7 +124,7 @@ int main(int argc, char** argv) {
   node4->AddChild(std::move(node5));
   node3->AddChild(std::move(node4));
   node2->AddChild(std::move(node3));
-  node1->AddChild(std::move(node2));
+  node1->AddChild(std::move(node2));*/
   env.GetUniverse()->AddChild(std::move(node1));
 
   // setup particle stack, and add primary particle

Framework/Cascade/Cascade.h

@@ -229,26 +229,24 @@ namespace corsika::cascade {
     		geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
     		auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / 
                 (velocity.GetNorm() * vParticle.GetEnergy() * 1_V);
-    		LengthType const steplength = min_distance / 
-                                      (directionBefore + directionBefore.cross(magneticfield) * k / 2).GetNorm();
     		// First Movement
     		//assuming magnetic field does not change during movement
-    		auto position = vParticle.GetPosition() + directionBefore * Steplength / 2;
+    		auto position = vParticle.GetPosition() + directionBefore * min_distance / 2;
     		// Change of direction by magnetic field
     		geometry::Vector<dimensionless_d> const directionAfter = directionBefore + directionBefore.cross(magneticfield) *
-                                                                Steplength * k; 
+                                                                min_distance * k; 
     		// Second Movement
-    		position = position + directionAfter * Steplength / 2;
-    		geometry::Vector<dimensionless_d> const direction = (position - vParticle.GetPosition()) / 
-    									                                      (position - vParticle.GetPosition()).GetNorm();
-        vParticle.SetMomentum(direction * vParticle.GetMomentum().GetNorm());
+    		position = position + directionAfter * min_distance / 2;
+        // here the particle is actually moved along the trajectory to new position:
+        // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
+        vParticle.SetMomentum(directionAfter.normalized() * vParticle.GetMomentum().GetNorm());
+        vParticle.SetPosition(position);
+      } else {
+        vParticle.SetPosition(step.PositionFromArclength(min_distance));
       }
-
-      // 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.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
+      std::cout << "New Position: " << vParticle.GetPosition().GetCoordinates() << std::endl;
 
       step.LimitEndTo(min_distance);
 

Processes/TrackingLine/TrackingLine_interpolation.h

@@ -115,6 +115,11 @@ namespace corsika::process {
 		std::cout << "TrackingLine   p: " << (direction * p.GetMomentum().GetNorm()).GetComponents() / 1_GeV
                   << " GeV " << std::endl;
 
+            auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (velocity.GetNorm() * p.GetEnergy() * 1_V);
+        geometry::Vector<dimensionless_d> const directionBefore = velocity.normalized();
+            double test =((directionBefore.cross(magneticfield)).dot(position-currentPosition) * k + 1) / (1_m * 1_m * (directionBefore.cross(magneticfield)).GetSquaredNorm() * k * k);
+            std::cout << "Test: " << test << k << std::endl;
+            
         } else {
         	std::cout << "TrackingLine   p: " << p.GetMomentum().GetComponents() / 1_GeV
                   << " GeV " << std::endl;
@@ -143,6 +148,8 @@ namespace corsika::process {
           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
+          
+          std::cout << "Mittelpunkt: " << sphere.GetCenter().GetCoordinates() << std::endl;
 
           if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
             auto const [t1, t2] = *opt;