From 19be0728b5d8bb1f2c624f6657f4a2ea144ba65e Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Fri, 20 Jan 2023 12:18:19 +0100
Subject: [PATCH] refractive index using Gladstone-Dale law & some first tests
---
.../media/GladstoneDaleRefractiveIndex.inl | 30 +++++++
.../media/GladstoneDaleRefractiveIndex.hpp | 55 ++++++++++++
tests/media/testRefractiveIndex.cpp | 84 ++++++++++++++++++-
3 files changed, 166 insertions(+), 3 deletions(-)
create mode 100644 corsika/detail/media/GladstoneDaleRefractiveIndex.inl
create mode 100644 corsika/media/GladstoneDaleRefractiveIndex.hpp
diff --git a/corsika/detail/media/GladstoneDaleRefractiveIndex.inl b/corsika/detail/media/GladstoneDaleRefractiveIndex.inl
new file mode 100644
index 000000000..d0a13a981
--- /dev/null
+++ b/corsika/detail/media/GladstoneDaleRefractiveIndex.inl
@@ -0,0 +1,30 @@
+/*
+* (c) Copyright 2023 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>
+ inline GladstoneDaleRefractiveIndex<T>::GladstoneDaleRefractiveIndex(double const seaLevelRefractiveIndex,
+ Point const point, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , seaLevelRefractiveIndex_(seaLevelRefractiveIndex)
+ , earthSurface_(point)
+ , seaLevelDensity_(this->getMassDensity(earthSurface_)) {}
+// , ratio_((1. - seaLevelRefractiveIndex_) / (seaLevelDensity_ * cube(1_m) / 1_kg)) {}
+
+ template <typename T>
+ inline double GladstoneDaleRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
+ return (seaLevelRefractiveIndex_ - 1.) * (this->getMassDensity(point) / seaLevelDensity_) + 1.;
+ }
+
+} // namespace corsika
diff --git a/corsika/media/GladstoneDaleRefractiveIndex.hpp b/corsika/media/GladstoneDaleRefractiveIndex.hpp
new file mode 100644
index 000000000..767351d8f
--- /dev/null
+++ b/corsika/media/GladstoneDaleRefractiveIndex.hpp
@@ -0,0 +1,55 @@
+/*
+* (c) Copyright 2023 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 tabulated refractive index.
+ *
+ * This class returns the value of refractive index
+ * for a specific height bin.
+ */
+ template <typename T>
+ class GladstoneDaleRefractiveIndex final : public T {
+
+ double const seaLevelRefractiveIndex_; ///< The constant refractive index at sea level.
+ MassDensityType const seaLevelDensity_; ///< The constant mass density at sea level.
+ Point const earthSurface_; ///< A point at earth's surface.
+// MassDensityType ratio_; ///< (1 - seaLevelRefractiveIndex_) / seaLevelDensity_.
+
+ public:
+ /**
+ * Construct a GladstoneDaleRefractiveIndex.
+ *
+ * This is initialized with a fixed refractive index
+ * at sea level and uses the Gladstone-Dale law to
+ * calculate the refractive index at a given point.
+ *
+ * @param seaLevelRefractiveIndex The refractive index at sea level.
+ * @param point AA point at earth's surface.
+ */
+ template <typename... Args>
+ GladstoneDaleRefractiveIndex(double const seaLevelRefractiveIndex, Point const point, 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& point) const override;
+
+ }; // END: class GladstoneDaleRefractiveIndex
+
+} // namespace corsika
+
+#include <corsika/detail/media/GladstoneDaleRefractiveIndex.inl>
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index 63d1310f3..3a7bc077b 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -20,6 +20,7 @@
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
#include <corsika/media/ExponentialRefractiveIndex.hpp>
+#include <corsika/media/GladstoneDaleRefractiveIndex.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <SetupTestTrajectory.hpp>
@@ -31,8 +32,11 @@ using namespace corsika;
template <typename TInterface>
using MyExtraEnv =
ExponentialRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
+template <typename TInterface2>
+using MyExtraEnv2 =
+ GladstoneDaleRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface2>>>;
-TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
+TEST_CASE("UniformRefractiveIndex w/ Homogeneous medium") {
logging::set_level(logging::level::info);
@@ -40,7 +44,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
- // setup our interface types
+ // set up our interface types
using IModelInterface = IRefractiveIndexModel<IMediumModel>;
using AtmModel = UniformRefractiveIndex<HomogeneousMedium<IModelInterface>>;
@@ -50,7 +54,7 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
// the composition we use for the homogenous medium
NuclearComposition const protonComposition({Code::Proton}, {1.});
- // the refrative index that we use
+ // the refractive index that we use
const double n{1.000327};
// create the atmospheric model
@@ -213,3 +217,77 @@ TEST_CASE("ExponentialRefractiveIndex w/ 5-layered atmosphere") {
CHECK(rIndex - n0 == Approx(0));
}
+
+TEST_CASE("GladstoneDaleRefractiveIndex w/ Homogeneous medium") {
+
+ logging::set_level(logging::level::info);
+
+ // get a CS and a point
+ CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+ Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+ // setup interface types
+ using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+ using AtmModel = GladstoneDaleRefractiveIndex<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({Code::Proton}, {1.});
+
+ // the refractive index at sea level
+ const double n0{1.000327};
+
+ // a point at the surface of the earth
+ Point const surface_{gCS, 0_m, 0_m, constants::EarthRadius::Mean};
+
+ // a random point in the atmosphere
+ Point const p1_{gCS, 1_km, 1_km, constants::EarthRadius::Mean + 10_km};
+
+ // create the atmospheric model and check refractive index
+ AtmModel medium(n0, surface_, density, protonComposition);
+
+ CHECK(n0 - medium.getRefractiveIndex(surface_) == Approx(0));
+ CHECK(n0 - medium.getRefractiveIndex(p1_) == Approx(0));
+
+}
+
+TEST_CASE("GladstoneDaleRefractiveIndex w/ 5-layered atmosphere") {
+
+ logging::set_level(logging::level::info);
+
+ // get a CS
+ CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+ // the center of the earth
+ Point const center_{gCS, 0_m, 0_m, 0_m};
+
+ // a point at the surface of the earth
+ Point const surface_{gCS, 0_m, 0_m, constants::EarthRadius::Mean};
+
+ // the refractive index at sea level
+ const double n0{1.000327};
+
+ // a reference point to calculate the refractive index there
+ Point const ref_{gCS, 0_km, 0_km, constants::EarthRadius::Mean + 10_km};
+
+ // setup a 5-layered environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType env;
+
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv2>(
+ env, AtmosphereId::LinsleyUSStd, center_, n0, surface_,
+ Medium::AirDry1Atm, MagneticFieldVector{gCS, 0_T, 50_uT, 0_T});
+
+ // get the universe for this environment
+ auto const* const universe{env.getUniverse().get()};
+ auto const* node{universe->getContainingNode(ref_)};
+ // get the refractive index
+ auto const rIndex{node->getModelProperties().getRefractiveIndex(ref_)};
+
+// CHECK(rIndex - n0 == Approx(0));
+}
\ No newline at end of file
--
GitLab