Trajectory interface

Framework/Geometry/BaseTrajectory.h

+#ifndef _include_BASETRAJECTORY_H
+#define _include_BASETRAJECTORY_H
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <string>
+
+namespace corsika::geometry {
+
+  /*!
+   * Base class for trajectories.
+   */
+  class BaseTrajectory {
+
+  public:
+    //!< t for \f$ t = 0 \f$, the starting Point shall be returned.
+    virtual Point GetPosition(corsika::units::TimeType t) const = 0;
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/Geometry/CMakeLists.txt

@@ -13,6 +13,9 @@ set (
   Helix.h
   BaseVector.h
   QuantityVector.h
+  BaseTrajectory.h
+  LineTrajectory.h
+  Trajectory.h
   )
 
 set (

Framework/Geometry/Helix.h

 #ifndef _include_HELIX_H_
 #define _include_HELIX_H_
 
+#include <corsika/geometry/BaseTrajectory.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <cmath>
 
 namespace corsika::geometry {
-
-  using corsika::units::frequency_d;
-  using corsika::units::FrequencyType;
-  using corsika::units::quantity;
-  using corsika::units::SpeedType;
-  using corsika::units::TimeType;
-
-  class Helix // TODO: inherit from to-be-implemented "Trajectory"
-  {
-    using SpeedVec = Vector<SpeedType::dimension_type>;
+  /*!
+   * A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
+   * Point r0 and
+   * the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
+   * denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
+   * and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
+   * \f{align*}{
+        \vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
+        \vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
+     \f}
+   */
+
+  class Helix : public BaseTrajectory {
+    using VelocityVec = Vector<corsika::units::SpeedType::dimension_type>;
 
     Point const r0;
-    FrequencyType const omegaC;
-    SpeedVec const vPar;
-    SpeedVec vPerp, uPerp;
+    corsika::units::FrequencyType const omegaC;
+    VelocityVec const vPar;
+    VelocityVec const vPerp, uPerp;
 
-    LengthType const radius;
+    corsika::units::LengthType const radius;
 
   public:
-    Helix(Point const& pR0, quantity<frequency_d> pOmegaC, SpeedVec const& pvPar,
-          SpeedVec const& pvPerp)
+    Helix(Point const& pR0, corsika::units::FrequencyType pOmegaC,
+          VelocityVec const& pvPar, VelocityVec const& pvPerp)
         : r0(pR0)
         , omegaC(pOmegaC)
         , vPar(pvPar)
@@ -35,7 +40,7 @@ namespace corsika::geometry {
         , uPerp(vPerp.cross(vPar.normalized()))
         , radius(pvPar.norm() / abs(pOmegaC)) {}
 
-    auto GetPosition(TimeType t) const {
+    Point GetPosition(corsika::units::TimeType t) const {
       return r0 + vPar * t +
              (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
     }

Framework/Geometry/LineTrajectory.h

 #ifndef _include_LINETRAJECTORY_H
 #define _include_LINETRAJECTORY_H
 
-#include <Units/PhysicalUnits.h>
-#include <corsika/Point.h>
-#include <corsika/Vector.h>
+#include <corsika/geometry/BaseTrajectory.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
 
-namesapce corsika {
+namespace corsika::geometry {
 
-  class LineTrajectory // TODO: inherit from Trajectory
-  {
-    using SpeedVec = Vector<Speed::dimension_type>;
+  class LineTrajectory : public BaseTrajectory {
+    using VelocityVec = Vector<corsika::units::SpeedType::dimension_type>;
 
     Point const r0;
-    SpeedVec const v0;
+    VelocityVec const v0;
 
-    LineTrajectory(Point const& pR0, SpeedVec const& pV0)
-        : r0(r0)
+  public:
+    LineTrajectory(Point const& pR0, VelocityVec const& pV0)
+        : r0(pR0)
         , v0(pV0) {}
 
-    auto GetPosition(Time t) const { return r0 + v0 * t; }
+    Point GetPosition(corsika::units::TimeType t) const override { return r0 + v0 * t; }
   };
 
-} // end namesapce
+} // namespace corsika::geometry
 
 #endif

Framework/Geometry/QuantityVector.h

@@ -101,6 +101,8 @@ namespace corsika {
     auto& operator-() const { return QuantityVector<dim>(-eVector); }
 
     auto normalized() const { return (*this) * (1 / norm()); }
+
+    auto operator==(QuantityVector<dim> const& p) const { return eVector == p.eVector; }
   };
 
 } // end namespace corsika

Framework/Geometry/Sphere.h

@@ -15,6 +15,7 @@ namespace corsika::geometry {
         : center(pCenter)
         , radius(pRadius) {}
 
+    //! returns true if the Point p is within the sphere
     auto isInside(Point const& p) const {
       return radius * radius > (center - p).squaredNorm();
     }

Framework/Geometry/Trajectory.h

+#ifndef _include_TRAJECTORY_H
+#define _include_TRAJECTORY_H
+
+#include <corsika/geometry/BaseTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::geometry {
+
+  class Trajectory {
+    corsika::units::TimeType const fTStart, fTEnd;
+    BaseTrajectory const& fTrajectory;
+
+  public:
+    Trajectory(corsika::units::TimeType pTStart, corsika::units::TimeType pTEnd,
+               BaseTrajectory const& pTrajectory)
+        : fTStart(pTStart)
+        , fTEnd(pTEnd)
+        , fTrajectory(pTrajectory) {}
+
+    Point GetPosition(corsika::units::TimeType t) const {
+      return fTrajectory.GetPosition(t + fTStart);
+    }
+
+    Point GetPosition(double u) const {
+      return GetPosition(fTEnd * u + fTStart * (1 - u));
+    }
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/Geometry/testGeometry.cc

@@ -2,14 +2,17 @@
                           // cpp file
 #include <catch2/catch.hpp>
 
-#include <Geometry/CoordinateSystem.h>
-#include <Geometry/Point.h>
-#include <Geometry/Sphere.h>
+#include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/Helix.h>
+#include <corsika/geometry/LineTrajectory.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
+#include <corsika/geometry/Trajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <cmath>
 
-using namespace phys::units;
-using namespace phys::units::literals;
+using namespace corsika::geometry;
+using namespace corsika::units;
 
 TEST_CASE("transformations between CoordinateSystems") {
   CoordinateSystem rootCS;
@@ -17,10 +20,11 @@ TEST_CASE("transformations between CoordinateSystems") {
   REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
               .isApprox(EigenTransform::Identity()));
 
-  QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+  corsika::QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
   Point p1(rootCS, coordinates);
 
-  QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla, 0. * tesla};
+  corsika::QuantityVector<magnetic_flux_density_d> components{1. * tesla, 0. * tesla,
+                                                              0. * tesla};
   Vector<magnetic_flux_density_d> v1(rootCS, components);
 
   REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() == Approx(0));
@@ -32,7 +36,7 @@ TEST_CASE("transformations between CoordinateSystems") {
   }
 
   SECTION("translations") {
-    QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
+    corsika::QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
 
     CoordinateSystem translatedCS = rootCS.translate(translationVector);
 
@@ -52,13 +56,13 @@ TEST_CASE("transformations between CoordinateSystems") {
   }
 
   SECTION("multiple translations") {
-    QuantityVector<length_d> const tv1{0_m, 5_m, 0_m};
+    corsika::QuantityVector<length_d> const tv1{0_m, 5_m, 0_m};
     CoordinateSystem cs2 = rootCS.translate(tv1);
 
-    QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
+    corsika::QuantityVector<length_d> const tv2{3_m, 0_m, 0_m};
     CoordinateSystem cs3 = rootCS.translate(tv2);
 
-    QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
+    corsika::QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
     CoordinateSystem cs4 = cs3.translate(tv3);
 
     REQUIRE(cs4.GetReference()->GetReference() == &rootCS);
@@ -70,7 +74,7 @@ TEST_CASE("transformations between CoordinateSystems") {
   }
 
   SECTION("rotations") {
-    QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
+    corsika::QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
     double const angle = 90. / 180. * M_PI;
 
     CoordinateSystem rotatedCS = rootCS.rotate(axis, angle);
@@ -85,9 +89,9 @@ TEST_CASE("transformations between CoordinateSystems") {
   }
 
   SECTION("multiple rotations") {
-    QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km};
-    QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m};
-    QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m};
+    corsika::QuantityVector<length_d> const zAxis{0_m, 0_m, 1_km};
+    corsika::QuantityVector<length_d> const yAxis{0_m, 7_nm, 0_m};
+    corsika::QuantityVector<length_d> const xAxis{2_m, 0_nm, 0_m};
 
     double const angle = 90. / 180. * M_PI;
 
@@ -110,7 +114,26 @@ TEST_CASE("Sphere") {
 
   SECTION("isInside") {
     REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
-
     REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
   }
 }
+
+TEST_CASE("Trajectories") {
+  CoordinateSystem rootCS;
+  Point r0(rootCS, {0_m, 0_m, 0_m});
+
+  SECTION("Line") {
+    Vector<SpeedType::dimension_type> v0(rootCS,
+                                         {1_m / second, 0_m / second, 0_m / second});
+
+    LineTrajectory lineTrajectory(r0, v0);
+    REQUIRE((lineTrajectory.GetPosition(2_s).GetCoordinates() -
+             corsika::QuantityVector<length_d>(2_m, 0_m, 0_m))
+                .norm()
+                .magnitude() == Approx(0));
+
+    BaseTrajectory* base = &lineTrajectory;
+    REQUIRE(lineTrajectory.GetPosition(2_s).GetCoordinates() ==
+            base->GetPosition(2_s).GetCoordinates());
+  }
+}