diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index d5e69be6dd383bd5c1ce126331a59c66bbdb3be6..e90160a67a2c74c55c4a40fbad479cfe86d05416 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -23,6 +23,8 @@ set (
ShowerAxis.h
IMagneticFieldModel.h
UniformMagneticField.h
+ IRefractiveIndexModel.h
+ UniformRefractiveIndex.h
)
set (
diff --git a/Environment/IRefractiveIndexModel.h b/Environment/IRefractiveIndexModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..ed8944e670e8f8654242378edcfdeb0e9b6a8120
--- /dev/null
+++ b/Environment/IRefractiveIndexModel.h
@@ -0,0 +1,45 @@
+/*
+ * (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.
+ */
+#pragma once
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::environment {
+
+ /**
+ * An interface for refractive index models.
+ *
+ * This is the base interface for refractive index mixins.
+ *
+ */
+ template <typename Model>
+ class IRefractiveIndexModel : public Model {
+
+ // a type-alias for a dimensionless refractive index
+ using RefractiveIndex = corsika::units::si::RefractiveIndexType;
+
+ public:
+ /**
+ * Evaluate the refractive index at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this point.
+ */
+ virtual RefractiveIndex GetRefractiveIndex(corsika::geometry::Point const&) const = 0;
+
+ /**
+ * A virtual default destructor.
+ */
+ virtual ~IRefractiveIndexModel() = default;
+
+ }; // END: class IRefractiveIndexModel
+
+} // namespace corsika::environment
diff --git a/Environment/UniformRefractiveIndex.h b/Environment/UniformRefractiveIndex.h
new file mode 100644
index 0000000000000000000000000000000000000000..a64f07c4f3a33293ff9b922b949deb2fd72292c2
--- /dev/null
+++ b/Environment/UniformRefractiveIndex.h
@@ -0,0 +1,66 @@
+/*
+ * (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.
+ */
+#pragma once
+
+#include <corsika/environment/IRefractiveIndexModel.h>
+
+namespace corsika::environment {
+
+ /**
+ * A uniform refractive index.
+ *
+ * This class returns the same refractive index
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class UniformRefractiveIndex : public T {
+
+ // a type-alias for a dimensionless refractive index
+ using RefractiveIndex = corsika::units::si::RefractiveIndexType;
+
+ RefractiveIndex 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(RefractiveIndex const n, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , n_(n) {}
+
+ /**
+ * Evaluate the refractive index at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this point.
+ */
+ RefractiveIndex GetRefractiveIndex(
+ corsika::geometry::Point const&) const final override {
+ return n_;
+ }
+
+ /**
+ * Set the refractive index returned by this instance.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this location.
+ */
+ void SetRefractiveIndex(RefractiveIndex const& n) { n_ = n; }
+
+ }; // END: class RefractiveIndex
+
+} // namespace corsika::environment
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 55ef5fba4f2473e6e59141953c1b191d2ceadaec..b30e5ffd42e283a14475dd708934372bc9a9df31 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -13,12 +13,17 @@
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/IRefractiveIndexModel.h>
#include <corsika/environment/InhomogeneousMedium.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/LinearApproximationIntegrator.h>
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/SlidingPlanarExponential.h>
+<<<<<<< HEAD
#include <corsika/environment/UniformMagneticField.h>
+=======
+#include <corsika/environment/UniformRefractiveIndex.h>
+>>>>>>> a9fca8f... Initial refractive index interface and models with unittests.
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Line.h>
#include <corsika/geometry/RootCoordinateSystem.h>
@@ -249,8 +254,14 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
const auto density{19.2_g / cube(1_cm)};
// create our atmospheric model
+<<<<<<< HEAD
AtmModel medium(B0, density, protonComposition);
+ // and test at several locations
+ REQUIRE(B0 == medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)));
+ REQUIRE(B0 == medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)));
+ REQUIRE(B0 == medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)));
+
// create a new magnetic field vector
QuantityVector B1(23_T, 57_T, -4_T);
@@ -279,4 +290,99 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
// and check the integrated grammage
REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+=======
+ AtmModel medium(n, 19.2_g / cube(1_cm), protonComposition);
+
+ // and require that it is constant
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+
+ // a new refractive index
+ RefractiveIndexType n2 = 2.3472123__;
+
+ // update the refractive index of this atmospheric model
+ medium.SetRefractiveIndex(n2);
+
+ // check that the returned refractive index is correct
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+>>>>>>> 12a0be6... Add test for SetRefractiveIndex
+}
+
+TEST_CASE("UniformRefractiveIndex w/ FlatExponential") {
+
+ // setup our interface types
+ using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+ using AtmModel = UniformRefractiveIndex<FlatExponential<IModelInterface>>;
+
+ // the composition we use for the homogenous medium
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ // define our axis vector
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+ // the parameters of our exponential model
+ LengthType const lambda = 3_m;
+ auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+
+ // the refractive index we want returned everywhere
+ RefractiveIndexType n = 1.000327__;
+
+ // create our atmospheric model
+ AtmModel medium(n, gOrigin, axis, rho0, lambda, protonComposition);
+
+ // check that the returned refractive index is correct
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+
+<<<<<<< HEAD
+ // the refractive index we want returned everywhere
+ RefractiveIndexType n = 1.1234__;
+
+ // create our atmospheric model
+ AtmModel const medium(n, 19.2_g / cube(1_cm), protonComposition);
+
+ // and require that it is constant
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+
+ // check the density and nuclear composition
+ 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
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ // and check the integrated grammage
+ REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+=======
+ // a new refractive index
+ RefractiveIndexType n2 = 1.123456__;
+
+ // update the refractive index of this atmospheric model
+ medium.SetRefractiveIndex(n2);
+
+ // check that the returned refractive index is correct
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+>>>>>>> 12a0be6... Add test for SetRefractiveIndex
}
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index e701928f8cc38b3fca6e91a2de3bbee12e02861a..35d4cd7a1b16fb218bb305b1edfbbd21779ae660 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -69,6 +69,12 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
using InverseGrammageType =
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
+ using RefractiveIndexType = phys::units::quantity<dimensionless_d, double>;
+
+ // define a new literal to construct dimensionless types
+ // This is a 'double underscore' literal, So 3__ is a
+ // dimensionless quantity = 3
+ QUANTITY_DEFINE_LITERALS(_, dimensionless_d)
namespace detail {
template <int N, typename T>