Gladstone-Dale law rerfactive index profile

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_));
   }
 

corsika/detail/media/GladstoneDaleRefractiveIndex.inl

+/*
+ * (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

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 =

corsika/media/GladstoneDaleRefractiveIndex.hpp

+/*
+ * (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>

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.
 

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