Environment/HomogeneousMedium.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_HomogeneousMedium_h_
-#define _include_HomogeneousMedium_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>
-
-#include <cassert>
-
-/**
- * a homogeneous medium
- */
-
-namespace corsika::environment {
-
-  template <class T>
-  class HomogeneousMedium : public T {
-    corsika::units::si::MassDensityType const fDensity;
-    NuclearComposition const fNuclComp;
-
-  public:
-    HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
-                      NuclearComposition pNuclComp)
-        : fDensity(pDensity)
-        , fNuclComp(pNuclComp) {}
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const override {
-      return fDensity;
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::LengthType pTo) const override {
-      using namespace corsika::units::si;
-      return pTo * fDensity;
-    }
-
-    corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return pGrammage / fDensity;
-    }
-  };
-
-} // namespace corsika::environment
-#endif

Environment/IMediumModel.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_IMediumModel_h
-#define _include_IMediumModel_h
-
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-#include <corsika/units/PhysicalUnits.h>
-
-namespace corsika::environment {
-
-  class IMediumModel {
-  public:
-    virtual ~IMediumModel() = default;
-
-    virtual corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const&) const = 0;
-
-    // todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
-    // approach for now, only lines are supported
-    virtual corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::LengthType) const = 0;
-
-    virtual corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const&,
-        corsika::units::si::GrammageType) const = 0;
-
-    virtual NuclearComposition const& GetNuclearComposition() const = 0;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/InhomogeneousMedium.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_environment_InhomogeneousMedium_h_
-#define _include_environment_InhomogeneousMedium_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>
-
-/**
- * A general inhomogeneous medium. The mass density distribution TDensityFunction must be
- * a \f$C^2\f$-function.
- */
-
-namespace corsika::environment {
-
-  template <class T, class TDensityFunction>
-  class InhomogeneousMedium : public T {
-    NuclearComposition const fNuclComp;
-    TDensityFunction const fDensityFunction;
-
-  public:
-    template <typename... Args>
-    InhomogeneousMedium(NuclearComposition pNuclComp, Args&&... rhoArgs)
-        : fNuclComp(pNuclComp)
-        , fDensityFunction(rhoArgs...) {}
-
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const& p) const override {
-      return fDensityFunction.EvaluateAt(p);
-    }
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    corsika::units::si::GrammageType IntegratedGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
-        corsika::units::si::LengthType pTo) const override {
-      return fDensityFunction.IntegrateGrammage(pLine, pTo);
-    }
-
-    corsika::units::si::LengthType ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
-        corsika::units::si::GrammageType pGrammage) const override {
-      return fDensityFunction.ArclengthFromGrammage(pLine, pGrammage);
-    }
-  };
-
-} // namespace corsika::environment
-#endif

Environment/LinearApproximationIntegrator.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_environment_LinearApproximationIntegrator_h_
-#define _include_environment_LinearApproximationIntegrator_h_
-
-#include <limits>
-
-#include <corsika/geometry/Line.h>
-#include <corsika/geometry/Point.h>
-#include <corsika/geometry/Trajectory.h>
-
-namespace corsika::environment {
-  template <class TDerived>
-  class LinearApproximationIntegrator {
-    auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
-
-  public:
-    auto IntegrateGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-        corsika::units::si::LengthType length) const {
-      auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-      return (c0 + 0.5 * c1 * length) * length;
-    }
-
-    auto ArclengthFromGrammage(
-        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-        corsika::units::si::GrammageType grammage) const {
-      auto const c0 = GetImplementation().fRho(line.GetPosition(0));
-      auto const c1 = GetImplementation().fRho.FirstDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
-    }
-
-    auto MaximumLength(corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-                       [[maybe_unused]] double relError) const {
-      using namespace corsika::units::si;
-      [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
-          line.GetPosition(0), line.NormalizedDirection());
-
-      // todo: provide a real, working implementation
-      return 1_m * std::numeric_limits<double>::infinity();
-    }
-  };
-} // namespace corsika::environment
-
-#endif

Environment/NameModel.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_NameModel_h
-#define _include_NameModel_h
-
-#include <string>
-#include <utility>
-
-namespace corsika::environment {
-
-  template <typename T>
-  struct NameModel : public T {
-    virtual std::string const& GetName() const = 0;
-    virtual ~NameModel() = default;
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/NuclearComposition.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_NuclearComposition_h
-#define _include_NuclearComposition_h
-
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-#include <cassert>
-#include <numeric>
-#include <random>
-#include <stdexcept>
-#include <vector>
-
-namespace corsika::environment {
-  class NuclearComposition {
-    std::vector<float> const fNumberFractions; //!< relative fractions of number density
-    std::vector<corsika::particles::Code> const
-        fComponents; //!< particle codes of consitutents
-
-    double const fAvgMassNumber;
-
-    template <class AConstIterator, class BConstIterator>
-    class WeightProviderIterator {
-      AConstIterator fAIter;
-      BConstIterator fBIter;
-
-    public:
-      using value_type = double;
-      using iterator_category = std::input_iterator_tag;
-      using pointer = value_type*;
-      using reference = value_type&;
-      using difference_type = ptrdiff_t;
-
-      WeightProviderIterator(AConstIterator a, BConstIterator b)
-          : fAIter(a)
-          , fBIter(b) {}
-
-      value_type operator*() const { return ((*fAIter) * (*fBIter)).magnitude(); }
-
-      WeightProviderIterator& operator++() { // prefix ++
-        ++fAIter;
-        ++fBIter;
-        return *this;
-      }
-
-      auto operator==(WeightProviderIterator other) { return fAIter == other.fAIter; }
-
-      auto operator!=(WeightProviderIterator other) { return !(*this == other); }
-    };
-
-  public:
-    NuclearComposition(std::vector<corsika::particles::Code> pComponents,
-                       std::vector<float> pFractions)
-        : fNumberFractions(pFractions)
-        , fComponents(pComponents)
-        , fAvgMassNumber(std::inner_product(
-              pComponents.cbegin(), pComponents.cend(), pFractions.cbegin(), 0.,
-              std::plus<double>(), [](auto const compID, auto const fraction) -> double {
-                if (particles::IsNucleus(compID)) {
-                  return particles::GetNucleusA(compID) * fraction;
-                } else {
-                  return particles::GetMass(compID) /
-                         units::si::ConvertSIToHEP(units::constants::u) * fraction;
-                }
-              })) {
-      assert(pComponents.size() == pFractions.size());
-      auto const sumFractions =
-          std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
-
-      if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
-        throw std::runtime_error("element fractions do not add up to 1");
-      }
-    }
-
-    template <typename TFunction>
-    auto WeightedSum(TFunction func) const {
-      using ResultQuantity = decltype(func(*fComponents.cbegin()));
-
-      auto const sum = [](auto x, auto y) { return x + y; };
-      auto const prod = [&](auto const compID, auto const fraction) {
-        return func(compID) * fraction;
-      };
-
-      if constexpr (phys::units::is_quantity_v<ResultQuantity>) {
-        return std::inner_product(
-            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
-            ResultQuantity::zero(), // .zero() is defined for quantity types only
-            sum, prod);
-      } else {
-        return std::inner_product(
-            fComponents.cbegin(), fComponents.cend(), fNumberFractions.cbegin(),
-            ResultQuantity(0), // in other cases we have to use a bare 0
-            sum, prod);
-      }
-    }
-
-    auto size() const { return fNumberFractions.size(); }
-
-    auto const& GetFractions() const { return fNumberFractions; }
-    auto const& GetComponents() const { return fComponents; }
-    auto const GetAverageMassNumber() const { return fAvgMassNumber; }
-
-    template <class TRNG>
-    corsika::particles::Code SampleTarget(
-        std::vector<corsika::units::si::CrossSectionType> const& sigma,
-        TRNG& randomStream) const {
-      using namespace corsika::units::si;
-
-      assert(sigma.size() == fNumberFractions.size());
-
-      std::discrete_distribution channelDist(
-          WeightProviderIterator<decltype(fNumberFractions.begin()),
-                                 decltype(sigma.begin())>(fNumberFractions.begin(),
-                                                          sigma.begin()),
-          WeightProviderIterator<decltype(fNumberFractions.begin()),
-                                 decltype(sigma.end())>(fNumberFractions.end(),
-                                                        sigma.end()));
-
-      auto const iChannel = channelDist(randomStream);
-      return fComponents[iChannel];
-    }
-  };
-
-} // namespace corsika::environment
-
-#endif

Environment/SlidingPlanarExponential.h

-
-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_Environment_SlidingPlanarExponential_h_
-#define _include_Environment_SlidingPlanarExponential_h_
-
-#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 fNuclComp;
-
-    using Base = BaseExponential<SlidingPlanarExponential<T>>;
-
-  public:
-    SlidingPlanarExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
-                             units::si::LengthType vLambda, NuclearComposition vNuclComp)
-        : Base(vP0, vRho, vLambda)
-        , fNuclComp(vNuclComp) {}
-
-    units::si::MassDensityType GetMassDensity(
-        geometry::Point const& vP) const override {
-      auto const height = (vP - Base::fP0).norm();
-      return Base::fRho0 * exp(Base::fInvLambda * height);
-    }
-
-    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
-
-    units::si::GrammageType IntegratedGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::LengthType vL) const override {
-      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
-      return Base::IntegratedGrammage(vLine, vL, axis);
-    }
-
-    units::si::LengthType ArclengthFromGrammage(
-        geometry::Trajectory<geometry::Line> const& vLine,
-        units::si::GrammageType vGrammage) const override {
-      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
-      return Base::ArclengthFromGrammage(vLine, vGrammage, axis);
-    }
-  };
-
-} // namespace corsika::environment
-#endif

Environment/VolumeTreeNode.h

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#ifndef _include_VolumeTreeNode_H
-#define _include_VolumeTreeNode_H
-
-#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
-
-#endif

Environment/testEnvironment.cc

-/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * See file AUTHORS for a list of contributors.
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
-                          // cpp file
-
-#include <corsika/environment/DensityFunction.h>
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/InhomogeneousMedium.h>
-#include <corsika/environment/LinearApproximationIntegrator.h>
-#include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/SlidingPlanarExponential.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::si::detail::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::si::detail::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::si::detail::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)));
-  }
-}

Framework/CMakeLists.txt

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