Environment/ShowerAxis.h

-/*
- * (c) Copyright 2020 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/environment/Environment.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cassert>
-#include <cstdlib>
-#include <fstream>
-#include <functional>
-#include <iterator>
-#include <memory>
-#include <stdexcept>
-#include <vector>
-
-#include <iostream>
-
-#include <boost/math/quadrature/gauss_kronrod.hpp>
-
-namespace corsika::environment {
-
-  /**
-   * \class ShowerAxis
-   *
-   * The environment::ShowerAxis is created from a geometry::Point and
-   * a geometry::Vector and inside an Environment. It internally uses
-   * a table with steps=10000 (default) rows for interpolation.
-   *
-   * The shower axis can convert location in the shower into a
-   * projected grammage along the shower axis.
-   *
-   **/
-
-  class ShowerAxis {
-  public:
-    template <typename TEnvModel>
-    ShowerAxis(geometry::Point const& pStart,
-               geometry::Vector<units::si::length_d> length,
-               environment::Environment<TEnvModel> const& env, int steps = 10'000)
-        : pointStart_(pStart)
-        , length_(length)
-        , max_length_(length_.norm())
-        , steplength_(max_length_ / steps)
-        , axis_normalized_(length / max_length_)
-        , X_(steps + 1) {
-      auto const* const universe = env.GetUniverse().get();
-
-      auto rho = [pStart, length, universe](double x) {
-        auto const p = pStart + length * x;
-        auto const* node = universe->GetContainingNode(p);
-        return node->GetModelProperties().GetMassDensity(p).magnitude();
-      };
-
-      double error;
-      int k = 0;
-      X_[0] = units::si::GrammageType::zero();
-      auto sum = units::si::GrammageType::zero();
-
-      for (int i = 1; i <= steps; ++i) {
-        auto const x_prev = (i - 1.) / steps;
-        auto const d_prev = max_length_ * x_prev;
-        auto const x = double(i) / steps;
-        auto const r = boost::math::quadrature::gauss_kronrod<double, 15>::integrate(
-            rho, x_prev, x, 15, 1e-9, &error);
-        auto const result =
-            units::si::MassDensityType(phys::units::detail::magnitude_tag, r) *
-            max_length_;
-
-        sum += result;
-        X_[i] = sum;
-
-        for (; sum > k * X_binning_; ++k) {
-          d_.emplace_back(d_prev + k * X_binning_ * steplength_ / result);
-        }
-      }
-
-      assert(std::is_sorted(X_.cbegin(), X_.cend()));
-      assert(std::is_sorted(d_.cbegin(), d_.cend()));
-    }
-
-    units::si::LengthType steplength() const;
-
-    units::si::GrammageType maximumX() const;
-
-    units::si::GrammageType minimumX() const;
-
-    units::si::GrammageType projectedX(geometry::Point const& p) const;
-
-    units::si::GrammageType X(units::si::LengthType) const;
-
-    geometry::Vector<units::si::dimensionless_d> const& GetDirection() const;
-
-    geometry::Point const& GetStart() const;
-
-  private:
-    geometry::Point const pointStart_;
-    geometry::Vector<units::si::length_d> const length_;
-    units::si::LengthType const max_length_, steplength_;
-    geometry::Vector<units::si::dimensionless_d> const axis_normalized_;
-    std::vector<units::si::GrammageType> X_;
-
-    // for storing the lengths corresponding to equidistant X values
-    units::si::GrammageType const X_binning_ = std::invoke([]() {
-      using namespace units::si;
-      return 1_g / 1_cm / 1_cm;
-    });
-    std::vector<units::si::LengthType> d_;
-  };
-} // namespace corsika::environment

Environment/SlidingPlanarExponential.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/environment/FlatExponential.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  // clang-format off
-  /**
-   * The SlidingPlanarExponential models mass density as
-   * \f[
-   *   \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right).
-   * \f]
-   * For grammage/length conversion, the density distribution is approximated as
-   * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the axis pointing
-   * from \f$ p_0 \f$ to \f$ r_0 \f$.
-   */
-   //clang-format on
-   
-  template <class T>
-  class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
-                                   public T {
-    NuclearComposition const nuclComp_;
-    units::si::LengthType const referenceHeight_;
-
-    using Base = BaseExponential<SlidingPlanarExponential<T>>;
-
-  public:
-    SlidingPlanarExponential(geometry::Point const& p0, units::si::MassDensityType rho0,
-                             units::si::LengthType lambda, NuclearComposition nuclComp, units::si::LengthType referenceHeight = units::si::LengthType::zero())
-        : Base(p0, rho0, lambda)
-        , nuclComp_(nuclComp),
-        referenceHeight_(referenceHeight) {}
-
-    units::si::MassDensityType GetMassDensity(
-        geometry::Point const& p) const override {
-      auto const height = (p - Base::fP0).norm() - referenceHeight_;
-      return Base::fRho0 * exp(Base::fInvLambda * height);
-    }
-
-    NuclearComposition const& GetNuclearComposition() const override { return nuclComp_; }
-
-    units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& line,
-        units::si::LengthType l) const override {
-      auto const axis = (line.GetR0() - Base::fP0).normalized();
-      return Base::IntegratedGrammage(line, l, axis);
-    }
-
-    units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& line,
-        units::si::GrammageType grammage) const override {
-      auto const axis = (line.GetR0() - Base::fP0).normalized();
-      return Base::ArclengthFromGrammage(line, grammage, axis);
-    }
-  };
-
-} // namespace corsika::environment

Environment/VolumeTreeNode.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/environment/IMediumModel.h>
-#include <corsika/geometry/Volume.h>
-#include <memory>
-#include <vector>
-
-namespace corsika::environment {
-
-  class Empty {}; //<! intended for usage as default template argument
-
-  template <typename TModelProperties = Empty>
-  class VolumeTreeNode {
-  public:
-    using IModelProperties = TModelProperties;
-    using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
-    using IMPSharedPtr = std::shared_ptr<IModelProperties>;
-    using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
-
-    VolumeTreeNode(VolUPtr pVolume = nullptr)
-        : fGeoVolume(std::move(pVolume)) {}
-
-    //! convenience function equivalent to Volume::Contains
-    bool Contains(corsika::geometry::Point const& p) const {
-      return fGeoVolume->Contains(p);
-    }
-
-    VolumeTreeNode<IModelProperties> const* Excludes(
-        corsika::geometry::Point const& p) const {
-      auto exclContainsIter =
-          std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
-                       [&](auto const& s) { return bool(s->Contains(p)); });
-
-      return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
-    }
-
-    /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
-     * \class Point \p p, or nullptr iff \p p is not contained in this volume.
-     */
-    VolumeTreeNode<IModelProperties> const* GetContainingNode(
-        corsika::geometry::Point const& p) const {
-      if (!Contains(p)) { return nullptr; }
-
-      if (auto const childContainsIter =
-              std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
-                           [&](auto const& s) { return bool(s->Contains(p)); });
-          childContainsIter == fChildNodes.cend()) // not contained in any of the children
-      {
-        if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
-        {
-          return exclContainsIter->GetContainingNode(p);
-        } else {
-          return this;
-        }
-      } else {
-        return (*childContainsIter)->GetContainingNode(p);
-      }
-    }
-
-    /**
-     * Traverses the VolumeTree pre- or post-order and calls the functor  \p func for each
-     * node. \p func takes a reference to VolumeTreeNode as argument. The return value \p
-     * func is ignored.
-     */
-    template <typename TCallable, bool preorder = true>
-    void walk(TCallable func) {
-      if constexpr (preorder) { func(*this); }
-
-      std::for_each(fChildNodes.begin(), fChildNodes.end(),
-                    [&](auto& v) { v->walk(func); });
-
-      if constexpr (!preorder) { func(*this); };
-    }
-
-    void AddChild(VTNUPtr pChild) {
-      pChild->fParentNode = this;
-      fChildNodes.push_back(std::move(pChild));
-      // It is a bad idea to return an iterator to the inserted element
-      // because it might get invalidated when the vector needs to grow
-      // later and the caller won't notice.
-    }
-
-    void ExcludeOverlapWith(VTNUPtr const& pNode) {
-      fExcludedNodes.push_back(pNode.get());
-    }
-
-    auto* GetParent() const { return fParentNode; };
-
-    auto const& GetChildNodes() const { return fChildNodes; }
-
-    auto const& GetExcludedNodes() const { return fExcludedNodes; }
-
-    auto const& GetVolume() const { return *fGeoVolume; }
-
-    auto const& GetModelProperties() const { return *fModelProperties; }
-
-    bool HasModelProperties() const { return fModelProperties.get() != nullptr; }
-
-    template <typename ModelProperties, typename... Args>
-    auto SetModelProperties(Args&&... args) {
-      static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
-                    "unusable type provided");
-
-      fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
-      return fModelProperties;
-    }
-
-    void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; }
-
-    template <class MediumType, typename... Args>
-    static auto CreateMedium(Args&&... args) {
-      static_assert(std::is_base_of_v<IMediumModel, MediumType>,
-                    "unusable type provided, needs to be derived from \"IMediumModel\"");
-
-      return std::make_shared<MediumType>(std::forward<Args>(args)...);
-    }
-
-  private:
-    std::vector<VTNUPtr> fChildNodes;
-    std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
-    VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
-    VolUPtr fGeoVolume;
-    IMPSharedPtr fModelProperties;
-  };
-
-} // namespace corsika::environment

Environment/testEnvironment.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/environment/DensityFunction.h>
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/IMagneticFieldModel.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/IRefractiveIndexModel.h>
-#include <corsika/environment/InhomogeneousMedium.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
-#include <corsika/environment/LinearApproximationIntegrator.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/SlidingPlanarExponential.h>
-#include <corsika/environment/UniformMagneticField.h>
-#include <corsika/environment/UniformRefractiveIndex.h>
-#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika::geometry;
-using namespace corsika::environment;
-using namespace corsika::particles;
-using namespace corsika::units::si;
-using namespace corsika;
-
-CoordinateSystem const& gCS =
-    RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-Point const gOrigin(gCS, {0_m, 0_m, 0_m});
-
-TEST_CASE("HomogeneousMedium") {
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-  HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
-}
-
-TEST_CASE("FlatExponential") {
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-
-  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
-  LengthType const lambda = 3_m;
-  auto const rho0 = 1_g / units::static_pow<3>(1_cm);
-  FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
-                                             protonComposition);
-  auto const tEnd = 5_s;
-
-  SECTION("horizontal") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {20_cm / second, 0_m / second, 0_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-
-    REQUIRE((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
-    REQUIRE((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
-  }
-
-  SECTION("vertical") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 0_m / second, 5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-    LengthType const length = 2 * lambda;
-    GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
-
-    REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
-    REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
-  }
-
-  SECTION("escape grammage") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 0_m / second, -5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-
-    GrammageType const escapeGrammage = rho0 * lambda;
-
-    REQUIRE(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
-    REQUIRE(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
-            std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
-  }
-
-  SECTION("inclined") {
-    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                                 gCS, {0_m / second, 5_m / second, 5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-    double const cosTheta = M_SQRT1_2;
-    LengthType const length = 2 * lambda;
-    GrammageType const exact =
-        rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
-    REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
-    REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
-  }
-}
-
-TEST_CASE("SlidingPlanarExponential") {
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-
-  LengthType const lambda = 3_m;
-  auto const rho0 = 1_g / units::static_pow<3>(1_cm);
-  auto const tEnd = 5_s;
-
-  SlidingPlanarExponential<IMediumModel> const medium(gOrigin, rho0, lambda,
-                                                      protonComposition);
-
-  SECTION("density") {
-    CHECK(medium.GetMassDensity({gCS, {0_m, 0_m, 3_m}}) /
-              medium.GetMassDensity({gCS, {0_m, 3_m, 0_m}}) ==
-          Approx(1));
-  }
-
-  SECTION("vertical") {
-    Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
-    FlatExponential<IMediumModel> const flat(gOrigin, axis, rho0, lambda,
-                                             protonComposition);
-    Line const line({gCS, {0_m, 0_m, 1_m}},
-                    Vector<SpeedType::dimension_type>(
-                        gCS, {0_m / second, 0_m / second, 5_m / second}));
-    Trajectory<Line> const trajectory(line, tEnd);
-
-    CHECK(medium.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
-          flat.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
-    CHECK(medium.IntegratedGrammage(trajectory, 2_m).magnitude() ==
-          flat.IntegratedGrammage(trajectory, 2_m).magnitude());
-    CHECK(medium.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
-          flat.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
-  }
-}
-
-auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
-
-struct Exponential {
-  auto operator()(corsika::geometry::Point const& p) const {
-    return exp(p.GetCoordinates()[0] / 1_m) * rho0;
-  }
-
-  template <int N>
-  auto Derivative(Point const& p, Vector<dimensionless_d> const& v) const {
-    return v.GetComponents()[0] * (*this)(p) / corsika::units::static_pow<N>(1_m);
-  }
-
-  auto FirstDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
-    return Derivative<1>(p, v);
-  }
-
-  auto SecondDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
-    return Derivative<2>(p, v);
-  }
-};
-
-TEST_CASE("InhomogeneousMedium") {
-  Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
-
-  Line line(gOrigin, Vector<SpeedType::dimension_type>(
-                         gCS, {20_m / second, 0_m / second, 0_m / second}));
-
-  auto const tEnd = 5_s;
-  Trajectory<Line> const trajectory(line, tEnd);
-
-  Exponential const e;
-  DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
-
-  SECTION("DensityFunction") {
-    REQUIRE(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
-            Approx(1));
-    REQUIRE(rho.EvaluateAt(gOrigin) == e(gOrigin));
-  }
-
-  auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
-  auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); };
-
-  auto constexpr l = 15_cm;
-
-  NuclearComposition const composition{{Code::Proton}, {1.f}};
-  InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
-
-  SECTION("Integration") {
-    REQUIRE(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
-            Approx(1).epsilon(1e-2));
-    REQUIRE(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
-                exactLength(exactGrammage(l)) ==
-            Approx(1).epsilon(1e-2));
-    REQUIRE(rho.MaximumLength(trajectory, 1e-2) >
-            l); // todo: write reasonable test when implementation is working
-
-    REQUIRE(rho.IntegrateGrammage(trajectory, l) ==
-            inhMedium.IntegratedGrammage(trajectory, l));
-    REQUIRE(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
-            inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
-  }
-}
-
-TEST_CASE("LayeredSphericalAtmosphereBuilder") {
-  LayeredSphericalAtmosphereBuilder builder(gOrigin);
-  builder.setNuclearComposition(
-      {{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
-
-  builder.addLinearLayer(1_km, 10_km);
-  builder.addLinearLayer(2_km, 20_km);
-  builder.addLinearLayer(3_km, 30_km);
-
-  REQUIRE(builder.size() == 3);
-
-  auto const builtEnv = builder.assemble();
-  auto const& univ = builtEnv.GetUniverse();
-
-  REQUIRE(builder.size() == 0);
-
-  auto const R = builder.getEarthRadius();
-
-  REQUIRE(univ->GetChildNodes().size() == 1);
-
-  REQUIRE(univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 35_km)) == univ.get());
-  REQUIRE(dynamic_cast<Sphere const&>(
-              univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 8_km))->GetVolume())
-              .GetRadius() == R + 10_km);
-  REQUIRE(dynamic_cast<Sphere const&>(
-              univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 12_km))->GetVolume())
-              .GetRadius() == R + 20_km);
-  REQUIRE(dynamic_cast<Sphere const&>(
-              univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 24_km))->GetVolume())
-              .GetRadius() == R + 30_km);
-}
-
-TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
-
-  // setup our interface types
-  using IModelInterface = IMagneticFieldModel<IMediumModel>;
-  using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;
-
-  // the composition we use for the homogenous medium
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-
-  // create a magnetic field vector
-  Vector B0(gCS, 0_T, 0_T, 0_T);
-
-  // the constant density
-  const auto density{19.2_g / cube(1_cm)};
-
-  // create our atmospheric model
-  AtmModel medium(B0, density, protonComposition);
-
-  // and test at several locations
-  REQUIRE(B0.GetComponents(gCS) ==
-          medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
-  REQUIRE(
-      B0.GetComponents(gCS) ==
-      medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).GetComponents(gCS));
-  REQUIRE(B0.GetComponents(gCS) ==
-          medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
-
-  // create a new magnetic field vector
-  Vector B1(gCS, 23_T, 57_T, -4_T);
-
-  // and update this atmospheric model
-  medium.SetMagneticField(B1);
-
-  // and test at several locations
-  REQUIRE(B1.GetComponents(gCS) ==
-          medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
-  REQUIRE(
-      B1.GetComponents(gCS) ==
-      medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).GetComponents(gCS));
-  REQUIRE(B1.GetComponents(gCS) ==
-          medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
-
-  // check the density and nuclear composition
-  REQUIRE(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
-  medium.GetNuclearComposition();
-
-  // create a line of length 1 m
-  Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                               gCS, {1_m / second, 0_m / second, 0_m / second}));
-
-  // the end time of our line
-  auto const tEnd = 1_s;
-
-  // and the associated trajectory
-  Trajectory<Line> const trajectory(line, tEnd);
-
-  // and check the integrated grammage
-  REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
-  REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
-}
-
-TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
-
-  // setup our interface types
-  using IModelInterface = IRefractiveIndexModel<IMediumModel>;
-  using AtmModel = UniformRefractiveIndex<HomogeneousMedium<IModelInterface>>;
-
-  // the constant density
-  const auto density{19.2_g / cube(1_cm)};
-
-  // the composition we use for the homogenous medium
-  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
-                                             std::vector<float>{1.f});
-
-  // the refrative index that we use
-  const double n{1.000327};
-
-  // create the atmospheric model
-  AtmModel medium(n, density, protonComposition);
-
-  // and require that it is constant
-  REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
-  REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
-  REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
-  REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
-
-  // a new refractive index
-  const double n2{2.3472123};
-
-  // update the refractive index of this atmospheric model
-  medium.SetRefractiveIndex(n2);
-
-  // check that the returned refractive index is correct
-  REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
-  REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
-  REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
-  REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
-
-  // define our axis vector
-  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
-
-  // check the density and nuclear composition
-  REQUIRE(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
-  medium.GetNuclearComposition();
-
-  // create a line of length 1 m
-  Line const line(gOrigin, Vector<SpeedType::dimension_type>(
-                               gCS, {1_m / second, 0_m / second, 0_m / second}));
-
-  // the end time of our line
-  auto const tEnd = 1_s;
-
-  // and the associated trajectory
-  Trajectory<Line> const trajectory(line, tEnd);
-
-  // and check the integrated grammage
-  REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
-  REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
-}

Environment/testShowerAxis.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/environment/DensityFunction.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/ShowerAxis.h>
-#include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika::geometry;
-using namespace corsika::environment;
-using namespace corsika::particles;
-using namespace corsika::units;
-using namespace corsika::units::si;
-using namespace corsika;
-
-const auto density = 1_kg / (1_m * 1_m * 1_m);
-
-auto setupEnvironment(particles::Code vTargetCode) {
-  // setup environment, geometry
-  auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
-  auto& universe = *(env->GetUniverse());
-  const geometry::CoordinateSystem& cs = env->GetCoordinateSystem();
-
-  auto theMedium =
-      environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
-          geometry::Point{cs, 0_m, 0_m, 0_m},
-          1_km * std::numeric_limits<double>::infinity());
-
-  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      density, environment::NuclearComposition(std::vector<particles::Code>{vTargetCode},
-                                               std::vector<float>{1.}));
-
-  auto const* nodePtr = theMedium.get();
-  universe.AddChild(std::move(theMedium));
-
-  return std::make_tuple(std::move(env), &cs, nodePtr);
-}
-
-TEST_CASE("Homogeneous Density") {
-  auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Nitrogen);
-  auto const& cs = *csPtr;
-  [[maybe_unused]] auto const& env_dummy = env;
-  [[maybe_unused]] auto const& node_dummy = nodePtr;
-
-  auto const observationHeight = 0_km;
-  auto const injectionHeight = 10_km;
-  auto const t = -observationHeight + injectionHeight;
-  Point const showerCore{cs, 0_m, 0_m, observationHeight};
-  Point const injectionPos = showerCore + Vector<dimensionless_d>{cs, {0, 0, 1}} * t;
-
-  environment::ShowerAxis const showerAxis{injectionPos, (showerCore - injectionPos),
-                                           *env, 20};
-
-  CHECK(showerAxis.steplength() == 500_m);
-
-  CHECK(showerAxis.maximumX() / (10_km * density) == Approx(1).epsilon(1e-8));
-
-  CHECK(showerAxis.minimumX() == 0_g / square(1_cm));
-
-  const Point p{cs, 10_km, 20_km, 8.3_km};
-  CHECK(showerAxis.projectedX(p) / (1.7_km * density) == Approx(1).epsilon(1e-8));
-
-  const units::si::LengthType d = 6.789_km;
-  CHECK(showerAxis.X(d) / (d * density) == Approx(1).epsilon(1e-8));
-
-  const Vector<dimensionless_d> dir{cs, {0, 0, -1}};
-  CHECK(showerAxis.GetDirection().GetComponents(cs) == dir.GetComponents(cs));
-  CHECK(showerAxis.GetStart().GetCoordinates() == injectionPos.GetCoordinates());
-}

Framework/Analytics/CMakeLists.txt

-# create the library
-add_library (CORSIKAanalytics INTERFACE)
-
-# namespace of library -> location of header files
-set (
-  CORSIKAanalytics_NAMESPACE
-  corsika/analytics
-  )
-
-# header files of this library
-set (
-  CORSIKAanalytics_HEADERS
-  ClassTimer.h
-  FunctionTimer.h
-  )
-
-# copy the headers into the namespace
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAanalytics
-  ${CORSIKAanalytics_NAMESPACE}
-  ${CORSIKAanalytics_HEADERS}
-  )
-
-# include directive for upstream code
-target_include_directories (
-  CORSIKAanalytics
-  INTERFACE
-  $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
-  $<INSTALL_INTERFACE:include/>
-  C8::ext:boost
-  )
-
-# and link against spdlog
-#target_link_libraries(
-#  CORSIKAanalytics
-#  INTERFACE
-#  spdlog::spdlog
-#)
-
-# install library
-install (
-  FILES ${CORSIKAanalytics_HEADERS}
-  DESTINATION include/${CORSIKAanalytics_NAMESPACE}
-  )
-
-# ----------------
-# code unit testing
-
-
-CORSIKA_ADD_TEST (testFunctionTimer)
-target_link_libraries (
- testFunctionTimer
- CORSIKAanalytics
- CORSIKAtesting
- )
-
- CORSIKA_ADD_TEST (testClassTimer)
-target_link_libraries (
- testClassTimer
- CORSIKAanalytics
- CORSIKAtesting
- )

Framework/Analytics/FunctionTimer.h

-/*
- * (c) Copyright 2020 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 <chrono>
-#include <utility>
-
-namespace corsika::analytics {
-
-  /// Wraps and measures the runtime of a single function type object
-  /**
-   *
-   * @tparam TFunc funtion pointer that should be wrapped
-   * @tparam TClock type of the clock that should be used for measurements
-   * @tparam TDuration type of std::duration to measure the elapsed time
-   */
-  template <typename TFunc, typename TClock = std::chrono::high_resolution_clock,
-            typename TDuration = std::chrono::microseconds>
-  class FunctionTimer {
-  private:
-    typename TClock::time_point start_;
-    TDuration timeDiff_;
-
-    TFunc function_;
-
-  public:
-    /// Constructs the wrapper with the given functionpointer
-    FunctionTimer(TFunc f)
-        : function_(f) {}
-
-    template <typename... TArgs>
-    auto operator()(TArgs&&... args) -> std::invoke_result_t<TFunc, TArgs...> {
-      start_ = TClock::now();
-      auto tmp = function_(std::forward<TArgs>(args)...);
-      timeDiff_ = std::chrono::duration_cast<TDuration>(TClock::now() - start_);
-      return tmp;
-    }
-
-    inline TDuration getTime() const { return timeDiff_; }
-  };
-
-} // namespace corsika::analytics

Framework/CMakeLists.txt

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

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
-  testCascade.h
-  )
-
-add_library (CORSIKAcascade INTERFACE)
-
-CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAcascade ${CORSIKAcascade_NAMESPACE} ${CORSIKAcascade_HEADERS})
-
-target_link_libraries(
-  CORSIKAcascade
-  INTERFACE
-  CORSIKAsetup
-  CORSIKArandom
-  CORSIKAstackinterface
-  CORSIKAparticles
-  CORSIKAgeometry
-  CORSIKAenvironment
-  CORSIKAprocesssequence
-  CORSIKAunits
-  CORSIKAlogging
-  )
-
-# 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
-CORSIKA_ADD_TEST(testCascade)
-target_link_libraries (
-  testCascade
-  CORSIKAcascade
-  ProcessStackInspector
-  ProcessTrackingLine
-  ProcessNullModel
-  CORSIKAtesting
-  )

Framework/Cascade/Cascade.dox

-/**
-
-
-  Here are have to explain the corsika::cascade::Cascade class and its
-  functionality.
-  
-
-*/
\ No newline at end of file

Framework/Cascade/Cascade.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/environment/Environment.h>
-#include <corsika/logging/Logging.h>
-#include <corsika/process/ProcessReturn.h>
-#include <corsika/random/ExponentialDistribution.h>
-#include <corsika/random/RNGManager.h>
-#include <corsika/random/UniformRealDistribution.h>
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/stack/history/EventType.hpp>
-#include <corsika/stack/history/HistorySecondaryProducer.hpp>
-
-#include <corsika/setup/SetupTrajectory.h>
-
-/*  see Issue 161, we need to include SetupStack only because we need
-    to globally define StackView. This is clearly not nice and should
-    be changed, when possible. It might be that StackView needs to be
-    templated in Cascade, but this would be even worse... so we don't
-    do that until it is really needed.
- */
-#include <corsika/setup/SetupStack.h>
-
-#include <cassert>
-#include <cmath>
-#include <limits>
-
-/**
- * The cascade namespace assembles all objects needed to simulate full particles cascades.
- */
-
-namespace corsika::cascade {
-
-  /**
-   * \class Cascade
-   *
-   * The Cascade class is constructed from template arguments making
-   * it very versatile. Via the template arguments physics models are
-   * plugged into the cascade simulation.
-   *
-   * <b>TTracking</b> must be a class according to the
-   * TrackingInterface providing the functions:
-   * <code>auto GetTrack(Particle const& p)</auto>,
-   * with the return type <code>geometry::Trajectory<corsika::geometry::Line>
-   * </code>
-   *
-   * <b>TProcessList</b> must be a ProcessSequence.   *
-   * <b>Stack</b> is the storage object for particle data, i.e. with
-   * Particle class type <code>Stack::ParticleType</code>
-   *
-   *
-   */
-
-  template <typename TTracking, typename TProcessList, typename TStack,
-            /*
-              TStackView is needed as explicit template parameter because
-              of issue 161 and the
-              inability of clang to understand "stack::MakeView" so far.
-             */
-            typename TStackView = corsika::setup::StackView>
-  class Cascade {
-    using Particle = typename TStack::ParticleType;
-    using VolumeTreeNode =
-        std::remove_pointer_t<decltype(((Particle*)nullptr)->GetNode())>;
-    using MediumInterface = typename VolumeTreeNode::IModelProperties;
-
-  private:
-    // Data members
-    corsika::environment::Environment<MediumInterface> const& environment_;
-    TTracking& tracking_;
-    TProcessList& process_sequence_;
-    TStack& stack_;
-    corsika::random::RNG& rng_ =
-        corsika::random::RNGManager::GetInstance().GetRandomStream("cascade");
-    unsigned int count_ = 0;
-
-  private:
-    // we only want fully configured objects
-    Cascade() = delete;
-
-  public:
-    /**
-     * Cascade class cannot be default constructed, but needs a valid
-     * list of physics processes for configuration at construct time.
-     */
-    Cascade(corsika::environment::Environment<MediumInterface> const& env, TTracking& tr,
-            TProcessList& pl, TStack& stack)
-        : environment_(env)
-        , tracking_(tr)
-        , process_sequence_(pl)
-        , stack_(stack)
-        , count_(0) {
-      C8LOG_INFO(c8_ascii_);
-      if constexpr (TStackView::has_event) {
-        C8LOG_INFO(" - With full cascade HISTORY.");
-      }
-    }
-
-    /**
-     * The Run function is the main simulation loop, which processes
-     * particles from the Stack until the Stack is empty.
-     */
-    void Run() {
-      setNodes();
-
-      while (!stack_.IsEmpty()) {
-        while (!stack_.IsEmpty()) {
-          C8LOG_TRACE("Stack: {}", stack_.as_string());
-          count_++;
-          auto pNext = stack_.GetNextParticle();
-          C8LOG_DEBUG(
-              "============== next particle : count={}, pid={}, "
-              ", stack entries={}"
-              ", stack deleted={}",
-              count_, pNext.GetPID(), stack_.getEntries(), stack_.getDeleted());
-          Step(pNext);
-          process_sequence_.DoStack(stack_);
-        }
-        // do cascade equations, which can put new particles on Stack,
-        // thus, the double loop
-        // DoCascadeEquations();
-      }
-    }
-
-    /**
-     * Force an interaction of the top particle of the stack at its current position.
-     * Note that SetNodes() or an equivalent procedure needs to be called first if you
-     * want to call forceInteraction() for the primary interaction.
-     */
-    void forceInteraction() {
-      C8LOG_DEBUG("forced interaction!");
-      setNodes();
-      auto vParticle = stack_.GetNextParticle();
-      TStackView secondaries(vParticle);
-      interaction(secondaries);
-      process_sequence_.DoSecondaries(secondaries);
-      vParticle.Delete(); // primary particle has interacted and is gone
-    }
-
-  private:
-    /**
-     * The Step function is executed for each particle from the
-     * stack. It will calcualte geometric transport of the particles,
-     * and apply continuous and stochastic processes to it, which may
-     * lead to energy losses, scattering, absorption, decays and the
-     * production of secondary particles.
-     *
-     * New particles produced in one step are subject to further
-     * processing, e.g. thinning, etc.
-     */
-    void Step(Particle& vParticle) {
-      using namespace corsika;
-      using namespace corsika::units::si;
-
-      // determine geometric tracking
-      auto [step, geomMaxLength, nextVol] = tracking_.GetTrack(vParticle);
-      [[maybe_unused]] auto const& dummy_nextVol = nextVol;
-
-      // determine combined total interaction length (inverse)
-      InverseGrammageType const total_inv_lambda =
-          process_sequence_.GetInverseInteractionLength(vParticle);
-
-      // sample random exponential step length in grammage
-      corsika::random::ExponentialDistribution expDist(1 / total_inv_lambda);
-      GrammageType const next_interact = expDist(rng_);
-
-      C8LOG_DEBUG(
-          "total_lambda={} g/cm2, "
-          ", next_interact={} g/cm2",
-          double((1. / total_inv_lambda) / 1_g * 1_cm * 1_cm),
-          double(next_interact / 1_g * 1_cm * 1_cm));
-
-      auto const* currentLogicalNode = vParticle.GetNode();
-
-      // assert that particle stays outside void Universe if it has no
-      // model properties set
-      assert(currentLogicalNode != &*environment_.GetUniverse() ||
-             environment_.GetUniverse()->HasModelProperties());
-
-      // convert next_step from grammage to length
-      LengthType const distance_interact =
-          currentLogicalNode->GetModelProperties().ArclengthFromGrammage(step,
-                                                                         next_interact);
-
-      // determine the maximum geometric step length from continuous processes
-      LengthType const distance_max = process_sequence_.MaxStepLength(vParticle, step);
-      C8LOG_DEBUG("distance_max={} m", distance_max / 1_m);
-
-      // determine combined total inverse decay time
-      InverseTimeType const total_inv_lifetime =
-          process_sequence_.GetInverseLifetime(vParticle);
-
-      // sample random exponential decay time
-      corsika::random::ExponentialDistribution expDistDecay(1 / total_inv_lifetime);
-      TimeType const next_decay = expDistDecay(rng_);
-      C8LOG_DEBUG(
-          "total_lifetime={} s"
-          ", next_decay={} s",
-          (1 / total_inv_lifetime) / 1_s, next_decay / 1_s);
-
-      // convert next_decay from time to length [m]
-      LengthType const distance_decay = next_decay * vParticle.GetMomentum().norm() /
-                                        vParticle.GetEnergy() * units::constants::c;
-
-      // take minimum of geometry, interaction, decay for next step
-      auto const min_distance =
-          std::min({distance_interact, distance_decay, distance_max, geomMaxLength});
-
-      C8LOG_DEBUG("transport particle by : {} m", min_distance / 1_m);
-
-      // here the particle is actually moved along the trajectory to new position:
-      // std::visit(setup::ParticleUpdate<Particle>{vParticle}, step);
-      vParticle.SetPosition(step.PositionFromArclength(min_distance));
-      // .... also update time, momentum, direction, ...
-      vParticle.SetTime(vParticle.GetTime() + min_distance / units::constants::c);
-
-      step.LimitEndTo(min_distance);
-
-      // apply all continuous processes on particle + track
-      if (process_sequence_.DoContinuous(vParticle, step) ==
-          process::EProcessReturn::eParticleAbsorbed) {
-        C8LOG_DEBUG("Cascade: delete absorbed particle PID={} E={} GeV",
-                    vParticle.GetPID(), vParticle.GetEnergy() / 1_GeV);
-        if (!vParticle.isDeleted()) vParticle.Delete();
-        return;
-      }
-
-      C8LOG_DEBUG("sth. happening before geometric limit ? {}",
-                  ((min_distance < geomMaxLength) ? "yes" : "no"));
-
-      if (min_distance < geomMaxLength) { // interaction to happen within geometric limit
-
-        // check whether decay or interaction limits this step the
-        // outcome of decay or interaction MAY be a) new particles in
-        // secondaries, b) the projectile particle deleted (or
-        // changed)
-
-        TStackView secondaries(vParticle);
-
-        if (min_distance != distance_max) {
-          /*
-            Create SecondaryView object on Stack. The data container
-            remains untouched and identical, and 'projectil' is identical
-            to 'vParticle' above this line. However,
-            projectil.AddSecondaries populate the SecondaryView, which can
-            then be used afterwards for further processing. Thus: it is
-            important to use projectle/view (and not vParticle) for Interaction,
-            and Decay!
-          */
-
-          [[maybe_unused]] auto projectile = secondaries.GetProjectile();
-
-          if (min_distance == distance_interact) {
-            interaction(secondaries);
-          } else {
-            assert(min_distance == distance_decay);
-            decay(secondaries);
-            // make sure particle actually did decay if it should have done so
-            if (secondaries.getSize() == 1 &&
-                projectile.GetPID() == secondaries.GetNextParticle().GetPID())
-              throw std::runtime_error(
-                  fmt::format("Cascade: {} decayed into itself!",
-                              particles::GetName(projectile.GetPID())));
-          }
-
-          process_sequence_.DoSecondaries(secondaries);
-          vParticle.Delete();
-
-        } else { // step-length limitation within volume
-          C8LOG_DEBUG("step-length limitation");
-          // no extra physics happens here. just proceed to next step.
-        }
-
-        [[maybe_unused]] auto const assertion = [&] {
-          auto const* numericalNodeAfterStep =
-              environment_.GetUniverse()->GetContainingNode(vParticle.GetPosition());
-          C8LOG_TRACE("Geometry check: numericalNodeAfterStep={} currentLogicalNode={}",
-                      fmt::ptr(numericalNodeAfterStep), fmt::ptr(currentLogicalNode));
-          return numericalNodeAfterStep == currentLogicalNode;
-        };
-
-        assert(assertion()); // numerical and logical nodes don't match
-      } else {               // boundary crossing, step is limited by volume boundary
-        vParticle.SetNode(nextVol);
-        /*
-          DoBoundary may delete the particle (or not)
-
-          caveat: any changes to vParticle, or even the production
-          of new secondaries is currently not passed to ParticleCut,
-          thus, particles outside the desired phase space may be produced.
-
-          todo: this must be fixed.
-        */
-        process_sequence_.DoBoundaryCrossing(vParticle, *currentLogicalNode, *nextVol);
-      }
-    }
-
-    process::EProcessReturn decay(TStackView& view) {
-      C8LOG_DEBUG("decay");
-      units::si::InverseTimeType const actual_decay_time =
-          process_sequence_.GetInverseLifetime(view.parent());
-
-      random::UniformRealDistribution<units::si::InverseTimeType> uniDist(
-          actual_decay_time);
-      const auto sample_process = uniDist(rng_);
-      auto const returnCode = process_sequence_.SelectDecay(view, sample_process);
-      if (returnCode != process::EProcessReturn::eDecayed) {
-        C8LOG_WARN("Particle did not decay!");
-      }
-      SetEventType(view, history::EventType::Decay);
-      return returnCode;
-    }
-
-    process::EProcessReturn interaction(TStackView& view) {
-      C8LOG_DEBUG("collide");
-
-      units::si::InverseGrammageType const current_inv_length =
-          process_sequence_.GetInverseInteractionLength(view.parent());
-
-      random::UniformRealDistribution<units::si::InverseGrammageType> uniDist(
-          current_inv_length);
-      const auto sample_process = uniDist(rng_);
-      auto const returnCode = process_sequence_.SelectInteraction(view, sample_process);
-      if (returnCode != process::EProcessReturn::eInteracted) {
-        C8LOG_WARN("Particle did not interace!");
-      }
-      SetEventType(view, history::EventType::Interaction);
-      return returnCode;
-    }
-
-    /**
-     * set the nodes for all particles on the stack according to their numerical
-     * position
-     */
-    void setNodes() {
-      std::for_each(stack_.begin(), stack_.end(), [&](auto& p) {
-        auto const* numericalNode =
-            environment_.GetUniverse()->GetContainingNode(p.GetPosition());
-        p.SetNode(numericalNode);
-      });
-    }
-
-    void SetEventType(TStackView& view, [[maybe_unused]] history::EventType eventType) {
-      if constexpr (TStackView::has_event) {
-        for (auto&& sec : view) { sec.GetEvent()->setEventType(eventType); }
-      }
-    }
-
-    // but this here temporarily. Should go into dedicated file later:
-    const char* c8_ascii_ =
-        R"V0G0N(
-  ,ad8888ba,     ,ad8888ba,    88888888ba    ad88888ba   88  88      a8P          db              ad88888ba   
- d8"'    `"8b   d8"'    `"8b   88      "8b  d8"     "8b  88  88    ,88'          d88b            d8"     "8b  
-d8'            d8'        `8b  88      ,8P  Y8,          88  88  ,88"           d8'`8b           Y8a     a8P  
-88             88          88  88aaaaaa8P'  `Y8aaaaa,    88  88,d88'           d8'  `8b           "Y8aaa8P"   
-88             88          88  88""""88'      `"""""8b,  88  8888"88,         d8YaaaaY8b          ,d8"""8b,   
-Y8,            Y8,        ,8P  88    `8b            `8b  88  88P   Y8b       d8""""""""8b        d8"     "8b  
- Y8a.    .a8P   Y8a.    .a8P   88     `8b   Y8a     a8P  88  88     "88,    d8'        `8b       Y8a     a8P  
-  `"Y8888Y"'     `"Y8888Y"'    88      `8b   "Y88888P"   88  88       Y8b  d8'          `8b       "Y88888P"
-	)V0G0N";
-  };
-
-} // namespace corsika::cascade

Framework/Cascade/testCascade.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/cascade/testCascade.h>
-
-#include <corsika/cascade/Cascade.h>
-
-#include <corsika/process/ProcessSequence.h>
-#include <corsika/process/null_model/NullModel.h>
-#include <corsika/process/stack_inspector/StackInspector.h>
-#include <corsika/process/tracking_line/TrackingLine.h>
-
-#include <corsika/particles/ParticleProperties.h>
-
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/RootCoordinateSystem.h>
-#include <corsika/geometry/Vector.h>
-
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/NuclearComposition.h>
-
-#include <catch2/catch.hpp>
-
-using namespace corsika;
-using namespace corsika::process;
-using namespace corsika::units;
-using namespace corsika::units::si;
-using namespace corsika::geometry;
-
-#include <limits>
-using namespace std;
-
-auto MakeDummyEnv() {
-  TestEnvironmentType env; // dummy environment
-  auto& universe = *(env.GetUniverse());
-
-  auto theMedium = TestEnvironmentType::CreateNode<Sphere>(
-      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
-      100_km * std::numeric_limits<double>::infinity());
-
-  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
-  theMedium->SetModelProperties<MyHomogeneousModel>(
-      1_g / (1_cm * 1_cm * 1_cm),
-      environment::NuclearComposition(
-          std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.}));
-
-  universe.AddChild(std::move(theMedium));
-
-  return env;
-}
-
-class ProcessSplit : public process::InteractionProcess<ProcessSplit> {
-
-  int fCalls = 0;
-  GrammageType fX0;
-
-public:
-  ProcessSplit(GrammageType const X0)
-      : fX0(X0) {}
-
-  template <typename Particle>
-  corsika::units::si::GrammageType GetInteractionLength(Particle const&) const {
-    return fX0;
-  }
-
-  template <typename TSecondaryView>
-  corsika::process::EProcessReturn DoInteraction(TSecondaryView& view) {
-    fCalls++;
-    auto const projectile = view.GetProjectile();
-    const HEPEnergyType E = projectile.GetEnergy();
-    view.AddSecondary(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            projectile.GetPID(), E / 2, projectile.GetMomentum(),
-            projectile.GetPosition(), projectile.GetTime()});
-    view.AddSecondary(
-        std::tuple<particles::Code, units::si::HEPEnergyType,
-                   corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-            projectile.GetPID(), E / 2, projectile.GetMomentum(),
-            projectile.GetPosition(), projectile.GetTime()});
-    return EProcessReturn::eInteracted;
-  }
-
-  int GetCalls() const { return fCalls; }
-};
-
-class ProcessCut : public process::SecondariesProcess<ProcessCut> {
-
-  int fCount = 0;
-  int fCalls = 0;
-  HEPEnergyType fEcrit;
-
-public:
-  ProcessCut(HEPEnergyType e)
-      : fEcrit(e) {}
-
-  template <typename TStack>
-  EProcessReturn DoSecondaries(TStack& vS) {
-    fCalls++;
-    auto p = vS.begin();
-    while (p != vS.end()) {
-      HEPEnergyType E = p.GetEnergy();
-      if (E < fEcrit) {
-        p.Delete();
-        fCount++;
-      }
-      ++p; // next particle
-    }
-    C8LOG_INFO(fmt::format("ProcessCut::DoSecondaries size={} count={}", vS.getEntries(),
-                           fCount));
-    return EProcessReturn::eOk;
-  }
-
-  int GetCount() const { return fCount; }
-  int GetCalls() const { return fCalls; }
-};
-
-TEST_CASE("Cascade", "[Cascade]") {
-
-  HEPEnergyType E0 = 100_GeV;
-
-  random::RNGManager& rmng = random::RNGManager::GetInstance();
-  rmng.RegisterRandomStream("cascade");
-
-  auto env = MakeDummyEnv();
-  auto const& rootCS = env.GetCoordinateSystem();
-  tracking_line::TrackingLine tracking;
-
-  stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0);
-  null_model::NullModel nullModel;
-
-  const GrammageType X0 = 20_g / square(1_cm);
-  const HEPEnergyType Ecrit = 85_MeV;
-  ProcessSplit split(X0);
-  ProcessCut cut(Ecrit);
-  auto sequence = nullModel % stackInspect % split % cut;
-  TestCascadeStack stack;
-  stack.Clear();
-  stack.AddParticle(
-      std::tuple<particles::Code, units::si::HEPEnergyType,
-                 corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
-          particles::Code::Electron, E0,
-          corsika::stack::MomentumVector(rootCS, {0_GeV, 0_GeV, -1_GeV}),
-          Point(rootCS, {0_m, 0_m, 10_km}), 0_ns});
-
-  cascade::Cascade<tracking_line::TrackingLine, decltype(sequence), TestCascadeStack,
-                   TestCascadeStackView>
-      EAS(env, tracking, sequence, stack);
-
-  SECTION("full cascade") {
-    EAS.Run();
-
-    CHECK(cut.GetCount() == 2048);
-    CHECK(cut.GetCalls() == 2047);
-    CHECK(split.GetCalls() == 2047);
-  }
-
-  SECTION("forced interaction") {
-    EAS.forceInteraction();
-    CHECK(stack.getEntries() == 2);
-    CHECK(split.GetCalls() == 1);
-  }
-}

Framework/Cascade/testCascade.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/environment/Environment.h>
-
-#include <corsika/stack/CombinedStack.h>
-#include <corsika/stack/SecondaryView.h>
-#include <corsika/stack/node/GeometryNodeStackExtension.h>
-#include <corsika/stack/nuclear_extension/NuclearStackExtension.h>
-
-using TestEnvironmentType =
-    corsika::environment::Environment<corsika::environment::IMediumModel>;
-
-template <typename T>
-using SetupGeometryDataInterface =
-    corsika::stack::node::GeometryDataInterface<T, TestEnvironmentType>;
-
-// combine particle data stack with geometry information for tracking
-template <typename StackIter>
-using StackWithGeometryInterface = corsika::stack::CombinedParticleInterface<
-    corsika::stack::nuclear_extension::ParticleDataStack::MPIType,
-    SetupGeometryDataInterface, StackIter>;
-
-using TestCascadeStack = corsika::stack::CombinedStack<
-    typename corsika::stack::nuclear_extension::ParticleDataStack::StackImpl,
-    corsika::stack::node::GeometryData<TestEnvironmentType>, StackWithGeometryInterface>;
-
-/*
-  See also Issue 161
-*/
-#if defined(__clang__)
-using TestCascadeStackView =
-    corsika::stack::SecondaryView<typename TestCascadeStack::StackImpl,
-                                  StackWithGeometryInterface>;
-#elif defined(__GNUC__) || defined(__GNUG__)
-using TestCascadeStackView = corsika::stack::MakeView<TestCascadeStack>::type;
-#endif

Framework/Geometry/BaseTrajectory.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>
-#include <string>
-
-namespace corsika::geometry {
-
-  /*!
-   * Interface / base class for trajectories.
-   */
-  class BaseTrajectory {
-
-    BaseTrajectory() = delete;
-
-  public:
-    BaseTrajectory(corsika::units::si::TimeType start, corsika::units::si::TimeType end)
-        : fTStart(start)
-        , fTEnd(end) {}
-
-    //!< for \f$ t = 0 \f$, the starting Point shall be returned.
-    virtual Point GetPosition(corsika::units::si::TimeType) const = 0;
-
-    //!< the Point is return from u=0 (start) to u=1 (end)
-    virtual Point GetPosition(double u) const = 0;
-
-    /*!
-     * returns the length between two points of the trajectory
-     * parameterized by \arg t1 and \arg t2. Requires \arg t2 > \arg t1.
-     */
-
-    virtual corsika::units::si::TimeType TimeFromArclength(
-        corsika::units::si::LengthType) const = 0;
-
-    virtual LengthType ArcLength(corsika::units::si::TimeType t1,
-                                 corsika::units::si::TimeType t2) const = 0;
-
-    virtual corsika::units::si::TimeType GetDuration(
-        corsika::units::si::TimeType t1, corsika::units::si::TimeType t2) const {
-      return t2 - t1;
-    }
-
-    virtual Point GetEndpoint() const { return GetPosition(fTEnd); }
-    virtual Point GetStartpoint() const { return GetPosition(fTStart); }
-
-  protected:
-    corsika::units::si::TimeType const fTStart, fTEnd;
-  };
-
-} // namespace corsika::geometry

Framework/Geometry/BaseVector.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/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

Framework/Geometry/CMakeLists.txt

-set (
-  GEOMETRY_SOURCES
-  CoordinateSystem.cc
-  )
-
-set (
-  GEOMETRY_HEADERS
-  Vector.h
-  Point.h
-  Line.h
-  Sphere.h
-  Plane.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
-  CORSIKAutilities
-  C8::ext::eigen3
-  )
-
-target_include_directories (
-  CORSIKAgeometry
-  INTERFACE
-  $<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
-CORSIKA_ADD_TEST(testGeometry)
-target_link_libraries (
-  testGeometry
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAtesting
-  )
-
-CORSIKA_ADD_TEST(testFourVector)
-target_link_libraries (
-  testFourVector
-  CORSIKAgeometry
-  CORSIKAunits
-  CORSIKAtesting
-  )

Framework/Geometry/CoordinateSystem.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/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
- *
- * 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/QuantityVector.h>
-#include <corsika/units/PhysicalUnits.h>
-#include <corsika/utl/sgn.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;
-  template <typename T>
-  class Vector;
-
-  using corsika::units::si::length_d;
-
-  class CoordinateSystem {
-    CoordinateSystem const* reference = nullptr;
-    EigenTransform 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);
-
-    CoordinateSystem(CoordinateSystem const& reference, EigenTransform const& transf)
-        : reference(&reference)
-        , transf(transf) {}
-
-    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);
-    }
-
-    /**
-     * creates a new CS in which vVec points in direction of the new z-axis
-     */
-    template <typename TDim>
-    auto RotateToZ(Vector<TDim> vVec) const {
-      auto const a = vVec.normalized().GetComponents(*this).eVector;
-      auto const a1 = a(0), a2 = a(1);
-
-      auto const s = utl::sgn(a(2));
-      auto const c = 1 / (1 + s * a(2));
-
-      Eigen::Matrix3d A, B;
-
-      if (s > 0) {
-        A << 1, 0, a1,                      // comment to prevent clang-format
-            0, 1, a2,                       // .
-            -a1, -a2, 1;                    // .
-        B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
-            -a1 * a2 * c, -a2 * a2 * c, 0,  // .
-            0, 0, -(a1 * a1 + a2 * a2) * c; // .
-
-      } else {
-        A << 1, 0, a1,                      // .
-            0, -1, a2,                      // .
-            a1, -a2, -1;                    // .
-        B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
-            -a1 * a2 * c, +a2 * a2 * c, 0,  // .
-            0, 0, (a1 * a1 + a2 * a2) * c;  // .
-      }
-
-      return CoordinateSystem(*this, EigenTransform(A + B));
-    }
-
-    template <typename TDim>
-    auto rotate(QuantityVector<TDim> 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);
-    }
-
-    template <typename TDim>
-    auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
-                            QuantityVector<TDim> 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; }
-
-    bool operator==(CoordinateSystem const& cs) const {
-      return reference == cs.reference && transf.matrix() == cs.transf.matrix();
-    }
-
-    bool operator!=(CoordinateSystem const& cs) const { return !(cs == *this); }
-  };
-
-} // namespace corsika::geometry

Framework/Geometry/FourVector.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/Vector.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <iostream>
-
-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

Framework/Geometry/Helix.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>
-#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