Framework/Geometry/Line.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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

Framework/Geometry/Plane.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::geometry {
-  class Plane {
-
-    using DimLessVec = Vector<corsika::units::si::dimensionless_d>;
-
-    Point const fCenter;
-    DimLessVec const fNormal;
-
-  public:
-    Plane(Point const& vCenter, DimLessVec const& vNormal)
-        : fCenter(vCenter)
-        , fNormal(vNormal.normalized()) {}
-
-    bool IsAbove(Point const& vP) const {
-      return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
-    }
-
-    units::si::LengthType DistanceTo(geometry::Point const& vP) const {
-      return (fNormal * (vP - fCenter).dot(fNormal)).norm();
-    }
-
-    Point const& GetCenter() const { return fCenter; }
-    DimLessVec const& GetNormal() const { return fNormal; }
-  };
-
-} // namespace corsika::geometry

Framework/Geometry/Point.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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; }
-    auto GetX() const { return BaseVector<length_d>::qVector.GetX(); }
-    auto GetY() const { return BaseVector<length_d>::qVector.GetY(); }
-    auto GetZ() const { return BaseVector<length_d>::qVector.GetZ(); }
-
-    /// 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

Framework/Geometry/QuantityVector.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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(double a, double b, double c)
-        : eVector{a, b, c} {
-      static_assert(
-          std::is_same_v<dim, phys::units::dimensionless_d>,
-          "initialization of dimensionful QuantityVector with pure numbers not allowed!");
-    }
-
-    QuantityVector(Eigen::Vector3d pBareVector)
-        : eVector(pBareVector) {}
-
-    auto operator[](size_t index) const {
-      return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
-    }
-
-    auto GetX() const { return (*this)[0]; }
-
-    auto GetY() const { return (*this)[1]; }
-
-    auto GetZ() const { return (*this)[2]; }
-
-    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 QuantityVector<dim>(eVector.normalized()); }
-
-    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

Framework/Geometry/RootCoordinateSystem.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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

Framework/Geometry/Sphere.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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

Framework/Geometry/Trajectory.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/geometry/Line.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 GetLength() const { return GetDistance(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 * corsika::units::si::second, t);
-    }
-
-    void LimitEndTo(corsika::units::si::LengthType limit) {
-      fTimeLength = T::TimeFromArclength(limit);
-    }
-
-    auto NormalizedDirection() const {
-      static_assert(std::is_same_v<T, corsika::geometry::Line>);
-      return T::GetV0().normalized();
-    }
-  };
-
-} // namespace corsika::geometry

Framework/Geometry/Vector.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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(); }
-    auto GetNorm() 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(); }
-    auto GetSquaredNorm() 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 ProdDim = phys::units::detail::product_d<dim, ScalarDim>;
-
-      return Vector<ProdDim>(*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 ProdDim = phys::units::detail::product_d<dim, dim2>;
-      return Vector<ProdDim>(*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 ProdDim = phys::units::detail::product_d<dim, dim2>;
-
-      return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
-                                                    bareResult);
-    }
-  };
-
-} // namespace corsika::geometry

Framework/Geometry/Volume.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#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

Framework/Geometry/testGeometry.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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 <catch2/catch.hpp>
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/geometry/Helix.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.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 =
-      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  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 = CoordinateSystem::CreateRootCS();
-    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};
-
-    QuantityVector<magnetic_flux_density_d> components{1. * tesla, 2. * tesla,
-                                                       3. * tesla};
-    Vector<magnetic_flux_density_d> v1(rootCS, components);
-
-    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(rotated3);
-    auto comp3 = v1.GetComponents(combined);
-    REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
-  }
-
-  SECTION("RotateToZ positive") {
-    Vector const v{rootCS, 0_m, 1_m, 1_m};
-    auto const csPrime = rootCS.RotateToZ(v);
-    Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
-    Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
-    Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
-
-    CHECK(xPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
-    CHECK(yPrime.dot(v).magnitude() == Approx(0).margin(absMargin));
-    CHECK((zPrime.dot(v) / 1_m).magnitude() == Approx(5 * sqrt(2)));
-
-    CHECK(zPrime.GetComponents(rootCS)[1].magnitude() ==
-          Approx(zPrime.GetComponents(rootCS)[2].magnitude()));
-    CHECK(zPrime.GetComponents(rootCS)[0].magnitude() == Approx(0));
-
-    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
-              yPrime.GetComponents(rootCS).eVector) == Approx(0));
-    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
-              xPrime.GetComponents(rootCS).eVector) == Approx(0));
-    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
-              zPrime.GetComponents(rootCS).eVector) == Approx(0));
-
-    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
-              yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
-    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
-              xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
-    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
-              zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
-  }
-
-  SECTION("RotateToZ negative") {
-    Vector const v{rootCS, 0_m, 0_m, -1_m};
-    auto const csPrime = rootCS.RotateToZ(v);
-    Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
-    Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
-    Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
-
-    CHECK(zPrime.dot(v).magnitude() > 0);
-    CHECK(xPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
-          Approx(0));
-    CHECK(yPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
-          Approx(0));
-
-    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
-              yPrime.GetComponents(rootCS).eVector) == Approx(0));
-    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
-              xPrime.GetComponents(rootCS).eVector) == Approx(0));
-    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
-              zPrime.GetComponents(rootCS).eVector) == Approx(0));
-
-    CHECK(yPrime.GetComponents(rootCS).eVector.dot(
-              yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
-    CHECK(xPrime.GetComponents(rootCS).eVector.dot(
-              xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
-    CHECK(zPrime.GetComponents(rootCS).eVector.dot(
-              zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
-  }
-}
-
-TEST_CASE("Sphere") {
-  CoordinateSystem& rootCS =
-      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  Point center(rootCS, {0_m, 3_m, 4_m});
-  Sphere sphere(center, 5_m);
-
-  SECTION("GetCenter") {
-    CHECK((sphere.GetCenter().GetCoordinates(rootCS) -
-           QuantityVector<length_d>(0_m, 3_m, 4_m))
-              .norm()
-              .magnitude() == Approx(0).margin(absMargin));
-    CHECK(sphere.GetRadius() / 5_m == Approx(1));
-  }
-
-  SECTION("Contains") {
-    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 =
-      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-  Point r0(rootCS, {0_m, 0_m, 0_m});
-
-  SECTION("Line") {
-    Vector<SpeedType::dimension_type> v0(rootCS,
-                                         {3_m / second, 0_m / second, 0_m / second});
-
-    Line const line(r0, v0);
-    CHECK(
-        (line.GetPosition(2_s).GetCoordinates() - QuantityVector<length_d>(6_m, 0_m, 0_m))
-            .norm()
-            .magnitude() == Approx(0).margin(absMargin));
-
-    CHECK((line.PositionFromArclength(4_m).GetCoordinates() -
-           QuantityVector<length_d>(4_m, 0_m, 0_m))
-              .norm()
-              .magnitude() == Approx(0).margin(absMargin));
-
-    CHECK((line.GetPosition(7_s) - line.PositionFromArclength(line.ArcLength(0_s, 7_s)))
-              .norm()
-              .magnitude() == Approx(0).margin(absMargin));
-
-    auto const t = 1_s;
-    Trajectory<Line> base(line, t);
-    CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
-
-    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(3));
-
-    CHECK((base.NormalizedDirection().GetComponents(rootCS) -
-           QuantityVector<dimensionless_d>{1, 0, 0})
-              .eVector.norm() == Approx(0).margin(absMargin));
-  }
-
-  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 T = 1_s;
-    auto const omegaC = 2 * M_PI / T;
-
-    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));
-
-    CHECK(
-        (helix.GetPosition(7_s) - helix.PositionFromArclength(helix.ArcLength(0_s, 7_s)))
-            .norm()
-            .magnitude() == Approx(0).margin(absMargin));
-
-    auto const t = 1234_s;
-    Trajectory<Helix> const base(helix, t);
-    CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
-
-    CHECK(base.ArcLength(0_s, 1_s) / 1_m == Approx(5));
-  }
-}

Framework/Logging/CMakeLists.txt

-# create the library
-add_library (CORSIKAlogging INTERFACE)
-
-# namespace of library -> location of header files
-set (
-  CORSIKAlogging_NAMESPACE
-  corsika/logging
-  )
-
-# header files of this library
-set (
-  CORSIKAlogging_HEADERS
-  Logging.h
-  )
-
-# copy the headers into the namespace
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAlogging
-  ${CORSIKAlogging_NAMESPACE}
-  ${CORSIKAlogging_HEADERS}
-  )
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAlogging
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-# and link against spdlog
-target_link_libraries(
-  CORSIKAlogging
-  INTERFACE
-  spdlog::spdlog
-)
-
-# install library
-install (
-  FILES ${CORSIKAlogging_HEADERS}
-  DESTINATION include/${CORSIKAlogging_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-
-CORSIKA_ADD_TEST (testLogging)
-target_link_libraries (
- testLogging
- CORSIKAlogging
- CORSIKAtesting
- )

Framework/Particles/CMakeLists.txt

-set(Python_ADDITIONAL_VERSIONS 3)
-find_package(PythonInterp 3 REQUIRED)
-
-add_custom_command (
-  OUTPUT  ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
-          ${PROJECT_BINARY_DIR}/Framework/Particles/particle_db.pkl
-  COMMAND ${PROJECT_SOURCE_DIR}/Framework/Particles/pdxml_reader.py 
-          ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleData.xml
-          ${PROJECT_SOURCE_DIR}/Framework/Particles/NuclearData.xml
-          ${PROJECT_SOURCE_DIR}/Framework/Particles/ParticleClassNames.xml
-  DEPENDS pdxml_reader.py
-          ParticleData.xml
-          NuclearData.xml
-          ParticleClassNames.xml
-  WORKING_DIRECTORY
-          ${PROJECT_BINARY_DIR}/Framework/Particles/
-  COMMENT "Read PYTHIA8 particle data and produce C++ source code GeneratedParticleProperties.inc"
-  VERBATIM
-  )
-
-set (
-  PARTICLE_SOURCES
-  ParticleProperties.cc
-  )
-
-# all public header files of library, includes automatic generated file(s)
-set (
-  PARTICLE_HEADERS
-  ParticleProperties.h
-  ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc # this one is auto-generated
-  )
-
-set (
-  PARTICLE_NAMESPACE
-  corsika/particles
-  )
-
-add_library (CORSIKAparticles STATIC ${PARTICLE_SOURCES})
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAparticles ${PARTICLE_NAMESPACE} ${PARTICLE_HEADERS})
-
-# ....................................................
-# since GeneratedParticleProperties.inc is an automatically produced file in the build directory,
-# create a symbolic link into the source tree, so that it can be found and edited more easily
-# this is not needed for the build to succeed! .......
-add_custom_command (
-  OUTPUT  ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc
-  COMMAND ${CMAKE_COMMAND} -E create_symlink ${PROJECT_BINARY_DIR}/include/corsika/particles/GeneratedParticleProperties.inc ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc
-  COMMENT "Generate link in source-dir: ${CMAKE_CURRENT_SOURCE_DIR}/GeneratedParticleProperties.inc"
-  )
-add_custom_target (SourceDirLink DEPENDS ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc)
-add_dependencies (CORSIKAparticles SourceDirLink)
-# .....................................................
-
-target_link_libraries (
-  CORSIKAparticles
-  PUBLIC
-  CORSIKAunits
-  CORSIKAlogging
-  )
-
-set_target_properties (
-  CORSIKAparticles
-  PROPERTIES
-  VERSION ${PROJECT_VERSION}
-  SOVERSION 1
-  PUBLIC_HEADER "${PARTICLE_HEADERS}"
-  )
-
-target_include_directories (
-  CORSIKAparticles
-  PUBLIC
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include>
-  )
-
-install (
-  FILES
-  ${PARTICLE_HEADERS}
-  DESTINATION
-  include/${PARTICLE_NAMESPACE}
-  )
-
-# --------------------
-# code unit testing
-CORSIKA_ADD_TEST(testParticles
-  SOURCES
-    testParticles.cc
-    ${PROJECT_BINARY_DIR}/Framework/Particles/GeneratedParticleProperties.inc
-)
-target_link_libraries (
-  testParticles
-  CORSIKAparticles
-  CORSIKAunits
-  CORSIKAtesting
-  )

Framework/Particles/ParticleProperties.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/particles/ParticleProperties.h>
-#include <iostream>
-
-namespace corsika::particles {
-
-  std::ostream& operator<<(std::ostream& stream, corsika::particles::Code const p) {
-    return stream << corsika::particles::GetName(p);
-  }
-
-  Code ConvertFromPDG(PDGCode p) {
-    static_assert(detail::conversionArray.size() % 2 == 1);
-    // this will fail, for the strange case where the maxPDG is negative...
-    unsigned int constexpr maxPDG{(detail::conversionArray.size() - 1) >> 1};
-    auto k = static_cast<PDGCodeType>(p);
-    if ((unsigned int)abs(k) <= maxPDG) {
-      return detail::conversionArray[k + maxPDG];
-    } else {
-      return detail::conversionMap.at(p);
-    }
-  }
-
-} // namespace corsika::particles

Framework/Particles/ParticleProperties.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-/**
-   @file Particles.h
-
-   Interface to particle properties
- */
-
-#pragma once
-
-#include <array>
-#include <cstdint>
-#include <iosfwd>
-
-#include <corsika/units/PhysicalUnits.h>
-
-/**
- *
- * The properties of all elementary particles is stored here. The data
- * are taken from the Pythia ParticleData.xml file.
- *
- */
-
-namespace corsika::particles {
-
-  /**
-   * @enum Code
-   * The Code enum is the actual place to define CORSIKA 8 particle codes.
-   */
-  enum class Code : int16_t;
-  enum class PDGCode : int32_t;
-  using CodeIntType = std::underlying_type<Code>::type;
-  using PDGCodeType = std::underlying_type<PDGCode>::type;
-
-  // forward declarations to be used in GeneratedParticleProperties
-  int16_t constexpr GetChargeNumber(Code const);
-  corsika::units::si::ElectricChargeType constexpr GetCharge(Code const);
-  corsika::units::si::HEPMassType constexpr GetMass(Code const);
-  PDGCode constexpr GetPDG(Code const);
-  constexpr std::string const& GetName(Code const);
-  corsika::units::si::TimeType constexpr GetLifetime(Code const);
-
-  bool constexpr IsNucleus(Code const);
-  bool constexpr IsHadron(Code const);
-  bool constexpr IsEM(Code const);
-  bool constexpr IsMuon(Code const);
-  bool constexpr IsNeutrino(Code const);
-  int constexpr GetNucleusA(Code const);
-  int constexpr GetNucleusZ(Code const);
-
-#include <corsika/particles/GeneratedParticleProperties.inc>
-
-  /*!
-   * returns mass of particle in natural units
-   */
-  corsika::units::si::HEPMassType constexpr GetMass(Code const p) {
-    if (p == Code::Nucleus)
-      throw std::runtime_error("Cannot GetMass() of particle::Nucleus -> unspecified");
-    return detail::masses[static_cast<CodeIntType>(p)];
-  }
-
-  /*!
-   * returns PDG id
-   */
-  PDGCode constexpr GetPDG(Code const p) {
-    return detail::pdg_codes[static_cast<CodeIntType>(p)];
-  }
-
-  /*!
-   * returns electric charge number of particle return 1 for a proton.
-   */
-  int16_t constexpr GetChargeNumber(Code const p) {
-    if (p == Code::Nucleus)
-      throw std::runtime_error(
-          "Cannot GetChargeNumber() of particle::Nucleus -> unspecified");
-    // electric_charges stores charges in units of (e/3), e.g. 3 for a proton
-    return detail::electric_charges[static_cast<CodeIntType>(p)] / 3;
-  }
-
-  /*!
-   * returns electric charge of particle, e.g. return 1.602e-19_C for a proton.
-   */
-  corsika::units::si::ElectricChargeType constexpr GetCharge(Code const p) {
-    if (p == Code::Nucleus)
-      throw std::runtime_error("Cannot GetCharge() of particle::Nucleus -> unspecified");
-    return GetChargeNumber(p) * corsika::units::constants::e;
-  }
-
-  constexpr std::string const& GetName(Code const p) {
-    return detail::names[static_cast<CodeIntType>(p)];
-  }
-
-  corsika::units::si::TimeType constexpr GetLifetime(Code const p) {
-    return detail::lifetime[static_cast<CodeIntType>(p)] * corsika::units::si::second;
-  }
-
-  bool constexpr IsHadron(Code const p) {
-    return detail::isHadron[static_cast<CodeIntType>(p)];
-  }
-
-  bool constexpr IsEM(Code c) {
-    return c == Code::Electron || c == Code::Positron || c == Code::Gamma;
-  }
-
-  bool constexpr IsMuon(Code c) { return c == Code::MuPlus || c == Code::MuMinus; }
-
-  bool constexpr IsNeutrino(Code c) {
-    return c == Code::NuE || c == Code::NuMu || c == Code::NuTau || c == Code::NuEBar ||
-           c == Code::NuMuBar || c == Code::NuTauBar;
-  }
-
-  int constexpr GetNucleusA(Code const p) {
-    if (p == Code::Nucleus) {
-      throw std::runtime_error("GetNucleusA(Code::Nucleus) is impossible!");
-    }
-    return detail::nucleusA[static_cast<CodeIntType>(p)];
-  }
-
-  int constexpr GetNucleusZ(Code const p) {
-    if (p == Code::Nucleus) {
-      throw std::runtime_error("GetNucleusZ(Code::Nucleus) is impossible!");
-    }
-    return detail::nucleusZ[static_cast<CodeIntType>(p)];
-  }
-
-  bool constexpr IsNucleus(Code const p) {
-    return (p == Code::Nucleus) || (GetNucleusA(p) != 0);
-  }
-
-  /**
-   * the output operator for humand-readable particle codes
-   **/
-
-  std::ostream& operator<<(std::ostream&, corsika::particles::Code);
-
-  Code ConvertFromPDG(PDGCode);
-
-  /**
-   * Get mass of nucleus
-   **/
-  corsika::units::si::HEPMassType constexpr GetNucleusMass(const int vA, const int vZ) {
-    return Proton::GetMass() * vZ + (vA - vZ) * Neutron::GetMass();
-  }
-
-  std::initializer_list<Code> constexpr getAllParticles() {
-    return detail::all_particles;
-  }
-
-} // namespace corsika::particles

Framework/Particles/particles.dox

-/**
-@page Particles Particle properties
-
-The properties of all particles are saved in static and flat
-arrays. There is a enum corsika::particles::Code to identify each
-particles, and each individual particles has its own static class,
-which can be used to retrieve its physical properties. 
-
-
-
-See namespace corsika::particles for all classes.
-
-
-
-*/
\ No newline at end of file

Framework/Particles/testParticles.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika::units;
-using namespace corsika::units::si;
-using namespace corsika::particles;
-
-TEST_CASE("ParticleProperties", "[Particles]") {
-
-  SECTION("Types") {
-    REQUIRE(Electron::GetCode() == Code::Electron);
-    REQUIRE(Positron::GetCode() == Code::Positron);
-    REQUIRE(Proton::GetCode() == Code::Proton);
-    REQUIRE(Neutron::GetCode() == Code::Neutron);
-    REQUIRE(Gamma::GetCode() == Code::Gamma);
-    REQUIRE(PiPlus::GetCode() == Code::PiPlus);
-  }
-
-  SECTION("Masses") {
-    REQUIRE(Electron::GetMass() / (511_keV) == Approx(1));
-    REQUIRE(Electron::GetMass() / GetMass(Code::Electron) == Approx(1));
-
-    REQUIRE((Proton::GetMass() + Neutron::GetMass()) /
-                corsika::units::constants::nucleonMass ==
-            Approx(2));
-  }
-
-  SECTION("Charges") {
-    REQUIRE(Electron::GetCharge() / constants::e == Approx(-1));
-    REQUIRE(Positron::GetCharge() / constants::e == Approx(+1));
-    REQUIRE(GetCharge(Positron::GetAntiParticle()) / constants::e == Approx(-1));
-  }
-
-  SECTION("Names") {
-    REQUIRE(Electron::GetName() == "e-");
-    REQUIRE(PiMinus::GetName() == "pi-");
-    REQUIRE(Nucleus::GetName() == "nucleus");
-    REQUIRE(Iron::GetName() == "iron");
-  }
-
-  SECTION("PDG") {
-    REQUIRE(GetPDG(Code::PiPlus) == PDGCode::PiPlus);
-    REQUIRE(GetPDG(Code::DPlus) == PDGCode::DPlus);
-    REQUIRE(GetPDG(Code::NuMu) == PDGCode::NuMu);
-    REQUIRE(GetPDG(Code::NuE) == PDGCode::NuE);
-    REQUIRE(GetPDG(Code::MuMinus) == PDGCode::MuMinus);
-
-    REQUIRE(static_cast<int>(GetPDG(Code::PiPlus)) == 211);
-    REQUIRE(static_cast<int>(GetPDG(Code::DPlus)) == 411);
-    REQUIRE(static_cast<int>(GetPDG(Code::NuMu)) == 14);
-    REQUIRE(static_cast<int>(GetPDG(Code::NuEBar)) == -12);
-    REQUIRE(static_cast<int>(GetPDG(Code::MuMinus)) == 13);
-  }
-
-  SECTION("Conversion PDG -> internal") {
-    REQUIRE(ConvertFromPDG(PDGCode::KStarMinus) == Code::KStarMinus);
-    REQUIRE(ConvertFromPDG(PDGCode::MuPlus) == Code::MuPlus);
-    REQUIRE(ConvertFromPDG(PDGCode::SigmaStarCMinusBar) == Code::SigmaStarCMinusBar);
-  }
-
-  SECTION("Lifetimes") {
-    REQUIRE(GetLifetime(Code::Electron) ==
-            std::numeric_limits<double>::infinity() * corsika::units::si::second);
-    REQUIRE(GetLifetime(Code::DPlus) < GetLifetime(Code::Gamma));
-    REQUIRE(GetLifetime(Code::RhoPlus) / corsika::units::si::second ==
-            (Approx(4.414566727909413e-24).epsilon(1e-3)));
-    REQUIRE(GetLifetime(Code::SigmaMinusBar) / corsika::units::si::second ==
-            (Approx(8.018880848563575e-11).epsilon(1e-5)));
-    REQUIRE(GetLifetime(Code::MuPlus) / corsika::units::si::second ==
-            (Approx(2.1970332555864364e-06).epsilon(1e-5)));
-  }
-
-  SECTION("Particle groups: electromagnetic") {
-    REQUIRE(IsEM(Code::Gamma));
-    REQUIRE(IsEM(Code::Electron));
-    REQUIRE_FALSE(IsEM(Code::MuPlus));
-    REQUIRE_FALSE(IsEM(Code::NuE));
-    REQUIRE_FALSE(IsEM(Code::Proton));
-    REQUIRE_FALSE(IsEM(Code::PiPlus));
-    REQUIRE_FALSE(IsEM(Code::Oxygen));
-  }
-
-  SECTION("Particle groups: hadrons") {
-    REQUIRE_FALSE(IsHadron(Code::Gamma));
-    REQUIRE_FALSE(IsHadron(Code::Electron));
-    REQUIRE_FALSE(IsHadron(Code::MuPlus));
-    REQUIRE_FALSE(IsHadron(Code::NuE));
-    REQUIRE(IsHadron(Code::Proton));
-    REQUIRE(IsHadron(Code::PiPlus));
-    REQUIRE(IsHadron(Code::Oxygen));
-    REQUIRE(IsHadron(Code::Nucleus));
-  }
-
-  SECTION("Particle groups: muons") {
-    REQUIRE_FALSE(IsMuon(Code::Gamma));
-    REQUIRE_FALSE(IsMuon(Code::Electron));
-    REQUIRE(IsMuon(Code::MuPlus));
-    REQUIRE_FALSE(IsMuon(Code::NuE));
-    REQUIRE_FALSE(IsMuon(Code::Proton));
-    REQUIRE_FALSE(IsMuon(Code::PiPlus));
-    REQUIRE_FALSE(IsMuon(Code::Oxygen));
-  }
-
-  SECTION("Particle groups: neutrinos") {
-    REQUIRE_FALSE(IsNeutrino(Code::Gamma));
-    REQUIRE_FALSE(IsNeutrino(Code::Electron));
-    REQUIRE_FALSE(IsNeutrino(Code::MuPlus));
-    REQUIRE(IsNeutrino(Code::NuE));
-    REQUIRE_FALSE(IsNeutrino(Code::Proton));
-    REQUIRE_FALSE(IsNeutrino(Code::PiPlus));
-    REQUIRE_FALSE(IsNeutrino(Code::Oxygen));
-  }
-
-  SECTION("Nuclei") {
-    REQUIRE_FALSE(IsNucleus(Code::Gamma));
-    REQUIRE(IsNucleus(Code::Argon));
-    REQUIRE_FALSE(IsNucleus(Code::Proton));
-    REQUIRE(IsNucleus(Code::Hydrogen));
-    REQUIRE(Argon::IsNucleus());
-    REQUIRE_FALSE(EtaC::IsNucleus());
-
-    REQUIRE(GetNucleusA(Code::Hydrogen) == 1);
-    REQUIRE(GetNucleusA(Code::Tritium) == 3);
-    REQUIRE(Hydrogen::GetNucleusZ() == 1);
-    REQUIRE(Tritium::GetNucleusA() == 3);
-
-    REQUIRE_THROWS(GetNucleusA(Code::Nucleus));
-    REQUIRE_THROWS(GetNucleusZ(Code::Nucleus));
-  }
-}

Framework/ProcessSequence/BaseProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/process/ProcessReturn.h> // for convenience
-
-namespace corsika::process {
-
-  class TDerived; // fwd decl
-
-  /**
-     \class BaseProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type BaseProcess<T>
-
-     \todo rename BaseProcess into just Process
-   */
-  class _BaseProcess {};
-
-  template <typename TDerived>
-  class BaseProcess : _BaseProcess {
-  protected:
-    friend TDerived;
-
-    BaseProcess() = default; // protected constructor will allow only
-                             // derived classes to be created, not
-                             // BaseProcess itself
-
-    TDerived& GetRef() { return static_cast<TDerived&>(*this); }
-    const TDerived& GetRef() const { return static_cast<const TDerived&>(*this); }
-
-  public:
-    // Base processor type for use in other template classes
-    using TProcessType = TDerived;
-  };
-
-} // namespace corsika::process

Framework/ProcessSequence/BoundaryCrossingProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/process/ProcessReturn.h>
-
-#include <type_traits>
-
-namespace corsika::process {
-
-  template <typename TDerived>
-  class BoundaryCrossingProcess : public BaseProcess<TDerived> {
-  private:
-  protected:
-  public:
-    /**
-     * This method is called when a particle crosses the boundary between the nodes
-     * \p from and \p to.
-     */
-    template <typename TParticle, typename TVolumeNode>
-    EProcessReturn DoBoundaryCrossing(TParticle&, TVolumeNode const& from,
-                                      TVolumeNode const& to);
-  };
-
-} // namespace corsika::process

Framework/ProcessSequence/CMakeLists.txt

-add_library (CORSIKAprocesssequence INTERFACE)
-
-#namespace of library->location of header files
-set (
-  CORSIKAprocesssequence_NAMESPACE
-  corsika/process
-  )
-
-#header files of this library
-set (
-  CORSIKAprocesssequence_HEADERS
-  BaseProcess.h
-  BoundaryCrossingProcess.h
-  ContinuousProcess.h
-  SecondariesProcess.h
-  InteractionProcess.h
-  StackProcess.h
-  DecayProcess.h
-  ProcessSequence.h
-  SwitchProcessSequence.h
-  ProcessReturn.h
-  ProcessTraits.h
-  )
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAprocesssequence ${CORSIKAprocesssequence_NAMESPACE} ${CORSIKAprocesssequence_HEADERS})
-
-#include directive for upstream code
-target_link_libraries (
-  CORSIKAprocesssequence
-  INTERFACE
-  CORSIKAsetup
-  CORSIKAlogging
-  )
-
-target_include_directories (
-  CORSIKAprocesssequence
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  )
-
-#install library
-install (
-  FILES ${CORSIKAprocesssequence_HEADERS}
-  DESTINATION include/${CORSIKAprocesssequence_NAMESPACE}
-  )
-  
-target_link_libraries (
-  CORSIKAprocesssequence
-  INTERFACE
-  CORSIKAenvironment
-)
-
-#-- -- -- -- -- -- -- --
-#code unit testing
-CORSIKA_ADD_TEST (testProcessSequence)
-target_link_libraries (
-  testProcessSequence
-  CORSIKAgeometry
-  CORSIKAprocesssequence
-  CORSIKAtesting
-  )
-
-# # CORSIKA_ADD_TEST (testSwitchProcessSequence)
-# # target_link_libraries (
-# #   testSwitchProcessSequence
-# #   CORSIKAgeometry
-# #   CORSIKAprocesssequence
-# #   CORSIKAtesting
-# #   )

Framework/ProcessSequence/ContinuousProcess.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * 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.
- */
-
-#pragma once
-
-#include <corsika/process/BaseProcess.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::process {
-
-  /**
-     \class ContinuousProcess
-
-     The structural base type of a process object in a
-     ProcessSequence. Both, the ProcessSequence and all its elements
-     are of type ContinuousProcess<T>
-
-   */
-
-  template <typename TDerived>
-  class ContinuousProcess : public BaseProcess<TDerived> {
-  private:
-  protected:
-  public:
-    // here starts the interface part
-    // -> enforce TDerived to implement DoContinuous...
-    template <typename TParticle, typename TTrack>
-    EProcessReturn DoContinuous(TParticle&, TTrack const&) const;
-
-    // -> enforce TDerived to implement MaxStepLength...
-    template <typename TParticle, typename TTrack>
-    units::si::LengthType MaxStepLength(TParticle const& p, TTrack const& track) const;
-  };
-
-} // namespace corsika::process