From e6ad6c0ff873c80b6ab26b24e27c035b5811cb2d Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Mon, 25 Jan 2021 16:16:41 +0100
Subject: [PATCH] added and tested Exponential Refractive Index class in media
---
.../media/ExponentialRefractiveIndex.inl | 31 +++++++
corsika/media/ExponentialRefractiveIndex.hpp | 60 ++++++++++++++
tests/media/testRefractiveIndex.cpp | 80 +++++++++++++++++++
3 files changed, 171 insertions(+)
create mode 100644 corsika/detail/media/ExponentialRefractiveIndex.inl
create mode 100644 corsika/media/ExponentialRefractiveIndex.hpp
diff --git a/corsika/detail/media/ExponentialRefractiveIndex.inl b/corsika/detail/media/ExponentialRefractiveIndex.inl
new file mode 100644
index 000000000..664b47b39
--- /dev/null
+++ b/corsika/detail/media/ExponentialRefractiveIndex.inl
@@ -0,0 +1,31 @@
+///*
+// * (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 <bits/stdc++.h>
+//#include <corsika/media/IRefractiveIndexModel.hpp>
+//
+//namespace corsika {
+//
+// template <typename T>
+// template <typename... Args>
+// ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(double const n0,
+// InverseLengthType const lambda, Args&&... args)
+// : T(std::forward<Args>(args)...)
+// , n_0(n0)
+// , lambda_(lambda) {}
+//
+// template <typename T>
+// double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
+// //TODO: THIS METHOD CURRENTLY ONLY USES THE Z-COORDINATE.
+// //NEED TO THINK IT FOR FUTURE WORK ON ARBITRARY GEOMETRIES.
+// return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+// }
+//
+//} // namespace corsika
diff --git a/corsika/media/ExponentialRefractiveIndex.hpp b/corsika/media/ExponentialRefractiveIndex.hpp
new file mode 100644
index 000000000..3ac778926
--- /dev/null
+++ b/corsika/media/ExponentialRefractiveIndex.hpp
@@ -0,0 +1,60 @@
+/*
+ * (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 <bits/stdc++.h>
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+ /**
+ * An exponential refractive index.
+ *
+ * This class returns the value of an exponential refractive index
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class ExponentialRefractiveIndex : public T {
+
+ double n_0; ///< n0 constant.
+ InverseLengthType lambda_; ///< lambda parameter.
+
+ public:
+ /**
+ * Construct an ExponentialRefractiveIndex.
+ *
+ * This is initialized with two parameters n_0 and lambda
+ * and returns the value of the exponential refractive index
+ * at a given point location.
+ *
+ * @param field The refractive index to return to a given point.
+ */
+ template <typename... Args>
+ ExponentialRefractiveIndex(double const n0,
+ InverseLengthType const lambda, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , n_0(n0)
+ , lambda_(lambda) {}
+ /**
+ * Evaluate the refractive index at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this point.
+ */
+ double getRefractiveIndex(Point const& point) const final override
+ {
+ //TODO: THIS METHOD CURRENTLY ONLY USES THE Z-COORDINATE.
+ //NEED TO THINK IT FOR FUTURE WORK ON ARBITRARY GEOMETRIES.
+ return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+ }
+
+ }; // END: class ExponentialRefractiveIndex
+
+} // namespace corsika
\ No newline at end of file
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index fb1e71fec..1888f6a08 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -16,6 +16,7 @@
#include <corsika/media/InhomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
+#include <corsika/media/ExponentialRefractiveIndex.hpp>
#include <SetupTestTrajectory.hpp>
@@ -89,3 +90,82 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
}
+
+
+TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
+
+ logging::set_level(logging::level::info);
+ corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+ // get a CS and a point
+ CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+ Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+ // setup our interface types
+ using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+ using AtmModel = ExponentialRefractiveIndex<HomogeneousMedium<IModelInterface>>;
+
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+
+ // the composition we use for the homogenous medium
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ // a new refractive index
+ const double n0{2};
+ const InverseLengthType lambda{56 / 1_m};
+
+ // create the atmospheric model and check refractive index
+ AtmModel medium(n0, lambda, density, protonComposition);
+ CHECK(n0 == medium.getRefractiveIndex(Point(gCS, 10_m, 14_m, 0_m)));
+
+ // another refractive index
+ const double n0_{1};
+ const InverseLengthType lambda_{1 / 1_km};
+
+ // create the atmospheric model and check refractive index
+ AtmModel medium_(n0_, lambda_, density, protonComposition);
+ CHECK(medium_.getRefractiveIndex(Point(gCS, 12_m, 56_m, 1_km)) == Approx(0.3678794));
+
+ // another refractive index
+ const double n0__{4};
+ const InverseLengthType lambda__{2 / 1_m};
+
+ // create the atmospheric model and check refractive index
+ AtmModel medium__(n0__, lambda__, density, protonComposition);
+ CHECK(medium__.getRefractiveIndex(Point(gCS, 0_m, 0_m, 35_km)) == Approx(0));
+
+ // define our axis vector
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+ // check the density and nuclear composition
+ REQUIRE(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+ medium.getNuclearComposition();
+ REQUIRE(density == medium_.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+ medium_.getNuclearComposition();
+ REQUIRE(density == medium__.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+ medium__.getNuclearComposition();
+
+ // create a line of length 1 m
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {1_m / second, 0_m / second, 0_m / second}));
+
+ // the end time of our line
+ auto const tEnd = 1_s;
+
+ // and the associated trajectory
+ setup::Trajectory const track =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd);
+// // and the associated trajectory
+// Trajectory<Line> const trajectory(line, tEnd);
+
+ // and check the integrated grammage
+ REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+ REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+ REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+}
--
GitLab