IAP GITLAB

Commit c3c55ba7 authored by Andre Schmidt's avatar Andre Schmidt

mal 1.0001

parent 0c125437
Pipeline #2024 failed with stages
......@@ -209,12 +209,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;
}
std::cout << " move particle by : " << min_distance << std::endl;
......@@ -242,11 +241,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);
......
......@@ -82,6 +82,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) {
......
......@@ -138,7 +138,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()) {
......@@ -238,9 +238,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),
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment