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/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/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/Cascade.dox

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