From b5796d7171198b38701e4671c1fc32631f9936a1 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 13 Aug 2018 18:36:06 +0200
Subject: [PATCH] added Helix
---
Framework/Geometry/CMakeLists.txt | 2 +-
Framework/Geometry/Helix.h | 31 ++++++++++++++++++++
Framework/Geometry/QuantityVector.h | 17 +++++++++++
Framework/Geometry/Vector.h | 44 ++++++++++++++++++++++-------
4 files changed, 83 insertions(+), 11 deletions(-)
create mode 100644 Framework/Geometry/Helix.h
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index eb7fb3650..16dce7305 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -1,6 +1,6 @@
set (GEOMETRY_SOURCES CoordinateSystem.cc)
-set (GEOMETRY_HEADERS Vector.h Point.h Sphere.h CoordinateSystem.h)
+set (GEOMETRY_HEADERS Vector.h Point.h Sphere.h CoordinateSystem.h Helix.h)
add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
diff --git a/Framework/Geometry/Helix.h b/Framework/Geometry/Helix.h
new file mode 100644
index 000000000..9dacfb9dd
--- /dev/null
+++ b/Framework/Geometry/Helix.h
@@ -0,0 +1,31 @@
+#ifndef _include_HELIX_H_
+#define _include_HELIX_H_
+
+#include <Geometry/Vector.h>
+#include <Geometry/Point.h>
+#include <Units/PhysicalUnits.h>
+
+class Helix // TODO: inherit from to-be-implemented "Trajectory"
+{
+ using SpeedVec = Vector<phys::units::speed_d>;
+
+ Point const r0;
+ phys::units::quantity<phys::units::frequency_d> const omegaC;
+ SpeedVec const vPar;
+ SpeedVec vPerp, uPerp;
+
+public:
+ Helix(Point const pR0, phys::units::quantity<phys::units::frequency_d> pOmegaC,
+ SpeedVec const pvPar, SpeedVec const pvPerp) :
+ r0(pR0), omegaC(pOmegaC), vPar(pvPar), vPerp(pvPerp), uPerp(vPar.normalized().cross(vPerp))
+ {
+
+ }
+
+ Point getPosition(phys::units::quantity<phys::units::time_interval_d> t) const
+ {
+ return r0 + vPar * t + (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
+ }
+};
+
+#endif
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 025469620..817671512 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -70,11 +70,28 @@ public:
}
}
+ template <typename ScalarDim>
+ auto operator/(phys::units::quantity<ScalarDim, double> const p) const
+ {
+ return (*this) * (1 / p);
+ }
+
auto operator*(double const p) const
{
return QuantityVector<dim>(eVector * p);
}
+ auto operator/(double const p) const
+ {
+ return QuantityVector<dim>(eVector / p);
+ }
+
+ auto& operator/=(double const p) const
+ {
+ eVector /= p;
+ return *this;
+ }
+
auto& operator*=(double const p)
{
eVector *= p;
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 5670aaaf1..2fe62e04b 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -92,23 +92,34 @@ public:
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const
{
- using ResQuantity = decltype(BaseVector<dim>::qVector * p);
- if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
+ using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
+
+ if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
{
return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
else
{
- using res_dim = typename decltype(BaseVector<dim>::qVector * p)::dimension;
- return Vector<res_dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
+ return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
}
+ template <typename ScalarDim>
+ auto operator/(phys::units::quantity<ScalarDim, double> const p) const
+ {
+ return (*this) * (1 / p);
+ }
+
auto operator*(double const p) const
{
return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
+ auto operator/(double const p) const
+ {
+ return Vector<dim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector / p);
+ }
+
auto& operator+=(Vector<dim> const& pVec)
{
BaseVector<dim>::qVector += pVec.getComponents(*BaseVector<dim>::cs);
@@ -131,12 +142,25 @@ public:
return (*this) * (1 / norm());
}
- //~ template <typename dim2>
- //~ auto operator*(Vector<dim2> const& pVec)
- //~ {
- //~ auto constexpr resulting_dim = dimension
- //~ return Vector<
- //~ }
+ template <typename dim2>
+ auto cross(Vector<dim2> pV) const
+ {
+ auto const c1 = getComponents().eVector;
+ auto const c2 = pV.getComponents(*BaseVector<dim>::cs).eVector;
+ auto const bareResult = c1.cross(c2);
+
+ using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+
+ if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless, not a "Quantity" anymore
+ {
+ return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
+ }
+ else
+ {
+ return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs, bareResult);
+ }
+ }
+
};
#endif
--
GitLab