IAP GITLAB

Commit 8bd31c91 authored by Matthieu Carrere's avatar Matthieu Carrere

added new commits of master branch

parents 5b5404fa e623e3ed
Pipeline #5017 passed with stages
in 26 minutes and 19 seconds
......@@ -106,11 +106,14 @@ include (conan) # self-provided in 'cmake' directory
#
# download and build all dependencies
message (STATUS "Installing dependencies from Conan (this may take some time)...")
if(NOT APPLE)
set(SETTINGS compiler.libcxx=libstdc++11) # not available on MacOS
endif()
conan_cmake_run (CONANFILE conanfile.txt
BASIC_SETUP CMAKE_TARGETS
BUILD missing
BUILD_TYPE "Release"
SETTINGS compiler.libcxx=libstdc++11)
SETTINGS ${SETTINGS})
#
# add cnpy temporarily. As long as SaveBoostHistogram does not save to parquet directly
#
......@@ -237,6 +240,7 @@ add_subdirectory (modules/sibyll)
add_subdirectory (modules/qgsjetII)
add_subdirectory (modules/urqmd)
add_subdirectory (modules/conex)
add_subdirectory (modules/epos)
#+++++++++++++++++++++++++++++++
# unit testing
......
......@@ -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 {
......@@ -35,27 +36,49 @@ namespace corsika {
}
inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
return initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
return getDirection(u) * initialVelocity_.getNorm();
}
inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
return getVelocity(u).normalized();
return (initialDirection_ +
initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
.normalized();
}
inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
return u * timeStep_ *
(double(getVelocity(u).getNorm() / initialVelocity_.getNorm()) + 1.0) / 2;
TimeType const step = timeStep_ * u;
double const correction = 1;
// the eventual (delta-L to L) correction factor is:
// (initialDirection_ + initialDirection_.cross(magneticfield_) * step *
// k_).getNorm();
return step / 2 * (correction + 1);
}
inline LengthType LeapFrogTrajectory::getLength(double const u) const {
return timeStep_ * initialVelocity_.getNorm() * u;
return getDuration(u) * initialVelocity_.getNorm();
}
inline void LeapFrogTrajectory::setLength(LengthType const limit) {
if (initialVelocity_.getNorm() == 0_m / 1_s) setDuration(0_s);
if (initialVelocity_.getNorm() == SpeedType::zero()) setDuration(0_s);
setDuration(limit / initialVelocity_.getNorm());
}
inline void LeapFrogTrajectory::setDuration(TimeType const limit) { timeStep_ = limit; }
inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
/*
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;
}
} // namespace corsika
......@@ -18,6 +18,8 @@ namespace corsika {
//---------------------------------------------------------------------------
// solve cubic equation A x^3 + B*x^2 + C*x + D = 0
// x^3 + a*x^2 + b*x + c = 0
// mainly along WolframAlpha formulas
inline std::vector<double> solve_cubic_real_analytic(long double A, long double B,
long double C, long double D,
double const epsilon) {
......@@ -29,8 +31,8 @@ namespace corsika {
long double c = D / A;
long double a2 = a * a;
long double q = (a2 - 3 * b) / 9;
long double r = (a * (2 * a2 - 9 * b) + 27 * c) / 54;
long double q = (3 * b - a2) / 9;
long double r = (a * (9 * b - 2 * a2) - 27 * c) / 54;
long double q3 = q * q * q;
// disc = q**3 + r**2
......@@ -38,8 +40,8 @@ namespace corsika {
long double w = r * r;
long double e = std::fma(r, r, -w);
// s:t = q*q exactly
long double s = -q * q;
long double t = std::fma(-q, q, -s);
long double s = q * q;
long double t = std::fma(q, q, -s);
// s:t * q + w:e = s*q + w + t*q +e = s*q+w + u:v + e = f + u:v + e
long double f = std::fma(s, q, w);
long double u = t * q;
......@@ -51,30 +53,34 @@ namespace corsika {
au = f - au;
ab = u - ab;
// sum all terms into final result
long double const disc = -(((e + uf) + au) + ab) + v;
long double const disc = (((e + uf) + au) + ab) + v;
if (disc >= 0) {
long double t = r / std::sqrt(q3);
CORSIKA_LOG_TRACE("disc={} {}", disc, q3 + r * r);
if (std::abs(disc) < epsilon) {
a /= 3;
long double const cbrtR = std::cbrt(r);
return {double(2 * cbrtR - a), double(-cbrtR - a)}; // 2nd solution is doublet
} else if (disc > 0) {
long double const S = std::cbrt(r + std::sqrt(disc));
long double const T = std::cbrt(r - std::sqrt(disc));
a /= 3;
return {double((S + T) - a)}; // plus two imaginary solution
} else { // disc < 0
long double t = r / std::sqrt(-q3);
if (t < -1) t = -1;
if (t > 1) t = 1;
t = std::acos(t);
a /= 3;
q = -2 * std::sqrt(q);
q = 2 * std::sqrt(-q);
return {double(q * std::cos(t / 3) - a),
double(q * std::cos((t + 2 * M_PI) / 3) - a),
double(q * std::cos((t - 2 * M_PI) / 3) - a)};
} else {
long double term1 = -cbrt(std::fabs(r) + std::sqrt(-disc));
if (r < 0) term1 = -term1;
long double term2 = (0 == term1 ? 0 : q / term1);
a /= 3;
long double test = 0.5 * std::sqrt(3.) * (term1 - term2);
if (std::fabs(test) < epsilon) {
return {double((term1 + term2) - 1), double(-0.5 * (term1 + term2) - a)};
}
return {double((term1 + term2) - a)};
double(q * std::cos((t + 4 * M_PI) / 3) - a)};
}
}
} // namespace andre
......@@ -228,67 +234,95 @@ namespace corsika {
}
/**
* Iterative approach.
* Iterative approach. https://en.wikipedia.org/wiki/Halley%27s_method
* Halley's method
*/
inline std::vector<double> solve_cubic_real(long double a, long double b, long double c,
long double d, double const epsilon) {
CORSIKA_LOG_TRACE("cubic_2: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}", a, b,
c, d, epsilon, (std::abs(a - 1) < epsilon),
CORSIKA_LOG_TRACE("cubic_iterative: a={:f}, b={:f}, c={:f}, d={:f}, epsilon={} {} {}",
a, b, c, d, epsilon, (std::abs(a - 1) < epsilon),
(std::abs(b) < epsilon));
if (std::abs(a) < epsilon) { // this is just a quadratic
return solve_quadratic_real(b, c, d, epsilon);
}
long double const dist = std::fma(b / a, b / a, -3 * c / a);
long double const xinfl = -b / (a * 3);
long double x1 = xinfl;
long double f_x1 = cubic_function(xinfl, a, b, c, d);
auto pre_opt = andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
long double x1 = 0; // start value
if (std::abs(f_x1) > epsilon) {
if (std::abs(dist) < epsilon) {
x1 = xinfl - std::cbrt(f_x1);
} else if (dist > 0) {
if (f_x1 > 0)
x1 = xinfl - 2 / 3 * std::sqrt(dist);
else
x1 = xinfl + 2 / 3 * std::sqrt(dist);
if (pre_opt.size()) {
x1 = pre_opt[0]; //*std::max_element(pre_opt.begin(), pre_opt.end());
#ifdef DEBUG
for (long double test_v : pre_opt) {
CORSIKA_LOG_TRACE("test,andre x={} f(x)={}", test_v,
cubic_function(test_v, a, b, c, d));
}
int niter = 0;
const int maxiter = 100;
do {
long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
// if (potential) saddle point... avoid
if (std::abs(f_prime_x1) < epsilon) {
x1 -= std::cbrt(f_x1);
} else {
x1 -=
f_x1 * f_prime_x1 / (static_pow<2>(f_prime_x1) - 0.5 * f_x1 * f_prime2_x1);
#endif
} else {
// this is only if the former solve_cubic_real_analytic would not result
// in any solution. We have no test case for this. This is excluded from tests:
// LCOV_EXCL_START
long double const dist = std::fma(b / a, b / a, -3 * c / a);
long double const xinfl = -b / (a * 3);
x1 = xinfl;
long double f_test = cubic_function(xinfl, a, b, c, d);
if (std::abs(f_test) > epsilon) {
if (std::abs(dist) < epsilon) {
x1 = xinfl - std::cbrt(f_test);
} else if (dist > 0) {
if (f_test > 0)
x1 = xinfl - 2 / 3 * std::sqrt(dist);
else
x1 = xinfl + 2 / 3 * std::sqrt(dist);
}
f_x1 = cubic_function(x1, a, b, c, d);
CORSIKA_LOG_TRACE("niter={} x1={} f_x1={} f_prime={} f_prime2={} eps={}", niter,
x1, f_x1, f_prime_x1, f_prime2_x1, epsilon);
} while ((++niter < maxiter) && (std::abs(f_x1) > epsilon));
}
// LCOV_EXCL_STOP
}
CORSIKA_LOG_TRACE("niter={}", niter);
long double f_x1 = cubic_function(x1, a, b, c, d);
long double dx1 = 0;
int niter = 0;
const int maxiter = 100;
do {
long double const f_prime_x1 = cubic_function_dfdx(x1, a, b, c);
long double const f_prime2_x1 = cubic_function_d2fd2x(x1, a, b);
// if (potential) saddle point... avoid
if (std::abs(f_prime_x1) < epsilon) {
dx1 = std::cbrt(f_x1);
} else {
dx1 = f_x1 * f_prime_x1 * 2 / (f_prime_x1 * f_prime_x1 * 2 - f_x1 * f_prime2_x1);
}
x1 -= dx1;
f_x1 = cubic_function(x1, a, b, c, d);
CORSIKA_LOG_TRACE(
"niter={} x1={:.20f} f_x1={:.20f} f_prime={:.20f} f_prime2={:.20f} dx1={}",
niter, x1, f_x1, f_prime_x1, f_prime2_x1,
f_x1 * f_prime_x1 / (f_prime_x1 * f_prime_x1 - f_x1 * f_prime2_x1 * 0.5));
} while ((++niter < maxiter) && (std::abs(f_x1) > epsilon * 1000) &&
(std::abs(dx1) > epsilon));
CORSIKA_LOG_TRACE("niter={}", niter);
if (niter >= maxiter) {
CORSIKA_LOG_DEBUG("niter reached max iterations {}", niter);
return andre::solve_cubic_real_analytic(a, b, c, d, epsilon);
}
CORSIKA_LOG_TRACE("x1={} f_x1={}", x1, f_x1);
double const b1 = x1 + b / a;
double const b0 = b1 * x1 + c / a;
std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, epsilon);
std::vector<double> quad_check = solve_quadratic_real(1, b1, b0, 1e-3);
CORSIKA_LOG_TRACE("quad_check=[{}], f(z)={}", fmt::join(quad_check, ", "),
cubic_function(x1, a, b, c, d));
quad_check.push_back(x1);
CORSIKA_LOG_TRACE("cubic: solve_cubic_real returns={}", fmt::join(quad_check, ", "));
return quad_check;
}
} // namespace corsika
} // namespace corsika
......@@ -62,7 +62,7 @@ namespace corsika {
long double f = std::fma(b, b, -w);
long double radicant = f + e;
CORSIKA_LOG_TRACE(" radicant={} {} ", radicant, b * b - a * c * 4);
CORSIKA_LOG_TRACE("radicant={} {} ", radicant, b * b - a * c * 4);
if (std::abs(radicant) < epsilon) { // just one real solution
return {double(-b / (2 * a))};
......
......@@ -35,6 +35,9 @@ namespace corsika {
// y^3 − c*y^2 + (bd−4e)*y − b^2*e−d^2+4*c*e = 0
std::vector<double> x3 = solve_cubic_real(1, a3, b3, c3, epsilon);
if (!x3.size()) {
return {}; // no solution, numeric problem (LCOV_EXCL_LINE)
}
long double y = x3[0]; // there is always at least one solution
// The essence - choosing Y with maximal absolute value.
if (x3.size() == 3) {
......@@ -74,8 +77,8 @@ namespace corsika {
// solving quadratic eqs. x^2 + p1*x + q1 = 0
// x^2 + p2*x + q2 = 0
std::vector<double> quad1 = solve_quadratic_real(1, p1, q1);
std::vector<double> quad2 = solve_quadratic_real(1, p2, q2);
std::vector<double> quad1 = solve_quadratic_real(1, p1, q1, 1e-5);
std::vector<double> quad2 = solve_quadratic_real(1, p2, q2, 1e-5);
if (quad2.size() > 0) {
for (auto val : quad2) quad1.push_back(val);
}
......@@ -109,17 +112,37 @@ namespace corsika {
}
CORSIKA_LOG_TRACE("check m={}", m);
if (m == 0) { return {0}; }
CORSIKA_LOG_TRACE("check m={}", m);
if (m < 0) {
// this is a rare numerical instability
// first: try analytical solution, second: discard (curved->straight tracking)
std::vector<double> const resolve_cubic_analytic =
andre::solve_cubic_real_analytic(1, p, p2 / 4 - r, -q2 / 8, epsilon);
CORSIKA_LOG_TRACE("andre::resolve_cubic_analytic: N={}, m=[{}]",
resolve_cubic_analytic.size(),
fmt::join(resolve_cubic_analytic, ", "));
if (!resolve_cubic_analytic.size()) return {};
for (auto const& v : resolve_cubic_analytic) {
CORSIKA_LOG_TRACE("check pol3(v)={}", (static_pow<3>(v) + static_pow<2>(v) * p +
v * (p2 / 4 - r) - q2 / 8));
if (std::abs(v) > epsilon && std::abs(v) > m) { m = v; }
}
CORSIKA_LOG_TRACE("check m={}", m);
if (m == 0) { return {0}; }
if (m < 0) {
return {}; // now we are out of options, cannot solve: curved->straight tracking
}
}
long double const quad_term1 = p / 2 + m;
long double const quad_term2 = std::sqrt(2 * m);
long double const quad_term3 = q / (2 * quad_term2);
std::vector<double> z_quad1 =
solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, epsilon);
solve_quadratic_real(1, quad_term2, quad_term1 - quad_term3, 1e-5);
std::vector<double> z_quad2 =
solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, epsilon);
solve_quadratic_real(1, -quad_term2, quad_term1 + quad_term3, 1e-5);
for (auto const& z : z_quad2) z_quad1.push_back(z);
return z_quad1;
}
......
......@@ -27,16 +27,22 @@ namespace corsika {
The current step did not yet reach the ObservationPlane, do nothing now and wait:
*/
if (!stepLimit) {
#ifdef DEBUG
// @todo this is actually needed to fix small instabilities of the leap-frog
// tracking: Note, this is NOT a general solution and should be clearly revised with
// a more robust tracking. #ifdef DEBUG
if (deleteOnHit_) {
// since this is basically a bug, it cannot be tested LCOV_EXCL_START
LengthType const check =
(particle.getPosition() - plane_.getCenter()).dot(plane_.getNormal());
if (check < 0_m) {
CORSIKA_LOG_DEBUG("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
}
}
#endif
return ProcessReturn::Ok;
CORSIKA_LOG_WARN("PARTICLE AVOIDED OBSERVATIONPLANE {}", check);
CORSIKA_LOG_WARN("Temporary fix: write and remove particle.");
} else
return ProcessReturn::Ok;
// LCOV_EXCL_STOP
} else
// #endif
return ProcessReturn::Ok;
}
HEPEnergyType const energy = particle.getEnergy();
......@@ -63,8 +69,9 @@ namespace corsika {
inline LengthType ObservationPlane<TTracking, TOutput>::getMaxStepLength(
TParticle const& particle, TTrajectory const& trajectory) {
CORSIKA_LOG_TRACE("particle={}, pos={}, dir={}, plane={}", particle.asString(),
particle.getPosition(), particle.getDirection(), plane_.asString());
CORSIKA_LOG_TRACE("getMaxStepLength, particle={}, pos={}, dir={}, plane={}",
particle.asString(), particle.getPosition(),
particle.getDirection(), plane_.asString());
auto const intersection = TTracking::intersect(particle, plane_);
......@@ -77,10 +84,10 @@ namespace corsika {
return std::numeric_limits<double>::infinity() * 1_m;
}
double const fractionOfIntersection = timeOfIntersection / trajectory.getDuration();
auto const pointOfIntersection = trajectory.getPosition(fractionOfIntersection);
auto dist = (trajectory.getPosition(0) - pointOfIntersection).getNorm();
CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength l={} m", dist / 1_m);
return dist;
CORSIKA_LOG_TRACE("ObservationPlane: getMaxStepLength dist={} m, pos={}",
trajectory.getLength(fractionOfIntersection) / 1_m,
trajectory.getPosition(fractionOfIntersection));
return trajectory.getLength(fractionOfIntersection);
}
template <typename TTracking, typename TOutput>
......
......@@ -73,7 +73,7 @@ namespace corsika {
// Sternheimer parameterization, density corrections towards high energies
// NOTE/TODO: when Cbar is 0 it needs to be approximated from parameterization ->
// MISSING
CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m{}=",
CORSIKA_LOG_TRACE("BetheBloch p.getMomentum().getNorm()/m={}",
p.getMomentum().getNorm() / m);
double const x = log10(p.getMomentum().getNorm() / m);
double delta = 0;
......
This diff is collapsed.
/*
* (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
*
* 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.
*/
#pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <epos.hpp>
namespace corsika::epos {
inline HEPMassType getEposMass(Code const pCode) {
if (pCode == Code::Nucleus)
throw std::runtime_error("Cannot getMass() of particle::Nucleus -> unspecified");
auto sCode = convertToEposRaw(pCode);
if (sCode == 0)
throw std::runtime_error("getEposMass: unknown particle!");
else {
float mass2 = 0;
::epos::idmass_(sCode, mass2);
return sqrt(mass2) * 1_GeV;
}
}
inline PDGCode getEposPDGId(Code const p) {
if (!is_nucleus(p)) {
int eid = corsika::epos::convertToEposRaw(p);
char nxs[4] = "nxs";
char pdg[4] = "pdg";
return static_cast<PDGCode>(::epos::idtrafo_(nxs, pdg, eid));
} else {
throw std::runtime_error("Epos id conversion not implemented for nuclei!");
}
}
} // namespace corsika::epos
......@@ -73,14 +73,18 @@ namespace corsika::proposal {
for (auto& it : rnd) it = distr(RNG_);
// calculate deflection based on particle energy, loss
[[maybe_unused]] auto deflection = (c->second).scatter->CalculateMultipleScattering(
auto deflection = (c->second).scatter->CalculateMultipleScattering(
grammage / 1_g * square(1_cm), vP.getEnergy() / 1_MeV, E_f / 1_MeV, rnd);
// TODO: multiple scattering is temporary deactivated !!!!!
[[maybe_unused]] auto [unused1, final_direction] =
PROPOSAL::multiple_scattering::ScatterInitialDirection(direction, deflection);
// update particle direction after continuous loss caused by multiple
// scattering
auto vec = QuantityVector(direction.GetX() * E_f, direction.GetY() * E_f,
direction.GetZ() * E_f);
auto particle_momentum = vP.getMomentum().getNorm();
auto vec = QuantityVector(final_direction.GetX() * particle_momentum,
final_direction.GetY() * particle_momentum,
final_direction.GetZ() * particle_momentum);
vP.setMomentum(MomentumVector(vP_dir.getCoordinateSystem(), vec));
}
......@@ -108,7 +112,9 @@ namespace corsika::proposal {
// if the particle has a charge take multiple scattering into account
if (vP.getChargeNumber() != 0) scatter(vP, dE, dX);
vP.setEnergy(final_energy);
vP.setMomentum(vP.getMomentum() * vP.getEnergy() / vP.getMomentum().getNorm());
auto new_momentum =
sqrt(vP.getEnergy() * vP.getEnergy() - vP.getMass() * vP.getMass());
vP.setMomentum(vP.getMomentum() * new_momentum / vP.getMomentum().getNorm());
return ProcessReturn::Ok;
}
......
......@@ -174,7 +174,7 @@ namespace corsika::pythia8 {
auto const t0 = projectile.getTime();
auto const& labMomentum = projectile.getMomentum();
CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
[[maybe_unused]] CoordinateSystemPtr const& labCS = labMomentum.getCoordinateSystem();
// define target kinematics in lab frame
// define boost to and from CoM frame
......
......@@ -44,9 +44,10 @@ namespace corsika {
fmt::ptr(&volumeNode), fmt::ptr(volumeNode.getParent()),
time_intersections_curr.hasIntersections());
if (time_intersections_curr.hasIntersections()) {
CORSIKA_LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
time_intersections_curr.getEntry() / 1_s,
time_intersections_curr.getExit() / 1_s);
CORSIKA_LOG_DEBUG(
"intersection times with currentLogicalVolumeNode: entry={} s and exit={} s",
time_intersections_curr.getEntry() / 1_s,
time_intersections_curr.getExit() / 1_s);
if (time_intersections_curr.getExit() <= minTime) {
minTime =
time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
......
......@@ -29,10 +29,10 @@ namespace corsika {
namespace tracking_leapfrog_curved {
template <typename TParticle>
inline auto make_LeapFrogStep(TParticle const& particle, LengthType steplength) {
inline auto Tracking::makeStep(TParticle const& particle, LengthType steplength) {
if (particle.getMomentum().getNorm() == 0_GeV) {
return std::make_tuple(particle.getPosition(), particle.getMomentum() / 1_GeV,
double(0));
double(0) * 1_m);
} // charge of the particle
ElectricChargeType const charge = particle.getCharge();
auto const* currentLogicalVolumeNode = particle.getNode();
......@@ -40,9 +40,14 @@ namespace corsika {
currentLogicalVolumeNode->getModelProperties().getMagneticField(
particle.getPosition());
VelocityVector velocity = particle.getVelocity();
decltype(meter / (second * volt)) k =
charge * constants::cSquared * 1_eV /
(velocity.getNorm() * particle.getEnergy() * 1_V);
auto const p_norm =
constants::c * convert_HEP_to_SI<MassType::dimension_type>(
particle.getMomentum().getNorm()); // kg *m /s
// k = q/|p|
decltype(1 / (tesla * meter)) const k =
charge / p_norm; // * initialVelocity.getNorm();
DirectionVector direction = velocity.normalized();
auto position = particle.getPosition(); // First Movement
// assuming magnetic field does not change during movement
......@@ -94,6 +99,8 @@ namespace corsika {
particle.getMomentum().getParallelProjectionOnto(magneticfield))
.getNorm();
CORSIKA_LOG_TRACE("p_perp={} eV", p_perp / 1_eV);
if (p_perp < 1_eV) {
// particle travel along, parallel to magnetic field. Rg is
// "0", but for purpose of step limit we return infinity here.
......@@ -101,10 +108,20 @@ namespace corsika {
return getLinearTrajectory(particle);
}
LengthType const gyroradius = (convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
constants::c / (abs(charge) * magnitudeB));
LengthType const gyroradius = convert_HEP_to_SI<MassType::dimension_type>(p_perp) *
constants::c / (abs(charge) * magnitudeB);
double const maxRadians = 0.01;
if (gyroradius > 1e9_m) {
// this cannot be really unit-tested. It is hidden. LCOV_EXCL_START
CORSIKA_LOG_WARN(
"CurvedLeapFrog is not very stable for extremely high gyroradius steps. "
"Rg={} -> straight tracking.",
gyroradius);
return getLinearTrajectory(particle);
// LCOV_EXCL_STOP
}
double const maxRadians = 0.01; // maximally allowed deflection
LengthType const steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius;
TimeType const steplimit_time = steplimit / initialVelocity.getNorm();
CORSIKA_LOG_DEBUG("gyroradius {}, steplimit: {} = {}", gyroradius, steplimit,
......@@ -121,19 +138,20 @@ namespace corsika {
decltype(1 / (tesla * second)) const k =
charge / p_norm *
initialVelocity.getNorm(); // since we use steps in time and not length
// units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1/(T s)
// units: C * s / m / kg * m/s = 1 / (T*m) * m/s = 1 / (T*s)
return std::make_tuple(
LeapFrogTrajectory(position, initialVelocity, magneticfield, k,
minTime), // trajectory
minNode); // next volume node
minTime), // --> trajectory
minNode); // --> next volume node
}