Framework/Geometry/BaseVector.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_BASEVECTOR_H_
-#define _include_BASEVECTOR_H_
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/QuantityVector.h>
-
-namespace corsika::geometry {
-
-  /*!
-   * Common base class for Vector and Point. Currently it does basically nothing.
-   */
-
-  template <typename dim>
-  class BaseVector {
-  protected:
-    QuantityVector<dim> qVector;
-    CoordinateSystem const* cs;
-
-  public:
-    BaseVector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
-        : qVector(pQVector)
-        , cs(&pCS) {}
-
-    auto const& GetCoordinateSystem() const { return *cs; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/CMakeLists.txt

-
-set (
-  GEOMETRY_SOURCES
-  CoordinateSystem.cc
-  )
-
-set (
-  GEOMETRY_HEADERS
-  Vector.h
-  Point.h
-  Line.h
-  Sphere.h
-  Volume.h
-  CoordinateSystem.h
-  RootCoordinateSystem.h
-  Helix.h
-  BaseVector.h
-  QuantityVector.h
-  Trajectory.h
-  FourVector.h
-  )
-
-set (
-  GEOMETRY_NAMESPACE
-  corsika/geometry
-  )
-
-add_library (CORSIKAgeometry STATIC ${GEOMETRY_SOURCES})
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAgeometry ${GEOMETRY_NAMESPACE} ${GEOMETRY_HEADERS})
-
-set_target_properties (
-  CORSIKAgeometry
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${GEOMETRY_HEADERS}"
-  )
-
-# target dependencies on other libraries (also the header onlys)
-target_link_libraries (
-  CORSIKAgeometry
-  CORSIKAunits
-  )
-
-target_include_directories (
-  CORSIKAgeometry
-  PUBLIC    ${EIGEN3_INCLUDE_DIR}
-  INTERFACE ${EIGEN3_INCLUDE_DIR}
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/include>
-  )
-
-install (
-  TARGETS CORSIKAgeometry
-  LIBRARY DESTINATION lib
-  ARCHIVE DESTINATION lib
-  PUBLIC_HEADER DESTINATION include/${GEOMETRY_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-add_executable (testGeometry testGeometry.cc)
-target_link_libraries (
-  testGeometry
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-CORSIKA_ADD_TEST(testGeometry)
-
-add_executable (testFourVector testFourVector.cc) 
-target_link_libraries (
-  testFourVector
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-CORSIKA_ADD_TEST(testFourVector)

Framework/Geometry/CoordinateSystem.cc

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <stdexcept>
-
-using namespace corsika::geometry;
-
-/**
- * returns the transformation matrix necessary to transform primitives with coordinates
- * in \a pFrom to \a pTo, e.g.
- * \f$ \vec{v}^{\text{(to)}} = \mathcal{M} \vec{v}^{\text{(from)}} \f$
- * (\f$ \vec{v}^{(.)} \f$ denotes the coordinates/components of the component in
- * the indicated CoordinateSystem).
- */
-EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom,
-                                                   CoordinateSystem const& pTo) {
-  CoordinateSystem const* a{&pFrom};
-  CoordinateSystem const* b{&pTo};
-  CoordinateSystem const* commonBase{nullptr};
-
-  while (a != b && b != nullptr) {
-    a = &pFrom;
-
-    while (a != b && a != nullptr) { a = a->GetReference(); }
-
-    if (a == b) break;
-
-    b = b->GetReference();
-  }
-
-  if (a == b && a != nullptr) {
-    commonBase = a;
-
-  } else {
-    throw std::runtime_error("no connection between coordinate systems found!");
-  }
-
-  EigenTransform t = EigenTransform::Identity();
-  auto* p = &pFrom;
-
-  while (p != commonBase) {
-    t = p->GetTransform() * t;
-    p = p->GetReference();
-  }
-
-  p = &pTo;
-
-  while (p != commonBase) {
-    t = t * p->GetTransform().inverse(Eigen::TransformTraits::Isometry);
-    p = p->GetReference();
-  }
-
-  return t;
-}

Framework/Geometry/CoordinateSystem.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_COORDINATESYSTEM_H_
-#define _include_COORDINATESYSTEM_H_
-
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <Eigen/Dense>
-#include <stdexcept>
-
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
-
-namespace corsika::geometry {
-
-  class RootCoordinateSystem;
-
-  using corsika::units::si::length_d;
-
-  class CoordinateSystem {
-    CoordinateSystem const* reference = nullptr;
-    EigenTransform transf;
-
-    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
-        : reference(&reference)
-        , transf(transf) {}
-
-    CoordinateSystem()
-        : // for creating the root CS
-        transf(EigenTransform::Identity()) {}
-
-  protected:
-    static auto CreateCS() { return CoordinateSystem(); }
-    friend corsika::geometry::RootCoordinateSystem; /// this is the only class that can
-                                                    /// create ONE unique root CS
-
-  public:
-    static EigenTransform GetTransformation(CoordinateSystem const& c1,
-                                            CoordinateSystem const& c2);
-
-    auto& operator=(const CoordinateSystem& pCS) {
-      reference = pCS.reference;
-      transf = pCS.transf;
-      return *this;
-    }
-
-    auto translate(QuantityVector<length_d> vector) const {
-      EigenTransform const translation{EigenTranslation(vector.eVector)};
-
-      return CoordinateSystem(*this, translation);
-    }
-
-    auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
-      if (axis.eVector.isZero()) {
-        throw std::runtime_error("null-vector given as axis parameter");
-      }
-
-      EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
-
-      return CoordinateSystem(*this, rotation);
-    }
-
-    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
-                            QuantityVector<phys::units::length_d> axis, double angle) {
-      if (axis.eVector.isZero()) {
-        throw std::runtime_error("null-vector given as axis parameter");
-      }
-
-      EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *
-                                  EigenTranslation(translation.eVector)};
-
-      return CoordinateSystem(*this, transf);
-    }
-
-    auto const* GetReference() const { return reference; }
-
-    auto const& GetTransform() const { return transf; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/FourVector.h

-#ifndef _include_corsika_framework_geometry_fourvector_h_
-#define _include_corsika_framework_geometry_fourvector_h_
-
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <iostream>
-#include <type_traits>
-
-namespace corsika::geometry {
-
-  /**
-     FourVector supports "full" units, e.g. E in [GeV/c] and p in [GeV],
-     or also t in [s] and r in [m], etc.
-
-     However, for HEP applications it is also possible to use E and p
-     both in [GeV].
-
-     The FourVector can return NormSqr and Norm, whereas Norm is
-     sqrt(abs(NormSqr)). The physical units are always calculated and
-     returned properly.
-
-     FourVector can also return if it is TimeLike, SpaceLike or PhotonLike.
-
-     When a FourVector is initialized with a lvalue reference, this is
-     also used for the internal storage, which should lead to complete
-     disappearance of the FourVector class during optimization.
-   */
-
-  template <typename TimeType, typename SpaceVecType>
-  class FourVector {
-
-  public:
-    using SpaceType = typename std::decay<SpaceVecType>::type::Quantity;
-
-    //! check the types and the physical units here:
-    static_assert(
-        std::is_same<typename std::decay<TimeType>::type, SpaceType>::value ||
-            std::is_same<typename std::decay<TimeType>::type,
-                         decltype(std::declval<SpaceType>() / corsika::units::si::meter *
-                                  corsika::units::si::second)>::value,
-        "Units of time-like and space-like coordinates must either be idential "
-        "(e.g. GeV) or [E/c]=[p]");
-
-  public:
-    FourVector(const TimeType& eT, const SpaceVecType& eS)
-        : fTimeLike(eT)
-        , fSpaceLike(eS) {}
-
-    TimeType GetTimeLikeComponent() const { return fTimeLike; }
-    SpaceVecType& GetSpaceLikeComponents() { return fSpaceLike; }
-    const SpaceVecType& GetSpaceLikeComponents() const { return fSpaceLike; }
-
-    auto GetNormSqr() const { return GetTimeSquared() - fSpaceLike.squaredNorm(); }
-
-    SpaceType GetNorm() const { return sqrt(abs(GetNormSqr())); }
-
-    bool IsTimelike() const {
-      return GetTimeSquared() < fSpaceLike.squaredNorm();
-    } //! Norm2 < 0
-
-    bool IsSpacelike() const {
-      return GetTimeSquared() > fSpaceLike.squaredNorm();
-    } //! Norm2 > 0
-
-    /* this is not numerically stable
-    bool IsPhotonlike() const {
-      return GetTimeSquared() == fSpaceLike.squaredNorm();
-    } //! Norm2 == 0
-    */
-
-    FourVector& operator+=(const FourVector& b) {
-      fTimeLike += b.fTimeLike;
-      fSpaceLike += b.fSpaceLike;
-      return *this;
-    }
-
-    FourVector& operator-=(const FourVector& b) {
-      fTimeLike -= b.fTimeLike;
-      fSpaceLike -= b.fSpaceLike;
-      return *this;
-    }
-
-    FourVector& operator*=(const double b) {
-      fTimeLike *= b;
-      fSpaceLike *= b;
-      return *this;
-    }
-
-    FourVector& operator/=(const double b) {
-      fTimeLike /= b;
-      fSpaceLike.GetComponents() /= b; // TODO: WHY IS THIS??????
-      return *this;
-    }
-
-    FourVector& operator/(const double b) {
-      *this /= b;
-      return *this;
-    }
-
-    /**
-       Note that the product between two 4-vectors assumes that you use
-       the same "c" convention for both. Only the LHS vector is checked
-       for this. You cannot mix different conventions due to
-       unit-checking.
-     */
-    SpaceType operator*(const FourVector& b) {
-      if constexpr (std::is_same<typename std::decay<TimeType>::type,
-                                 decltype(std::declval<SpaceType>() /
-                                          corsika::units::si::meter *
-                                          corsika::units::si::second)>::value)
-        return fTimeLike * b.fTimeLike *
-                   (corsika::units::constants::c * corsika::units::constants::c) -
-               fSpaceLike.norm();
-      else
-        return fTimeLike * fTimeLike - fSpaceLike.norm();
-    }
-
-  private:
-    /**
-       This function is automatically compiled to use of ignore the
-       extra factor of "c" for the time-like quantity
-     */
-    auto GetTimeSquared() const {
-      if constexpr (std::is_same<typename std::decay<TimeType>::type,
-                                 decltype(std::declval<SpaceType>() /
-                                          corsika::units::si::meter *
-                                          corsika::units::si::second)>::value)
-        return fTimeLike * fTimeLike *
-               (corsika::units::constants::c * corsika::units::constants::c);
-      else
-        return fTimeLike * fTimeLike;
-    }
-
-  protected:
-    //! the data members
-    TimeType fTimeLike;
-    SpaceVecType fSpaceLike;
-
-    //! the friends: math operators
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator+(const FourVector<T, U>&, const FourVector<T, U>&);
-
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator-(const FourVector<T, U>&, const FourVector<T, U>&);
-
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator*(const FourVector<T, U>&, const double);
-
-    template <typename T, typename U>
-    friend FourVector<typename std::decay<T>::type, typename std::decay<U>::type>
-    operator/(const FourVector<T, U>&, const double);
-  };
-
-  /**
-      The math operator+
-   */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator+(const FourVector<TimeType, SpaceVecType>& a,
-            const FourVector<TimeType, SpaceVecType>& b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(
-        a.fTimeLike + b.fTimeLike, a.fSpaceLike + b.fSpaceLike);
-  }
-
-  /**
-     The math operator-
-  */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator-(const FourVector<TimeType, SpaceVecType>& a,
-            const FourVector<TimeType, SpaceVecType>& b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(
-        a.fTimeLike - b.fTimeLike, a.fSpaceLike - b.fSpaceLike);
-  }
-
-  /**
-     The math operator*
-  */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator*(const FourVector<TimeType, SpaceVecType>& a, const double b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(a.fTimeLike * b,
-                                                               a.fSpaceLike * b);
-  }
-
-  /**
-      The math operator/
-   */
-  template <typename TimeType, typename SpaceVecType>
-  inline FourVector<typename std::decay<TimeType>::type,
-                    typename std::decay<SpaceVecType>::type>
-  operator/(const FourVector<TimeType, SpaceVecType>& a, const double b) {
-    return FourVector<typename std::decay<TimeType>::type,
-                      typename std::decay<SpaceVecType>::type>(a.fTimeLike / b,
-                                                               a.fSpaceLike / b);
-  }
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Helix.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_HELIX_H_
-#define _include_HELIX_H_
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <cmath>
-
-namespace corsika::geometry {
-  /*!
-   * 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 {
-
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    corsika::units::si::FrequencyType const omegaC;
-    VelocityVec const vPar;
-    VelocityVec const vPerp, uPerp;
-
-    corsika::units::si::LengthType const radius;
-
-  public:
-    Helix(Point const& pR0, corsika::units::si::FrequencyType pOmegaC,
-          VelocityVec const& pvPar, VelocityVec const& pvPerp)
-        : r0(pR0)
-        , omegaC(pOmegaC)
-        , vPar(pvPar)
-        , vPerp(pvPerp)
-        , uPerp(vPerp.cross(vPar.normalized()))
-        , radius(pvPar.norm() / abs(pOmegaC)) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const {
-      return r0 + vPar * t +
-             (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
-    }
-
-    Point PositionFromArclength(corsika::units::si::LengthType l) const {
-      return GetPosition(TimeFromArclength(l));
-    }
-
-    auto GetRadius() const { return radius; }
-
-    corsika::units::si::LengthType ArcLength(corsika::units::si::TimeType t1,
-                                             corsika::units::si::TimeType t2) const {
-      return (vPar + vPerp).norm() * (t2 - t1);
-    }
-
-    corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType l) const {
-      return l / (vPar + vPerp).norm();
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Line.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_LINETRAJECTORY_H
-#define _include_LINETRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  class Line {
-
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    VelocityVec const v0;
-
-  public:
-    Line(Point const& pR0, VelocityVec const& pV0)
-        : r0(pR0)
-        , v0(pV0) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
-
-    Point PositionFromArclength(corsika::units::si::LengthType l) const {
-      return r0 + v0.normalized() * l;
-    }
-
-    LengthType ArcLength(corsika::units::si::TimeType t1,
-                         corsika::units::si::TimeType t2) const {
-      return v0.norm() * (t2 - t1);
-    }
-
-    corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType t) const {
-      return t / v0.norm();
-    }
-
-    auto GetR0() const { return r0; }
-    auto GetV0() const { return v0; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Point.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_POINT_H_
-#define _include_POINT_H_
-
-#include <corsika/geometry/BaseVector.h>
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  using corsika::units::si::length_d;
-  using corsika::units::si::LengthType;
-
-  /*!
-   * A Point represents a point in position space. It is defined by its
-   * coordinates with respect to some CoordinateSystem.
-   */
-  class Point : public BaseVector<length_d> {
-  public:
-    Point(CoordinateSystem const& pCS, QuantityVector<length_d> pQVector)
-        : BaseVector<length_d>(pCS, pQVector) {}
-
-    Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
-        : BaseVector<length_d>(cs, {x, y, z}) {}
-
-    // TODO: this should be private or protected, we don NOT want to expose numbers
-    // without reference to outside:
-    auto GetCoordinates() const { return BaseVector<length_d>::qVector; }
-
-    /// this always returns a QuantityVector as triple
-    auto GetCoordinates(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<length_d>::cs) {
-        return BaseVector<length_d>::qVector;
-      } else {
-        return QuantityVector<length_d>(
-            CoordinateSystem::GetTransformation(*BaseVector<length_d>::cs, pCS) *
-            BaseVector<length_d>::qVector.eVector);
-      }
-    }
-
-    /*!
-     * transforms the Point into another CoordinateSystem by changing its
-     * coordinates interally
-     */
-    void rebase(CoordinateSystem const& pCS) {
-      BaseVector<length_d>::qVector = GetCoordinates(pCS);
-      BaseVector<length_d>::cs = &pCS;
-    }
-
-    Point operator+(Vector<length_d> const& pVec) const {
-      return Point(*BaseVector<length_d>::cs,
-                   GetCoordinates() + pVec.GetComponents(*BaseVector<length_d>::cs));
-    }
-
-    /*!
-     * returns the distance Vector between two points
-     */
-    Vector<length_d> operator-(Point const& pB) const {
-      auto& cs = *BaseVector<length_d>::cs;
-      return Vector<length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/QuantityVector.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_QUANTITYVECTOR_H_
-#define _include_QUANTITYVECTOR_H_
-
-#include <corsika/units/PhysicalUnits.h>
-
-#include <Eigen/Dense>
-
-#include <iostream>
-#include <utility>
-
-namespace corsika::geometry {
-
-  /*!
-   * A QuantityVector is a three-component container based on Eigen::Vector3d
-   * with a phys::units::si::dimension. Arithmethic operators are defined that
-   * propagate the dimensions by dimensional analysis.
-   */
-
-  template <typename dim>
-  class QuantityVector {
-  public:
-    using Quantity = phys::units::quantity<dim, double>; //< the phys::units::quantity
-                                                         // corresponding to the dimension
-
-  public:
-    Eigen::Vector3d eVector; //!< the actual container where the raw numbers are stored
-
-    typedef dim dimension; //!< should be a phys::units::dimension
-
-    QuantityVector(Quantity a, Quantity b, Quantity c)
-        : eVector{a.magnitude(), b.magnitude(), c.magnitude()} {}
-
-    QuantityVector(Eigen::Vector3d pBareVector)
-        : eVector(pBareVector) {}
-
-    auto operator[](size_t index) const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
-    }
-
-    auto norm() const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
-    }
-
-    auto squaredNorm() const {
-      using QuantitySquared =
-          decltype(std::declval<Quantity>() * std::declval<Quantity>());
-      return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
-    }
-
-    auto operator+(QuantityVector<dim> const& pQVec) const {
-      return QuantityVector<dim>(eVector + pQVec.eVector);
-    }
-
-    auto operator-(QuantityVector<dim> const& pQVec) const {
-      return QuantityVector<dim>(eVector - pQVec.eVector);
-    }
-
-    template <typename ScalarDim>
-    auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
-      using ResQuantity = phys::units::detail::Product<ScalarDim, dim, double, double>;
-
-      if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
-                                                              // a "Quantity" anymore
-      {
-        return QuantityVector<phys::units::dimensionless_d>(eVector * p.magnitude());
-      } else {
-        return QuantityVector<typename ResQuantity::dimension_type>(eVector *
-                                                                    p.magnitude());
-      }
-    }
-
-    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) {
-      eVector /= p;
-      return *this;
-    }
-
-    auto& operator*=(double const p) {
-      eVector *= p;
-      return *this;
-    }
-
-    auto& operator+=(QuantityVector<dim> const& pQVec) {
-      eVector += pQVec.eVector;
-      return *this;
-    }
-
-    auto& operator-=(QuantityVector<dim> const& pQVec) {
-      eVector -= pQVec.eVector;
-      return *this;
-    }
-
-    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; }
-  };
-
-  template <typename dim>
-  auto& operator<<(std::ostream& os, corsika::geometry::QuantityVector<dim> qv) {
-    using Quantity = phys::units::quantity<dim, double>;
-
-    os << '(' << qv.eVector(0) << ' ' << qv.eVector(1) << ' ' << qv.eVector(2) << ") "
-       << phys::units::to_unit_symbol<dim, double>(
-              Quantity(phys::units::detail::magnitude_tag, 1));
-    return os;
-  }
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/RootCoordinateSystem.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_corsika_geometry_rootcoordinatesystem_h_
-#define _include_corsika_geometry_rootcoordinatesystem_h_
-
-#include <corsika/utl/Singleton.h>
-
-#include <corsika/geometry/CoordinateSystem.h>
-
-/*!
- * This is the only way to get a root-coordinate system, and it is a
- * singleton. All other CoordinateSystems must be relative to the
- * RootCoordinateSystem
- */
-
-namespace corsika::geometry {
-
-  class RootCoordinateSystem : public corsika::utl::Singleton<RootCoordinateSystem> {
-
-    friend class corsika::utl::Singleton<RootCoordinateSystem>;
-
-  protected:
-    RootCoordinateSystem() {}
-
-  public:
-    corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() { return fRootCS; }
-    const corsika::geometry::CoordinateSystem& GetRootCoordinateSystem() const {
-      return fRootCS;
-    }
-
-  private:
-    corsika::geometry::CoordinateSystem fRootCS; // THIS IS IT
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Sphere.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_SPHERE_H_
-#define _include_SPHERE_H_
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Volume.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  class Sphere : public Volume {
-    Point const fCenter;
-    LengthType const fRadius;
-
-  public:
-    Sphere(Point const& pCenter, LengthType const pRadius)
-        : fCenter(pCenter)
-        , fRadius(pRadius) {}
-
-    //! returns true if the Point p is within the sphere
-    bool Contains(Point const& p) const override {
-      return fRadius * fRadius > (fCenter - p).squaredNorm();
-    }
-
-    auto& GetCenter() const { return fCenter; }
-    auto GetRadius() const { return fRadius; }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Trajectory.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_TRAJECTORY_H
-#define _include_TRAJECTORY_H
-
-#include <corsika/geometry/Point.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  template <typename T>
-  class Trajectory : public T {
-
-    corsika::units::si::TimeType fTimeLength;
-
-  public:
-    using T::ArcLength;
-    using T::GetPosition;
-
-    Trajectory(T const& theT, corsika::units::si::TimeType timeLength)
-        : T(theT)
-        , fTimeLength(timeLength) {}
-
-    /*Point GetPosition(corsika::units::si::TimeType t) const {
-      return fTraj.GetPosition(t + fTStart);
-      }*/
-
-    Point GetPosition(double u) const { return T::GetPosition(fTimeLength * u); }
-
-    corsika::units::si::TimeType GetDuration() const { return fTimeLength; }
-
-    corsika::units::si::LengthType GetDistance(corsika::units::si::TimeType t) const {
-      assert(t > fTimeLength);
-      assert(t >= 0 * corsika::units::si::second);
-      return T::ArcLength(0, t);
-    }
-
-    void LimitEndTo(corsika::units::si::LengthType limit) {
-      fTimeLength = T::TimeFromArclength(limit);
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Vector.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_VECTOR_H_
-#define _include_VECTOR_H_
-
-#include <corsika/geometry/BaseVector.h>
-#include <corsika/geometry/QuantityVector.h>
-
-#include <corsika/units/PhysicalUnits.h>
-
-/*!
- * A Vector represents a 3-vector in Euclidean space. It is defined by components
- * given in a specific CoordinateSystem. It has a physical dimension ("unit")
- * as part of its type, so you cannot mix up e.g. electric with magnetic fields
- * (but you could calculate their cross-product to get an energy flux vector).
- *
- * When transforming coordinate systems, a Vector is subject to the rotational
- * part only and invariant under translations.
- */
-
-namespace corsika::geometry {
-
-  template <typename dim>
-  class Vector : public BaseVector<dim> {
-  public:
-    using Quantity = phys::units::quantity<dim, double>;
-
-  public:
-    Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector)
-        : BaseVector<dim>(pCS, pQVector) {}
-
-    Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z)
-        : BaseVector<dim>(cs, QuantityVector<dim>(x, y, z)) {}
-
-    /*!
-     * returns a QuantityVector with the components given in the "home"
-     * CoordinateSystem of the Vector
-     */
-    auto GetComponents() const { return BaseVector<dim>::qVector; }
-
-    /*!
-     * returns a QuantityVector with the components given in an arbitrary
-     * CoordinateSystem
-     */
-    auto GetComponents(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<dim>::cs) {
-        return BaseVector<dim>::qVector;
-      } else {
-        return QuantityVector<dim>(
-            CoordinateSystem::GetTransformation(*BaseVector<dim>::cs, pCS).linear() *
-            BaseVector<dim>::qVector.eVector);
-      }
-    }
-
-    /*!
-     * transforms the Vector into another CoordinateSystem by changing
-     * its components internally
-     */
-    void rebase(CoordinateSystem const& pCS) {
-      BaseVector<dim>::qVector = GetComponents(pCS);
-      BaseVector<dim>::cs = &pCS;
-    }
-
-    /*!
-     * returns the norm/length of the Vector. Before using this method,
-     * think about whether squaredNorm() might be cheaper for your computation.
-     */
-    auto norm() const { return BaseVector<dim>::qVector.norm(); }
-
-    /*!
-     * returns the squared norm of the Vector. Before using this method,
-     * think about whether norm() might be cheaper for your computation.
-     */
-    auto squaredNorm() const { return BaseVector<dim>::qVector.squaredNorm(); }
-
-    /*!
-     * returns a Vector \f$ \vec{v}_{\parallel} \f$ which is the parallel projection
-     * of this vector \f$ \vec{v}_1 \f$ along another Vector \f$ \vec{v}_2 \f$ given by
-     *   \f[
-     *     \vec{v}_{\parallel} = \frac{\vec{v}_1 \cdot \vec{v}_2}{\vec{v}_2^2} \vec{v}_2
-     *   \f]
-     */
-    template <typename dim2>
-    auto parallelProjectionOnto(Vector<dim2> const& pVec,
-                                CoordinateSystem const& pCS) const {
-      auto const ourCompVec = GetComponents(pCS);
-      auto const otherCompVec = pVec.GetComponents(pCS);
-      auto const& a = ourCompVec.eVector;
-      auto const& b = otherCompVec.eVector;
-
-      return Vector<dim>(pCS, QuantityVector<dim>(b * ((a.dot(b)) / b.squaredNorm())));
-    }
-
-    template <typename dim2>
-    auto parallelProjectionOnto(Vector<dim2> const& pVec) const {
-      return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
-    }
-
-    auto operator+(Vector<dim> const& pVec) const {
-      auto const components =
-          GetComponents(*BaseVector<dim>::cs) + pVec.GetComponents(*BaseVector<dim>::cs);
-      return Vector<dim>(*BaseVector<dim>::cs, components);
-    }
-
-    auto operator-(Vector<dim> const& pVec) const {
-      auto const components = GetComponents() - pVec.GetComponents(*BaseVector<dim>::cs);
-      return Vector<dim>(*BaseVector<dim>::cs, components);
-    }
-
-    auto& operator*=(double const p) {
-      BaseVector<dim>::qVector *= p;
-      return *this;
-    }
-
-    template <typename ScalarDim>
-    auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
-      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 {
-        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);
-      return *this;
-    }
-
-    auto& operator-=(Vector<dim> const& pVec) {
-      BaseVector<dim>::qVector -= pVec.GetComponents(*BaseVector<dim>::cs);
-      return *this;
-    }
-
-    auto& operator-() const {
-      return Vector<dim>(*BaseVector<dim>::cs, -BaseVector<dim>::qVector);
-    }
-
-    auto normalized() const { return (*this) * (1 / norm()); }
-
-    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);
-      }
-    }
-
-    template <typename dim2>
-    auto dot(Vector<dim2> pV) const {
-      auto const c1 = GetComponents().eVector;
-      auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
-      auto const bareResult = c1.dot(c2);
-
-      using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
-
-      return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Volume.h

-
-/**
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_VOLUME_H_
-#define _include_VOLUME_H_
-
-#include <corsika/geometry/Point.h>
-
-namespace corsika::geometry {
-
-  class Volume {
-
-  public:
-    //! returns true if the Point p is within the volume
-    virtual bool Contains(Point const& p) const = 0;
-
-    virtual ~Volume() = default;
-  };
-
-} // namespace corsika::geometry
-
-#endif