Framework/CMakeLists.txt

-
-add_subdirectory (Units)
-add_subdirectory (Geometry)
-add_subdirectory (Particles)
-add_subdirectory (Logging)
-add_subdirectory (StackInterface)
-add_subdirectory (ProcessSequence)
-add_subdirectory (Cascade)
-add_subdirectory (Random)

Framework/Cascade/CMakeLists.txt

-
-# namespace of library -> location of header files
-set (
-  CORSIKAcascade_NAMESPACE
-  corsika/cascade
-  )
-
-# header files of this library
-set (
-  CORSIKAcascade_HEADERS
-  Cascade.h
-  )
-
-#set (
-#  CORSIKAcascade_SOURCES
-#  Cascade.cc
-#  )
-
-#add_library (CORSIKAcascade STATIC ${CORSIKAcascade_SOURCES})
-add_library (CORSIKAcascade INTERFACE)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
-
-#target_link_libraries (
-#  CORSIKAcascade
-#  CORSIKAparticles
-#  CORSIKAunits
-#  CORSIKAthirdparty # for catch2
-#  )
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAcascade
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# install library
-install (
-  FILES ${CORSIKAcascade_HEADERS}
-  DESTINATION include/${CORSIKAcascade_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-add_executable (
-  testCascade
-  testCascade.cc
-  )
-
-target_link_libraries (
-  testCascade
-  CORSIKAcascade
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAprocesssequence
-  SuperStupidStack
-  CORSIKAunits
-  CORSIKAthirdparty # for catch2
-  )
-
-add_test (
-  NAME testCascade
-  COMMAND testCascade
-  )

Framework/Cascade/Cascade.cc

-#include <corsika/cascade/Cascade.h>
-
-using namespace corsika::cascade;
-
-template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
-Cascade<ProcessList, Particle, Trajectory, Stack>::Cascade() {
-  //  kkk;
-  //  kk;
-}
-
-template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
-void Cascade::Init() {
-  fStack.Init();
-  fProcesseList.Init();
-}
-
-template <typename ProcessList, typename Particle, typename Trajectory, typename Stack>
-void Cascade::Run() {
-  if (!fStack.IsEmpty()) {
-    if (!fStack.IsEmpty()) {
-      Particle& p = fStack.GetNextParticle();
-      Step(p);
-    }
-    // do cascade equations, which can put new particles on Stack,
-    // thus, the double loop
-    // DoCascadeEquations(); //
-  }
-}
-
-template <typename Sequence, typename Trajectory>
-void Cascade::Step(Particle& particle) {
-  double nextStep = fProcesseList.MinStepLength(particle);
-  Trajectory trajectory = fProcesseList.Transport(particle, nextStep);
-  sequence.DoContinuous(particle, trajectory);
-  sequence.DoDiscrete(particle);
-}

Framework/Cascade/Cascade.h

-#ifndef _include_Cascade_h_
-#define _include_Cascade_h_
-
-#include <corsika/geometry/LineTrajectory.h> // to be removed
-#include <corsika/geometry/Point.h>          // to be removed
-#include <corsika/process/ProcessReturn.h>
-#include <corsika/units/PhysicalUnits.h>
-
-using namespace corsika::units::si;
-
-namespace corsika::cascade {
-
-  template <typename Trajectory, typename ProcessList, typename Stack>
-  class Cascade {
-
-    typedef typename Stack::ParticleType Particle;
-
-  public:
-    Cascade(ProcessList& pl, Stack& stack)
-        : fProcesseList(pl)
-        , fStack(stack) {}
-
-    void Init() {
-      fStack.Init();
-      fProcesseList.Init();
-    }
-
-    void Run() {
-      while (!fStack.IsEmpty()) {
-        while (!fStack.IsEmpty()) {
-          // Particle& p = *fStack.GetNextParticle();
-          EnergyType Emin;
-          typename Stack::StackIterator pMin(fStack, 0);
-          bool first = true;
-          for (typename Stack::StackIterator ip = fStack.begin(); ip != fStack.end();
-               ++ip) {
-            if (first || ip.GetEnergy() < Emin) {
-              first = false;
-              pMin = ip;
-              Emin = pMin.GetEnergy();
-            }
-          }
-
-          Step(pMin);
-        }
-        // do cascade equations, which can put new particles on Stack,
-        // thus, the double loop
-        // DoCascadeEquations(); //
-      }
-    }
-
-    void Step(Particle& particle) {
-      double nextStep = fProcesseList.MinStepLength(particle);
-      corsika::geometry::CoordinateSystem root;
-      Trajectory trajectory(
-          corsika::geometry::Point(root, {0_m, 0_m, 0_m}),
-          corsika::geometry::Vector<corsika::units::si::SpeedType::dimension_type>(
-              root, 0 * 1_m / second, 0 * 1_m / second, 1 * 1_m / second));
-      corsika::process::EProcessReturn status =
-          fProcesseList.DoContinuous(particle, trajectory, fStack);
-      if (status == corsika::process::EProcessReturn::eParticleAbsorbed) {
-        fStack.Delete(particle);
-      } else {
-        fProcesseList.DoDiscrete(particle, fStack);
-      }
-    }
-
-  private:
-    Stack& fStack;
-    ProcessList& fProcesseList;
-  };
-
-} // namespace corsika::cascade
-
-#endif

Framework/Cascade/Step.cc

-
-
-namespace cascade;
-
-void Cascade::Step(auto& sequence, Particle& particle) {
-  double nextStep = sequence.MinStepLength(particle);
-  Trajectory trajectory = sequence.Transport(particle, nextStep);
-  sequence.DoContinuous(particle, trajectory);
-  sequence.DoDiscrete(particle);
-}

Framework/Cascade/testCascade.cc

-#include <corsika/cascade/Cascade.h>
-#include <corsika/geometry/LineTrajectory.h>
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/stack/super_stupid/SuperStupidStack.h>
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-using namespace corsika::process;
-using namespace corsika::units;
-
-#include <iostream>
-using namespace std;
-
-static int fCount = 0;
-
-class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
-public:
-  ProcessSplit() {}
-
-  template <typename Particle>
-  double MinStepLength(Particle&) const {
-    return 0;
-  }
-
-  template <typename Particle, typename Trajectory, typename Stack>
-  EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
-    return EProcessReturn::eOk;
-  }
-
-  template <typename Particle, typename Stack>
-  void DoDiscrete(Particle& p, Stack& s) const {
-    EnergyType E = p.GetEnergy();
-    if (E < 85_MeV) {
-      p.Delete();
-      fCount++;
-    } else {
-      p.SetEnergy(E / 2);
-      s.NewParticle().SetEnergy(E / 2);
-    }
-  }
-
-  void Init() { fCount = 0; }
-
-  int GetCount() { return fCount; }
-
-private:
-};
-
-class ProcessReport : public corsika::process::BaseProcess<ProcessReport> {
-  bool fReport = false;
-
-public:
-  ProcessReport(bool v)
-      : fReport(v) {}
-
-  template <typename Particle>
-  double MinStepLength(Particle&) const {
-    return 0;
-  }
-
-  template <typename Particle, typename Trajectory, typename Stack>
-  EProcessReturn DoContinuous(Particle& p, Trajectory& t, Stack& s) const {
-    static int countStep = 0;
-    if (!fReport) return EProcessReturn::eOk;
-    // std::cout << "generation  " << countStep << std::endl;
-    int i = 0;
-    EnergyType Etot = 0_GeV;
-    for (auto& iterP : s) {
-      EnergyType E = iterP.GetEnergy();
-      Etot += E;
-      /*      std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()
-                << ", E=" << double(E / 1_GeV) << " GeV "
-                << " | " << std::endl;
-      */
-    }
-    countStep++;
-    // cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
-    cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << fCount
-         << endl;
-    return EProcessReturn::eOk;
-  }
-
-  template <typename Particle, typename Stack>
-  void DoDiscrete(Particle& p, Stack& s) const {}
-  void Init() {}
-};
-
-TEST_CASE("Cascade", "[Cascade]") {
-
-  ProcessReport p0(true);
-  ProcessSplit p1;
-  const auto sequence = p0 + p1;
-  corsika::stack::super_stupid::SuperStupidStack stack;
-
-  corsika::cascade::Cascade<corsika::geometry::LineTrajectory, decltype(sequence),
-                            decltype(stack)>
-      EAS(sequence, stack);
-
-  stack.Clear();
-  auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
-  particle.SetEnergy(E0);
-  EAS.Init();
-  EAS.Run();
-
-  SECTION("sectionTwo") {
-    for (int i = 0; i < 0; ++i) {
-      stack.Clear();
-      auto particle = stack.NewParticle();
-      EnergyType E0 = 100_GeV * pow(10, i);
-      particle.SetEnergy(E0);
-      EAS.Init();
-      EAS.Run();
-
-      // cout << "Result: E0=" << E0 / 1_GeV << "GeV, count=" << p1.GetCount() << endl;
-    }
-  }
-}

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 {
-
-  /*!
-   * Interface / base class for trajectories.
-   */
-  class BaseTrajectory {
-
-  public:
-    //!< 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
-
-#endif

Framework/Geometry/BaseVector.h

-#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) {}
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/CMakeLists.txt

-
-set (
-  GEOMETRY_SOURCES
-  CoordinateSystem.cc
-  )
-
-set (
-  GEOMETRY_HEADERS
-  Vector.h
-  Point.h
-  Sphere.h
-  Volume.h
-  CoordinateSystem.h
-  Helix.h
-  BaseVector.h
-  QuantityVector.h
-  BaseTrajectory.h
-  LineTrajectory.h
-  Trajectory.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
-  PRIVATE   ${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
-  )
-
-add_test (NAME testGeometry COMMAND testGeometry)
-

Framework/Geometry/CoordinateSystem.cc

-#include <corsika/geometry/CoordinateSystem.h>
-
-using namespace corsika::geometry;
-
-EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& c1,
-                                                   CoordinateSystem const& c2) {
-  CoordinateSystem const* a{&c1};
-  CoordinateSystem const* b{&c2};
-  CoordinateSystem const* commonBase{nullptr};
-
-  while (a != b && b != nullptr) {
-    a = &c1;
-
-    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::string("no connection between coordinate systems found!");
-  }
-
-  EigenTransform t = EigenTransform::Identity();
-
-  auto* p = &c1;
-
-  while (p != commonBase) {
-    t = p->GetTransform() * t;
-    p = p->GetReference();
-  }
-
-  p = &c2;
-
-  while (p != commonBase) {
-    t = p->GetTransform().inverse(Eigen::TransformTraits::Isometry) * t;
-    p = p->GetReference();
-  }
-
-  return t;
-}

Framework/Geometry/CoordinateSystem.h

-#ifndef _include_COORDINATESYSTEM_H_
-#define _include_COORDINATESYSTEM_H_
-
-#include <corsika/geometry/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <Eigen/Dense>
-
-typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
-typedef Eigen::Translation<double, 3> EigenTranslation;
-
-namespace corsika::geometry {
-
-  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) {}
-
-  public:
-    static EigenTransform GetTransformation(CoordinateSystem const& c1,
-                                            CoordinateSystem const& c2);
-
-    CoordinateSystem()
-        : // for creating the root CS
-        transf(EigenTransform::Identity()) {}
-
-    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 {
-      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) {
-      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/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 {
-  /*!
-   * 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::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 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
-
-#endif

Framework/Geometry/LineTrajectory.h

-#ifndef _include_LINETRAJECTORY_H
-#define _include_LINETRAJECTORY_H
-
-#include <corsika/geometry/BaseTrajectory.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-
-  class LineTrajectory : public BaseTrajectory {
-    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;
-
-    Point const r0;
-    VelocityVec const v0;
-
-  public:
-    LineTrajectory(Point const& pR0, VelocityVec const& pV0)
-        : r0(pR0)
-        , v0(pV0) {}
-
-    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
-
-#endif

Framework/Geometry/Point.h

-#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<phys::units::length_d> {
-  public:
-    Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector)
-        : BaseVector<phys::units::length_d>(pCS, pQVector) {}
-
-    Point(CoordinateSystem const& cs, LengthType x, LengthType y, LengthType z)
-        : BaseVector<phys::units::length_d>(cs, {x, y, z}) {}
-
-    auto GetCoordinates() const { return BaseVector<phys::units::length_d>::qVector; }
-
-    auto GetCoordinates(CoordinateSystem const& pCS) const {
-      if (&pCS == BaseVector<phys::units::length_d>::cs) {
-        return BaseVector<phys::units::length_d>::qVector;
-      } else {
-        return QuantityVector<phys::units::length_d>(
-            CoordinateSystem::GetTransformation(*BaseVector<phys::units::length_d>::cs,
-                                                pCS) *
-            BaseVector<phys::units::length_d>::qVector.eVector);
-      }
-    }
-
-    /*!
-     * transforms the Point into another CoordinateSystem by changing its
-     * coordinates interally
-     */
-    void rebase(CoordinateSystem const& pCS) {
-      BaseVector<phys::units::length_d>::qVector = GetCoordinates(pCS);
-      BaseVector<phys::units::length_d>::cs = &pCS;
-    }
-
-    Point operator+(Vector<phys::units::length_d> const& pVec) const {
-      return Point(
-          *BaseVector<phys::units::length_d>::cs,
-          GetCoordinates() + pVec.GetComponents(*BaseVector<phys::units::length_d>::cs));
-    }
-
-    /*!
-     * returns the distance Vector between two points
-     */
-    Vector<phys::units::length_d> operator-(Point const& pB) const {
-      auto& cs = *BaseVector<phys::units::length_d>::cs;
-      return Vector<phys::units::length_d>(cs, GetCoordinates() - pB.GetCoordinates(cs));
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/QuantityVector.h

-#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 {
-  protected:
-    // todo: check if we need to move "quantity" into namespace corsika::units
-    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; }
-  };
-
-} // namespace corsika::geometry
-
-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;
-}
-
-#endif

Framework/Geometry/Sphere.h

-#ifndef _include_SPHERE_H_
-#define _include_SPHERE_H_
-
-#include <corsika/geometry/Volume.h>
-#include <corsika/geometry/Point.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();
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

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::si::TimeType const fTStart, fTEnd;
-    BaseTrajectory const& fTrajectory;
-
-  public:
-    Trajectory(corsika::units::si::TimeType pTStart, corsika::units::si::TimeType pTEnd,
-               BaseTrajectory const& pTrajectory)
-        : fTStart(pTStart)
-        , fTEnd(pTEnd)
-        , fTrajectory(pTrajectory) {}
-
-    Point GetPosition(corsika::units::si::TimeType t) const {
-      return fTrajectory.GetPosition(t + fTStart);
-    }
-
-    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
-
-#endif

Framework/Geometry/Vector.h

-#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> {
-    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);
-      }
-    }
-  };
-
-} // namespace corsika::geometry
-
-#endif

Framework/Geometry/Volume.h

-#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

Framework/Geometry/testGeometry.cc

-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-#include <catch2/catch.hpp>
-
-#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 corsika::geometry;
-using namespace corsika::units::si;
-
-double constexpr absMargin = 1.0e-8;
-
-TEST_CASE("transformations between CoordinateSystems") {
-  CoordinateSystem rootCS;
-
-  REQUIRE(CoordinateSystem::GetTransformation(rootCS, rootCS)
-              .isApprox(EigenTransform::Identity()));
-
-  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};
-  Vector<magnetic_flux_density_d> v1(rootCS, components);
-
-  REQUIRE((p1.GetCoordinates() - coordinates).norm().magnitude() ==
-          Approx(0).margin(absMargin));
-  REQUIRE((p1.GetCoordinates(rootCS) - coordinates).norm().magnitude() ==
-          Approx(0).margin(absMargin));
-
-  SECTION("unconnected CoordinateSystems") {
-    CoordinateSystem rootCS2;
-    REQUIRE_THROWS(CoordinateSystem::GetTransformation(rootCS, rootCS2));
-  }
-
-  SECTION("translations") {
-    QuantityVector<length_d> const translationVector{0_m, 4_m, 0_m};
-
-    CoordinateSystem translatedCS = rootCS.translate(translationVector);
-
-    REQUIRE(translatedCS.GetReference() == &rootCS);
-
-    REQUIRE((p1.GetCoordinates(translatedCS) + translationVector).norm().magnitude() ==
-            Approx(0).margin(absMargin));
-
-    // Vectors are not subject to translations
-    REQUIRE(
-        (v1.GetComponents(rootCS) - v1.GetComponents(translatedCS)).norm().magnitude() ==
-        Approx(0).margin(absMargin));
-
-    Point p2(translatedCS, {0_m, 0_m, 0_m});
-    REQUIRE(((p2 - p1).GetComponents() - translationVector).norm().magnitude() ==
-            Approx(0).margin(absMargin));
-  }
-
-  SECTION("multiple translations") {
-    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};
-    CoordinateSystem cs3 = rootCS.translate(tv2);
-
-    QuantityVector<length_d> const tv3{0_m, 0_m, 2_m};
-    CoordinateSystem cs4 = cs3.translate(tv3);
-
-    REQUIRE(cs4.GetReference()->GetReference() == &rootCS);
-
-    REQUIRE(CoordinateSystem::GetTransformation(cs3, cs2).isApprox(
-        rootCS.translate({3_m, -5_m, 0_m}).GetTransform()));
-    REQUIRE(CoordinateSystem::GetTransformation(cs2, cs3).isApprox(
-        rootCS.translate({-3_m, +5_m, 0_m}).GetTransform()));
-  }
-
-  SECTION("rotations") {
-    QuantityVector<length_d> const axis{0_m, 0_m, 1_km};
-    double const angle = 90. / 180. * M_PI;
-
-    CoordinateSystem rotatedCS = rootCS.rotate(axis, angle);
-    REQUIRE(rotatedCS.GetReference() == &rootCS);
-
-    REQUIRE(v1.GetComponents(rotatedCS)[1].magnitude() ==
-            Approx((-1. * tesla).magnitude()));
-
-    // vector norm invariant under rotation
-    REQUIRE(v1.GetComponents(rotatedCS).norm().magnitude() ==
-            Approx(v1.GetComponents(rootCS).norm().magnitude()));
-  }
-
-  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};
-
-    double const angle = 90. / 180. * M_PI;
-
-    CoordinateSystem rotated1 = rootCS.rotate(zAxis, angle);
-    CoordinateSystem rotated2 = rotated1.rotate(yAxis, angle);
-    CoordinateSystem rotated3 = rotated2.rotate(zAxis, -angle);
-
-    CoordinateSystem combined = rootCS.rotate(xAxis, -angle);
-
-    auto comp1 = v1.GetComponents(rootCS);
-    auto comp3 = v1.GetComponents(combined);
-    REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
-  }
-}
-
-TEST_CASE("Sphere") {
-  CoordinateSystem rootCS;
-  Point center(rootCS, {0_m, 3_m, 4_m});
-  Sphere sphere(center, 5_m);
-
-  SECTION("isInside") {
-    REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
-    REQUIRE(sphere.Contains(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 const lineTrajectory(r0, v0);
-    CHECK((lineTrajectory.GetPosition(2_s).GetCoordinates() -
-           QuantityVector<length_d>(2_m, 0_m, 0_m))
-              .norm()
-              .magnitude() == Approx(0).margin(absMargin));
-
-    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});
-
-    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);
-
-    CHECK((helix.GetPosition(1_s).GetCoordinates() -
-           QuantityVector<length_d>(0_m, 0_m, 4_m))
-              .norm()
-              .magnitude() == Approx(0).margin(absMargin));
-
-    CHECK((helix.GetPosition(0.25_s).GetCoordinates() -
-           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));
-  }
-}