diff --git a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl index af31f8209eedad47bd9de9cbf1344c17c9ef758e..2c5390de3293e1e2560c9811c5c0062123ef2a54 100644 --- a/corsika/detail/framework/geometry/LeapFrogTrajectory.inl +++ b/corsika/detail/framework/geometry/LeapFrogTrajectory.inl @@ -54,6 +54,12 @@ namespace corsika { return step / 2 * (correction + 1); } + template <typename Particle> + inline TimeType LeapFrogTrajectory::getTime(Particle const& particle, + double const u) const { + return particle.getTime() + getDuration(u); + } + inline LengthType LeapFrogTrajectory::getLength(double const u) const { return getDuration(u) * initialVelocity_.getNorm(); } diff --git a/corsika/detail/framework/geometry/StraightTrajectory.inl b/corsika/detail/framework/geometry/StraightTrajectory.inl index a8fef9b035a38745095976a17f6a0644b720a2a5..38add4e50257e957847beff81c3f54c448d326a7 100644 --- a/corsika/detail/framework/geometry/StraightTrajectory.inl +++ b/corsika/detail/framework/geometry/StraightTrajectory.inl @@ -23,6 +23,12 @@ namespace corsika { return u * timeStep_; } + template <typename Particle> + inline TimeType StraightTrajectory::getTime(Particle const& particle, + double const u) const { + return particle.getTime() + getDuration(u); // timeStep_ * u; + } + inline LengthType StraightTrajectory::getLength(double const u) const { if (timeLength_ == 0_s) return 0_m; if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s) diff --git a/corsika/framework/geometry/LeapFrogTrajectory.hpp b/corsika/framework/geometry/LeapFrogTrajectory.hpp index 64bd5015916a9d953d0a689f1ba22ceef667d603..984097ac3ab764d3398ae144a7cfd3fccddfe170 100644 --- a/corsika/framework/geometry/LeapFrogTrajectory.hpp +++ b/corsika/framework/geometry/LeapFrogTrajectory.hpp @@ -72,6 +72,10 @@ namespace corsika { ///! duration along potentially bend trajectory TimeType getDuration(double const u = 1) const; + ///! time at the start (u=0) or at the end (u=1) of the track of a particle + template <typename Particle> + TimeType getTime(Particle const& particle, double const u) const; + ///! total length along potentially bend trajectory LengthType getLength(double const u = 1) const; diff --git a/corsika/framework/geometry/StraightTrajectory.hpp b/corsika/framework/geometry/StraightTrajectory.hpp index 4bed6b7f6ead06789a5712b11148f509b7eb5e9c..a47da0db5b7adc82281dbd473678c27f38be31b4 100644 --- a/corsika/framework/geometry/StraightTrajectory.hpp +++ b/corsika/framework/geometry/StraightTrajectory.hpp @@ -84,6 +84,10 @@ namespace corsika { ///! duration along potentially bend trajectory TimeType getDuration(double const u = 1) const; + ///! time at the start (u=0) or at the end (u=1) of the track of a particle + template <typename Particle> + TimeType getTime(Particle const& particle, double const u) const; + ///! total length along potentially bend trajectory LengthType getLength(double const u = 1) const; diff --git a/tests/framework/testGeometry.cpp b/tests/framework/testGeometry.cpp index d9f2e7c29340bed32dfbc8f4a74158e8fea490ff..b820d946f6a3131d304fe4ba25076e1e263da8e0 100644 --- a/tests/framework/testGeometry.cpp +++ b/tests/framework/testGeometry.cpp @@ -21,6 +21,8 @@ #include <corsika/framework/geometry/Sphere.hpp> #include <corsika/framework/geometry/StraightTrajectory.hpp> +#include <corsika/setup/SetupStack.hpp> + #include <PhysicalUnitsCatch2.hpp> // namespace corsike::testing using namespace corsika; @@ -331,6 +333,23 @@ TEST_CASE("Geometry Box") { TEST_CASE("Geometry Trajectories") { CoordinateSystemPtr rootCS = get_root_CoordinateSystem(); Point r0(rootCS, {0_m, 0_m, 0_m}); + // Create a particle and and a stack so we can test .getTime() method + const Code particle{Code::Electron}; + setup::Stack stack; + // the mass of the particle + const auto pmass{get_mass(particle)}; + // set an arbitrary energy value + const HEPEnergyType E0{1_TeV}; + // compute the corresponding momentum + const HEPMomentumType P0{sqrt(E0 * E0 - pmass * pmass)}; + // create the momentum vector + const auto plab{MomentumVector(rootCS, {0_GeV, 0_GeV, P0})}; + // create an arbitrary location of the particle + const Point pos(rootCS, 50_m, 10_m, 80_m); + // add it finally to the stack + auto const particle1{stack.addParticle( + std::make_tuple(particle, calculate_kinetic_energy(plab.getNorm(), pmass), + plab.normalized(), pos, 0_ns))}; SECTION("Line") { SpeedType const V0 = 3_m / second; @@ -357,6 +376,8 @@ TEST_CASE("Geometry Trajectories") { auto const t = 1_s; StraightTrajectory base(line, t); CHECK(line.getPosition(t).getCoordinates() == base.getPosition(1.).getCoordinates()); + // test the getTime() method for straight trajectory + CHECK(base.getTime(particle1, 1) / 1_s == Approx(1)); CHECK((base.getDirection(0).getComponents(rootCS) - QuantityVector<dimensionless_d>{1, 0, 0}) @@ -375,6 +396,8 @@ TEST_CASE("Geometry Trajectories") { std::numeric_limits<TimeType::value_type>::infinity() * 1_s); base2.setDuration(10_s); CHECK(base2.getDuration() / 1_s == Approx(10)); + // test the getTime() method for straight trajectory + CHECK(base2.getTime(particle1, 0) / 1_s == Approx(0)); base2.setLength(1.3_m); CHECK(base2.getDuration() * V0 / meter == Approx(1.3)); @@ -406,6 +429,25 @@ TEST_CASE("Geometry Trajectories") { .getNorm() .magnitude() == Approx(0).margin(absMargin)); } + + SECTION("LeapFrog Trajectory") { + + // Create a velocity Vector + VelocityVector v0(rootCS, {5e+2_m / second, 5e+2_m / second, 5e+2_m / second}); + // Create a magnetic filed Vector + Vector B0(rootCS, 5_T, 5_T, 5_T); + // provide a k constant and a random time for the LeapFrog Trajectory + auto const k{1_m * ((1_m) / ((1_s * 1_s) * 1_V))}; + auto const t = 1e-12_s; + + // construct a LeapFrog Trajectory + LeapFrogTrajectory base(pos, v0, B0, k, t); + + // test the getTime() method for trajectories + CHECK((base.getTime(particle1, 1) - t) / 1_s == Approx(0)); + CHECK(base.getTime(particle1, 0) / 1_s == Approx(0)); + CHECK((base.getTime(particle1, 0) + t) / 1_s == Approx(1e-12)); + } } TEST_CASE("Distance between points") {