From 676c72e695213db960647dd1374933b30950b9f8 Mon Sep 17 00:00:00 2001
From: Remy Prechelt <prechelt@hawaii.edu>
Date: Mon, 16 Nov 2020 23:17:52 -1000
Subject: [PATCH] Port over refractive index models.
---
.../detail/media/UniformRefractiveIndex.inl | 31 +++++++
corsika/media/IRefractiveIndexModel.hpp | 40 +++++++++
corsika/media/UniformRefractiveIndex.hpp | 59 +++++++++++++
tests/media/testRefractiveIndex.cpp | 87 +++++++++++++++++++
4 files changed, 217 insertions(+)
create mode 100644 corsika/detail/media/UniformRefractiveIndex.inl
create mode 100644 corsika/media/IRefractiveIndexModel.hpp
create mode 100644 corsika/media/UniformRefractiveIndex.hpp
create mode 100644 tests/media/testRefractiveIndex.cpp
diff --git a/corsika/detail/media/UniformRefractiveIndex.inl b/corsika/detail/media/UniformRefractiveIndex.inl
new file mode 100644
index 000000000..1abea2470
--- /dev/null
+++ b/corsika/detail/media/UniformRefractiveIndex.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 <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+ template <typename T>
+ template <typename... Args>
+ UniformRefractiveIndex<T>::UniformRefractiveIndex(double const n, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , n_(n) {}
+
+ template <typename T>
+ double UniformRefractiveIndex<T>::getRefractiveIndex(Point const&) const {
+ return n_;
+ }
+
+ template <typename T>
+ void UniformRefractiveIndex<T>::setRefractiveIndex(double const& n) {
+ n_ = n;
+ }
+
+} // namespace corsika
diff --git a/corsika/media/IRefractiveIndexModel.hpp b/corsika/media/IRefractiveIndexModel.hpp
new file mode 100644
index 000000000..b042540fb
--- /dev/null
+++ b/corsika/media/IRefractiveIndexModel.hpp
@@ -0,0 +1,40 @@
+/*
+ * (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/framework/geometry/Point.hpp>
+
+namespace corsika {
+
+ /**
+ * An interface for refractive index models.
+ *
+ * This is the base interface for refractive index mixins.
+ *
+ */
+ template <typename TModel>
+ class IRefractiveIndexModel : public TModel {
+
+ public:
+ /**
+ * Evaluate the refractive index at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this point.
+ */
+ virtual double getRefractiveIndex(Point const&) const = 0;
+
+ /**
+ * A virtual default destructor.
+ */
+ virtual ~IRefractiveIndexModel() = default;
+
+ }; // END: class IRefractiveIndexModel
+
+} // namespace corsika
diff --git a/corsika/media/UniformRefractiveIndex.hpp b/corsika/media/UniformRefractiveIndex.hpp
new file mode 100644
index 000000000..2dc5c790a
--- /dev/null
+++ b/corsika/media/UniformRefractiveIndex.hpp
@@ -0,0 +1,59 @@
+/*
+ * (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/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+ /**
+ * A uniform refractive index.
+ *
+ * This class returns the same refractive index
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class UniformRefractiveIndex final : public T {
+
+ double n_; ///< The constant refractive index that we use.
+
+ public:
+ /**
+ * Construct a UniformRefractiveIndex.
+ *
+ * This is initialized with a fixed refractive index
+ * and returns this refractive index at all locations.
+ *
+ * @param field The refractive index to return everywhere.
+ */
+ template <typename... Args>
+ UniformRefractiveIndex(double const n, Args&&... args);
+
+ /**
+ * 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&) const override;
+
+ /**
+ * Set the refractive index returned by this instance.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this location.
+ */
+ void setRefractiveIndex(double const& n);
+
+ }; // END: class RefractiveIndex
+
+} // namespace corsika
+
+#include <corsika/detail/media/UniformRefractiveIndex.inl>
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
new file mode 100644
index 000000000..871ec4148
--- /dev/null
+++ b/tests/media/testRefractiveIndex.cpp
@@ -0,0 +1,87 @@
+/*
+ * (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/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/InhomogeneousMedium.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::units::si;
+
+TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
+
+ CoordinateSystem const& gCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+ // setup our interface types
+ using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+ using AtmModel = UniformRefractiveIndex<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});
+
+ // the refrative index that we use
+ const double n{1.000327};
+
+ // create the atmospheric model
+ AtmModel medium(n, density, protonComposition);
+
+ // and require that it is constant
+ CHECK(n == medium.getRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ CHECK(n == medium.getRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ CHECK(n == medium.getRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ CHECK(n == medium.getRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+
+ // a new refractive index
+ const double n2{2.3472123};
+
+ // update the refractive index of this atmospheric model
+ medium.setRefractiveIndex(n2);
+
+ // check that the returned refractive index is correct
+ CHECK(n2 == medium.getRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ CHECK(n2 == medium.getRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ CHECK(n2 == medium.getRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ CHECK(n2 == medium.getRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+
+ // define our axis vector
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+ // check the density and nuclear composition
+ CHECK(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
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ // and check the integrated grammage
+ CHECK((medium.integratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.arclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+}
--
GitLab