diff --git a/Framework/Geometry/BaseTrajectory.h b/Framework/Geometry/BaseTrajectory.h
new file mode 100644
index 0000000000000000000000000000000000000000..ec5e4062cfe2be63d1f7be83bdb4a3e02cc1d596
--- /dev/null
+++ b/Framework/Geometry/BaseTrajectory.h
@@ -0,0 +1,23 @@
+#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
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index c55295c897dcf0253293b3342d5042aa77cc8e26..a4fc43cb3364881d7f3e42939b0985b69df4afab 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -13,6 +13,9 @@ set (
Helix.h
BaseVector.h
QuantityVector.h
+ BaseTrajectory.h
+ LineTrajectory.h
+ Trajectory.h
)
set (
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
index c392356d46e64d63d156ca76f848845e916cde36..29bb49438f23387413094c8319486db8eea67e70 100644
--- a/Framework/Geometry/Helix.h
+++ b/Framework/Geometry/Helix.h
@@ -1,33 +1,38 @@
#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;
}
diff --git a/Framework/Geometry/LineTrajectory.h b/Framework/Geometry/LineTrajectory.h
index 4d56a598f58b8a5073b583654b1377c09d4d92f6..556da4642c52689c492fbd619863a3bfc0729d99 100644
--- a/Framework/Geometry/LineTrajectory.h
+++ b/Framework/Geometry/LineTrajectory.h
@@ -1,26 +1,27 @@
#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
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index c5af419ee4df97ec731a7fb0c95123fc2c843de5..671a6856092a0811aaf5d407ef68b0f7fbd37e4f 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/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
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 552dfde6f1341a7c0bbfb3903e27c8fbab5e1e05..28b4e7dc63d5ded1ecd6d047df442af4557e4f20 100644
--- a/Framework/Geometry/Sphere.h
+++ b/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();
}
diff --git a/Framework/Geometry/Trajectory.h b/Framework/Geometry/Trajectory.h
new file mode 100644
index 0000000000000000000000000000000000000000..8a640cacd027510a64d0c574c48e66dbe11fd1d9
--- /dev/null
+++ b/Framework/Geometry/Trajectory.h
@@ -0,0 +1,31 @@
+#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
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index af824df1c37a31182c8a31dd2d2fbe7fb64afd6f..d453b8f7b9caef7703a09971e6974ab36d9b80bc 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/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());
+ }
+}