mal 1.0001

Framework/Cascade/Cascade.h

@@ -227,12 +227,11 @@ namespace corsika::cascade {
       std::cout << "distance_max=" << distance_max << std::endl;
 
       // take minimum of geometry, interaction, decay for next step
-      std::cout << "Interaction: " << distance_interact << std::endl;
-      std::cout << "Decay: " << distance_decay << std::endl;
-      std::cout << "ObsPlane: " << distance_max << std::endl;
-      std::cout << "Transition: " << geomMaxLength << std::endl;
-      auto const min_distance = std::min(
+      auto min_distance = std::min(
           {distance_interact, distance_decay, distance_max, geomMaxLength});
+      if (min_distance == geomMaxLength) {
+        min_distance = 1.001 * geomMaxLength;
+      } 
 
       C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
 
@@ -260,11 +259,11 @@ namespace corsika::cascade {
   		// Second Movement
   		position = position + directionAfter * min_distance / 2;
       auto distance = position - vParticle.GetPosition();
-      //distance.norm() != min_distance if q != 0
-      //small error can be neglected
+      // distance.norm() != min_distance if q != 0
+      // small error can be neglected
       if (distance.norm() != 0_m) {
       velocity = distance.normalized() * velocity.norm();
-      } //no velocity update for very small steps
+      } // no velocity update for very small steps
       
       // here the particle is actually moved along the trajectory to new position:
       // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);

Processes/ObservationPlane/ObservationPlane.cc

@@ -88,6 +88,7 @@ LengthType ObservationPlane::MaxStepLength(setup::Stack::ParticleType const& vPa
       plane_.GetNormal().dot(velocity.cross(magneticfield)) * 2 * k)) - 
       velocity.dot(plane_.GetNormal()) / velocity.GetNorm() ) / 
       (plane_.GetNormal().dot(velocity.cross(magneticfield)) * k);
+      
     if (MaxStepLength1 <= 0_m && MaxStepLength2 <= 0_m) {
       return std::numeric_limits<double>::infinity() * 1_m;
     } else if (MaxStepLength1 <= 0_m || MaxStepLength2 < MaxStepLength1) {

Processes/TrackingLine/TrackingLine.h

@@ -135,7 +135,7 @@ namespace corsika::process {
   									                                            (position - currentPosition).GetNorm();
   		      velocity1 = direction * velocity.GetNorm();
           } // instead of changing the line with magnetic field, the TimeOfIntersection() could be changed
-          // line has some errors
+          // using line has some errors for huge steps
           geometry::Line line(currentPosition, velocity1);
 
           if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
@@ -232,9 +232,10 @@ namespace corsika::process {
                                                                   min * k;
         // Second Movement
         position = position + directionAfter * velocity.norm() * min / 2;
-        geometry::Vector<dimensionless_d> const direction = (position - currentPosition) / 
-                                                            (position - currentPosition).GetNorm();
-        velocity = direction * velocity.GetNorm();
+        if ((position - currentPosition).GetNorm() != 0_m) {
+          geometry::Vector<dimensionless_d> const direction = (position - currentPosition).normalized();
+          velocity = direction * velocity.norm();
+        } // no velocity update for very small steps
         geometry::Line lineWithB(currentPosition, velocity);
 
         return std::make_tuple(geometry::Trajectory<geometry::Line>(lineWithoutB, min),