From 7bbfe4c71087aab08fa692202e97880a624459a7 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Wed, 30 Jan 2019 20:16:14 +0100
Subject: [PATCH] inhomogeneous medium + test (still WIP)
---
Environment/CMakeLists.txt | 4 ++
Environment/HomogeneousMedium.h | 2 +-
Environment/InhomogeneousMedium.h | 60 ++++++++++++++++++++++++++
Environment/testEnvironment.cc | 71 ++++++++++++++++++-------------
4 files changed, 107 insertions(+), 30 deletions(-)
create mode 100644 Environment/InhomogeneousMedium.h
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index a476b89e..b04a58bd 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -4,6 +4,10 @@ set (
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.h
+ InhomogeneousMedium.h
+ HomogeneousMedium.h
+ LinearIntegrator.h
+ DensityFunction.h
Environment.h
)
diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
index d058f67c..02f7aabc 100644
--- a/Environment/HomogeneousMedium.h
+++ b/Environment/HomogeneousMedium.h
@@ -36,7 +36,7 @@ namespace corsika::environment {
HomogeneousMedium(corsika::units::si::MassDensityType pDensity,
NuclearComposition pNuclComp)
: fDensity(pDensity)
- , fNuclComp(pNuclComp){};
+ , fNuclComp(pNuclComp) {}
corsika::units::si::MassDensityType GetMassDensity(
corsika::geometry::Point const&) const override {
diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h
new file mode 100644
index 00000000..c506fb17
--- /dev/null
+++ b/Environment/InhomogeneousMedium.h
@@ -0,0 +1,60 @@
+/**
+ * (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>
+
+/**
+ * An inhomogeneous medium. The mass density distribution TDensityFunction must be a
+ * \f$C^1\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.IntegratedGrammage(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
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 0695d02e..01edec05 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -11,14 +11,18 @@
#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/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/InhomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.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>
-#include <random>
-#include <vector>
using namespace corsika::geometry;
using namespace corsika::environment;
@@ -31,32 +35,41 @@ TEST_CASE("HomogeneousMedium") {
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
}
-TEST_CASE("NuclearComposition") {
- NuclearComposition const composition(
- std::vector<corsika::particles::Code>{corsika::particles::Code::Proton,
- corsika::particles::Code::Neutron},
- std::vector<float>{2.f / 3.f, 1.f / 3.f});
- SECTION("SampleTarget") {
- std::vector<CrossSectionType> crossSections{50_mbarn, 100_mbarn};
-
- std::mt19937 rng;
-
- int proton{0}, neutron{0};
-
- for (int i = 0; i < 1'000'000; ++i) {
- corsika::particles::Code p = composition.SampleTarget(crossSections, rng);
- switch (p) {
- case corsika::particles::Code::Proton:
- proton++;
- break;
- case corsika::particles::Code::Neutron:
- neutron++;
- break;
- default:
- throw std::runtime_error("");
- }
- }
-
- REQUIRE(static_cast<double>(proton) / neutron == Approx(1).epsilon(1e-2));
+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);
}
+};
+
+TEST_CASE("DensityFunction") {
+ CoordinateSystem& cs = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ Point const origin(cs, {0_m, 0_m, 0_m});
+
+ Vector direction(cs, QuantityVector<dimensionless_d>(1, 0, 0));
+
+ Line line(origin, Vector<SpeedType::dimension_type>(
+ cs, {200_m / second, 0_m / second, 0_m / second}));
+
+ auto const tEnd = 5_s;
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ Exponential const e;
+ REQUIRE(e.Derivative<1>(origin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+ Approx(1));
+
+ DensityFunction const rho(e);
+ REQUIRE(rho.EvaluateAt(origin) == e(origin));
+
+ auto const exactGrammage = [](auto x) { return rho0 * exp(x / 1_m); };
+ REQUIRE(rho.IntegrateGrammage(trajectory, 5_cm) / exactGrammage(5_cm) ==
+ Approx(1).epsilon(1e-2));
}
--
GitLab