From baa8946af7fe60429a2d5bfd8d34efba5c09e3ec Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Tue, 2 Jun 2020 20:12:18 +0200
Subject: [PATCH] implemented ShowerAxis
---
Environment/CMakeLists.txt | 2 +
Environment/ShowerAxis.cc | 41 +++++++++++++++++
Environment/ShowerAxis.h | 94 ++++++++++++++++++++++++++++++++++++++
3 files changed, 137 insertions(+)
create mode 100644 Environment/ShowerAxis.cc
create mode 100644 Environment/ShowerAxis.h
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index f4b93c759..e01188125 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -1,6 +1,7 @@
set (
ENVIRONMENT_SOURCES
LayeredSphericalAtmosphereBuilder.cc
+ ShowerAxis.cc
)
set (
@@ -19,6 +20,7 @@ set (
FlatExponential.h
SlidingPlanarExponential.h
LayeredSphericalAtmosphereBuilder.h
+ ShowerAxis.h
)
set (
diff --git a/Environment/ShowerAxis.cc b/Environment/ShowerAxis.cc
new file mode 100644
index 000000000..e7dd237d5
--- /dev/null
+++ b/Environment/ShowerAxis.cc
@@ -0,0 +1,41 @@
+/*
+ * (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/ShowerAxis.h>
+
+using namespace corsika::environment;
+using namespace corsika::units::si;
+
+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;
+ int const upper = lower + 1;
+
+ if (upper >= steps) {
+ throw std::runtime_error("shower axis too short, cannot extrapolate");
+ } else if (lower < 0) {
+ throw std::runtime_error("cannot extrapolate to points behind point of injection");
+ }
+
+ assert(0 <= lambda && lambda <= 1.);
+ return (*X_)[lower] * lambda + (*X_)[upper] * (1 - lambda);
+}
+
+LengthType ShowerAxis::steplength() const { return steplength_; }
+
+GrammageType ShowerAxis::maximumX() const { return *(X_->rbegin()); }
+
+GrammageType ShowerAxis::minimumX() const { return *(X_->cbegin()); }
+
+GrammageType ShowerAxis::projectedX(geometry::Point const& p) const {
+ auto const projectedLength = (p - pointStart_).dot(axis_normalized_);
+ return X(projectedLength);
+}
diff --git a/Environment/ShowerAxis.h b/Environment/ShowerAxis.h
new file mode 100644
index 000000000..4b7c27346
--- /dev/null
+++ b/Environment/ShowerAxis.h
@@ -0,0 +1,94 @@
+/*
+ * (c) Copyright 2020 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/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 <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)
+ : pointStart_(pStart)
+ , length_(length)
+ , max_length_(length_.norm())
+ , steplength_(max_length_ / steps)
+ , axis_normalized_(length / max_length_) {
+ 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;
+ for (int i = 1; i < steps; ++i) {
+ auto const x_prev = (i - 1) / (steps - 1.);
+ auto const d_prev = max_length_ * x_prev;
+ auto const x = i / (steps - 1.);
+ 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_;
+ auto const resultTotal = result + (*X_)[i - 1];
+ (*X_)[i] = resultTotal;
+
+ for (; resultTotal > 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;
+
+ private:
+ static constexpr int steps = 10'000;
+ units::si::GrammageType const X_binning_ = std::invoke([]() {
+ using namespace units::si;
+ return 1_g / 1_cm / 1_cm;
+ });
+ std::unique_ptr<std::array<units::si::GrammageType, steps>> const X_{
+ new std::array<units::si::GrammageType, steps>};
+ 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::LengthType> d_;
+ };
+} // namespace corsika::environment
--
GitLab