diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h index 485e4b734d960e77f332a3b09146215af0e7ac6e..e23ffc828165e5c7b041a38ea91f7157c051f9d6 100644 --- a/Framework/Geometry/BaseTrajectory.h +++ b/Framework/Geometry/BaseTrajectory.h @@ -9,13 +9,20 @@ namespace corsika::geometry { /*! - * Base class for trajectories. + * Interface / base class for trajectories. */ class BaseTrajectory { public: - //!< t for \f$ t = 0 \f$, the starting Point shall be returned. - virtual Point GetPosition(corsika::units::si::TimeType t) const = 0; + //!< for \f$ t = 0 \f$, the starting Point shall be returned. + virtual Point GetPosition(corsika::units::si::TimeType) const = 0; + + /*! + * returns the arc length between two points of the trajectory + * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1. + + virtual LengthType DistanceBetween(corsika::units::si::TimeType t1, + corsika::units::si::TimeType t2) const = 0; }; } // namespace corsika::geometry diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h index af92a68a5bb1e0bfae39390e46e14eafb7b68bf7..bd62d1919b5500ab22f2bad1eb2968a93dc90c0f 100644 --- a/Framework/Geometry/Helix.h +++ b/Framework/Geometry/Helix.h @@ -40,12 +40,17 @@ namespace corsika::geometry { , uPerp(vPerp.cross(vPar.normalized())) , radius(pvPar.norm() / abs(pOmegaC)) {} - Point GetPosition(corsika::units::si::TimeType t) const { + Point GetPosition(corsika::units::si::TimeType t) const override { return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC; } auto GetRadius() const { return radius; } + + LengthType DistanceBetween(corsika::units::si::TimeType t1, + corsika::units::si::TimeType t2) const override { + return (vPar + vPerp).norm() * (t2 - t1); + } }; } // namespace corsika::geometry diff --git a/Framework/Geometry/LineTrajectory.h b/Framework/Geometry/LineTrajectory.h index fcdc8c16ccdc24921d245ed0c6b529cd8d800876..ee65b5e54e9ec3743b12b6aa8a14a6965e6d4eeb 100644 --- a/Framework/Geometry/LineTrajectory.h +++ b/Framework/Geometry/LineTrajectory.h @@ -22,6 +22,12 @@ namespace corsika::geometry { Point GetPosition(corsika::units::si::TimeType t) const override { return r0 + v0 * t; } + + LengthType DistanceBetween(corsika::units::si::TimeType t1, + corsika::units::si::TimeType t2) const override { + assert(t2 >= t1); + return v0.norm() * (t2 - t1); + } }; } // namespace corsika::geometry diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h index 773fc2425d44e2fa9ee74f6862ca594367f78505..0ce2b536659d1662e08aff8a4ea233ade3f430f1 100644 --- a/Framework/Geometry/Trajectory.h +++ b/Framework/Geometry/Trajectory.h @@ -24,6 +24,10 @@ namespace corsika::geometry { Point GetPosition(double u) const { return GetPosition(fTEnd * u + fTStart * (1 - u)); } + + auto GetEndpoint() const { return GetPosition(fTEnd); } + + auto GetStartpoint() const { return GetPosition(fTStart); } }; } // namespace corsika::geometry diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc index 5419ad9a4c9092f78b034dbf9e13c2424912371d..3c7f63af445c847541fd2c3fc665aa8ae4d3198d 100644 --- a/Framework/Geometry/testGeometry.cc +++ b/Framework/Geometry/testGeometry.cc @@ -138,12 +138,17 @@ TEST_CASE("Trajectories") { BaseTrajectory const* base = &lineTrajectory; CHECK(lineTrajectory.GetPosition(2_s).GetCoordinates() == base->GetPosition(2_s).GetCoordinates()); + + CHECK(base->DistanceBetween(1_s, 2_s) / 1_m == Approx(1)); } SECTION("Helix") { Vector<SpeedType::dimension_type> const vPar( - rootCS, {0_m / second, 0_m / second, 4_m / second}), - vPerp(rootCS, {1_m / second, 0_m / second, 0_m / second}); + rootCS, {0_m / second, 0_m / second, 4_m / second}); + + Vector<SpeedType::dimension_type> const vPerp( + rootCS, {3_m / second, 0_m / second, 0_m / second}); + auto const omegaC = 2 * M_PI / 1_s; Helix const helix(r0, omegaC, vPar, vPerp); @@ -154,12 +159,14 @@ TEST_CASE("Trajectories") { .magnitude() == Approx(0).margin(absMargin)); CHECK((helix.GetPosition(0.25_s).GetCoordinates() - - QuantityVector<length_d>(-1_m / (2 * M_PI), -1_m / (2 * M_PI), 1_m)) + QuantityVector<length_d>(-3_m / (2 * M_PI), -3_m / (2 * M_PI), 1_m)) .norm() .magnitude() == Approx(0).margin(absMargin)); BaseTrajectory const* base = &helix; CHECK(helix.GetPosition(1234_s).GetCoordinates() == base->GetPosition(1234_s).GetCoordinates()); + + CHECK(base->DistanceBetween(1_s, 2_s) / 1_m == Approx(5)); } }