Environment/ShowerAxis.cc

-/*
- * (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.
- */
-
-#include <corsika/environment/ShowerAxis.h>
-#include <sstream>
-
-using namespace corsika::environment;
-using namespace corsika::units::si;
-using namespace corsika;
-
-GrammageType ShowerAxis::X(LengthType l) const {
-  auto const fractionalBin = l / steplength_;
-  int const lower = fractionalBin; // indices of nearest X support points
-  auto const lambda = fractionalBin - lower;
-  decltype(X_.size()) const upper = lower + 1;
-
-  if (lower < 0) {
-    throw std::runtime_error("cannot extrapolate to points behind point of injection");
-  }
-
-  if (upper >= X_.size()) {
-    std::stringstream errormsg;
-    errormsg << "shower axis too short, cannot extrapolate (l / max_length_ = "
-             << l / max_length_ << ")";
-    throw std::runtime_error(errormsg.str().c_str());
-  }
-
-  assert(0 <= lambda && lambda <= 1.);
-
-  std::cout << l << ": " << lower << " " << lambda << " " << upper << std::endl;
-
-  // linear interpolation between X[lower] and X[upper]
-  return X_[upper] * lambda + X_[lower] * (1 - lambda);
-}
-
-LengthType ShowerAxis::steplength() const { return steplength_; }
-
-GrammageType ShowerAxis::maximumX() const { return *X_.rbegin(); }
-
-GrammageType ShowerAxis::minimumX() const { return GrammageType::zero(); }
-
-GrammageType ShowerAxis::projectedX(geometry::Point const& p) const {
-  auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
-  return X(projectedLength);
-}
-
-geometry::Vector<units::si::dimensionless_d> const& ShowerAxis::GetDirection() const {
-  return axis_normalized_;
-}
-
-geometry::Point const& ShowerAxis::GetStart() const { return pointStart_; }

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 {
-  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/testShowerAxis.cc

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

-add_subdirectory (Testing)
-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)

Framework/Cascade/Cascade.dox

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