initial attempt to move from delta-L to L

corsika/detail/framework/geometry/LeapFrogTrajectory.inl

@@ -12,6 +12,7 @@
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
 
 namespace corsika {
 
@@ -45,7 +46,11 @@ namespace corsika {
   }
 
   inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
-    return u * timeStep_;
+    TimeType const step = timeStep_ * u;
+    double const correction = 1;
+    //    (initialDirection_ + initialDirection_.cross(magneticfield_) * step *
+    //    k_).getNorm();
+    return step / 2 * (correction + 1);
   }
 
   inline LengthType LeapFrogTrajectory::getLength(double const u) const {
@@ -58,9 +63,20 @@ namespace corsika {
   }
 
   inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
-    // double const correction =
-    //    (1. + fabs(0.5 * initialDirection_.cross(magneticfield_).getNorm() * limit *
-    //    k_));
+    /*
+    initial attempt to calculate delta-L from assumed full-leap-frog-length L:
+
+  Note: often return 0. Not good enough yet.
+
+    LengthType const L = initialVelocity_.getNorm() * limit; // distance
+    double const a = (initialVelocity_.cross(magneticfield_) * k_).getSquaredNorm() / 4 /
+                     square(1_m) * static_pow<4>(1_s);
+    double const d = L * initialVelocity_.getNorm() / square(1_m) * 1_s;
+    double const e = -square(L) / square(1_m);
+    std::vector<double> solutions = solve_quartic_real(a, 0, 0, d, e);
+    CORSIKA_LOG_DEBUG("setDuration limit={} L={} solution={}", limit, L,
+                      fmt::join(solutions, ", "));
+    */
     timeStep_ = limit;
   }