From eb1c270d55c7b80ad132e02b71619d5bd73091a1 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 20 Sep 2018 17:14:11 +0200
Subject: [PATCH] added DistanceBetween() for trajectories

---
 Framework/Geometry/BaseTrajectory.h | 13 ++++++++++---
 Framework/Geometry/Helix.h          |  7 ++++++-
 Framework/Geometry/LineTrajectory.h |  6 ++++++
 Framework/Geometry/Trajectory.h     |  4 ++++
 Framework/Geometry/testGeometry.cc  | 13 ++++++++++---
 5 files changed, 36 insertions(+), 7 deletions(-)

diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
index 485e4b73..e23ffc82 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 af92a68a..bd62d191 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 fcdc8c16..ee65b5e5 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 773fc242..0ce2b536 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 5419ad9a..3c7f63af 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));
   }
 }
-- 
GitLab