From 12397763d3ce7a1e481b85d5ca9b3919e9666d28 Mon Sep 17 00:00:00 2001
From: Nikos Karastathis <n.karastathis@kit.edu>
Date: Tue, 25 Apr 2023 18:31:49 +0000
Subject: [PATCH] Gladstone-Dale law rerfactive index profile
---
.../media/ExponentialRefractiveIndex.inl | 3 +-
.../media/GladstoneDaleRefractiveIndex.inl | 30 +++++++
corsika/framework/core/PhysicalUnits.hpp | 2 +
.../media/GladstoneDaleRefractiveIndex.hpp | 56 ++++++++++++
corsika/media/UniformRefractiveIndex.hpp | 2 +-
tests/media/testRefractiveIndex.cpp | 85 ++++++++++++++++++-
6 files changed, 173 insertions(+), 5 deletions(-)
create mode 100644 corsika/detail/media/GladstoneDaleRefractiveIndex.inl
create mode 100644 corsika/media/GladstoneDaleRefractiveIndex.hpp
diff --git a/corsika/detail/media/ExponentialRefractiveIndex.inl b/corsika/detail/media/ExponentialRefractiveIndex.inl
index 172631250..5a9d490fd 100644
--- a/corsika/detail/media/ExponentialRefractiveIndex.inl
+++ b/corsika/detail/media/ExponentialRefractiveIndex.inl
@@ -24,7 +24,8 @@ namespace corsika {
, radius_(radius) {}
template <typename T>
- double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
+ inline double ExponentialRefractiveIndex<T>::getRefractiveIndex(
+ Point const& point) const {
return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
}
diff --git a/corsika/detail/media/GladstoneDaleRefractiveIndex.inl b/corsika/detail/media/GladstoneDaleRefractiveIndex.inl
new file mode 100644
index 000000000..5510ead9f
--- /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 referenceRefractiveIndex, Point const point, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , referenceRefractivity_(referenceRefractiveIndex - 1)
+ , referenceInvDensity_(1 / this->getMassDensity(point)) {}
+
+ template <typename T>
+ inline double GladstoneDaleRefractiveIndex<T>::getRefractiveIndex(
+ Point const& point) const {
+ return referenceRefractivity_ * (this->getMassDensity(point) * referenceInvDensity_) +
+ 1.;
+ }
+
+} // namespace corsika
diff --git a/corsika/framework/core/PhysicalUnits.hpp b/corsika/framework/core/PhysicalUnits.hpp
index ba04c8199..5a84d64ca 100644
--- a/corsika/framework/core/PhysicalUnits.hpp
+++ b/corsika/framework/core/PhysicalUnits.hpp
@@ -92,6 +92,8 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::dimensions<-1, 0, 0>, double>;
using InverseTimeType =
phys::units::quantity<phys::units::dimensions<0, 0, -1>, double>;
+ using InverseMassDensityType =
+ phys::units::quantity<phys::units::dimensions<3, -1, 0>, double>;
using InverseGrammageType =
phys::units::quantity<phys::units::dimensions<2, -1, 0>, double>;
using MagneticFluxType =
diff --git a/corsika/media/GladstoneDaleRefractiveIndex.hpp b/corsika/media/GladstoneDaleRefractiveIndex.hpp
new file mode 100644
index 000000000..97f76e774
--- /dev/null
+++ b/corsika/media/GladstoneDaleRefractiveIndex.hpp
@@ -0,0 +1,56 @@
+/*
+ * (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 : public T {
+
+ double const
+ referenceRefractivity_; ///< The constant reference refractive index minus one.
+ InverseMassDensityType const
+ referenceInvDensity_; ///< The constant inverse mass density at reference point.
+
+ 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/corsika/media/UniformRefractiveIndex.hpp b/corsika/media/UniformRefractiveIndex.hpp
index a86afa05f..c5bcebd6c 100644
--- a/corsika/media/UniformRefractiveIndex.hpp
+++ b/corsika/media/UniformRefractiveIndex.hpp
@@ -19,7 +19,7 @@ namespace corsika {
* for all evaluated locations.
*/
template <typename T>
- class UniformRefractiveIndex final : public T {
+ class UniformRefractiveIndex : public T {
double n_; ///< The constant refractive index that we use.
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index 63d1310f3..65e2c749d 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,78 @@ 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 rIndex1{node->getModelProperties().getRefractiveIndex(ref_)};
+ auto const rIndex2{node->getModelProperties().getRefractiveIndex(surface_)};
+
+ CHECK(rIndex1 - n0 == Approx(-0.0002591034));
+ CHECK(rIndex2 - n0 == Approx(0));
+}
\ No newline at end of file
--
GitLab