From c772f06c1bae9a4699777cde26f661abcc6a172c Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Fri, 11 Jun 2021 10:40:47 +0200
Subject: [PATCH] initial attempt to move from delta-L to L
---
.../framework/geometry/LeapFrogTrajectory.inl | 24 +++++++++++++++----
1 file changed, 20 insertions(+), 4 deletions(-)
diff --git a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
index eed4221ce..7293a7f59 100644
--- a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl
+++ b/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;
}
--
GitLab