added and tested Exponential Refractive Index class in media

corsika/detail/media/ExponentialRefractiveIndex.inl

+///*
+// * (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 <bits/stdc++.h>
+//#include <corsika/media/IRefractiveIndexModel.hpp>
+//
+//namespace corsika {
+//
+//  template <typename T>
+//  template <typename... Args>
+//  ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(double const n0,
+//                             InverseLengthType const lambda, Args&&... args)
+//      : T(std::forward<Args>(args)...)
+//      , n_0(n0)
+//      , lambda_(lambda) {}
+//
+//  template <typename T>
+//  double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
+//  //TODO: THIS METHOD CURRENTLY ONLY USES THE Z-COORDINATE.
+//  //NEED TO THINK IT FOR FUTURE WORK ON ARBITRARY GEOMETRIES.
+//  return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+//  }
+//
+//} // namespace corsika

corsika/media/ExponentialRefractiveIndex.hpp

+/*
+ * (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 <bits/stdc++.h>
+#include <corsika/media/IRefractiveIndexModel.hpp>
+
+namespace corsika {
+
+  /**
+   * An exponential refractive index.
+   *
+   * This class returns the value of an exponential refractive index
+   * for all evaluated locations.
+   *
+   */
+  template <typename T>
+  class ExponentialRefractiveIndex : public T {
+
+    double n_0; ///< n0 constant.
+    InverseLengthType lambda_; ///< lambda parameter.
+
+  public:
+    /**
+     * Construct an ExponentialRefractiveIndex.
+     *
+     * This is initialized with two parameters n_0 and lambda
+     * and returns the value of the exponential refractive index
+     * at a given point location.
+     *
+     * @param field    The refractive index to return to a given point.
+     */
+    template <typename... Args>
+    ExponentialRefractiveIndex(double const n0,
+                               InverseLengthType const lambda, Args&&... args)
+        : T(std::forward<Args>(args)...)
+        , n_0(n0)
+        , lambda_(lambda) {}
+    /**
+     * 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 final override
+    {
+      //TODO: THIS METHOD CURRENTLY ONLY USES THE Z-COORDINATE.
+      //NEED TO THINK IT FOR FUTURE WORK ON ARBITRARY GEOMETRIES.
+      return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+    }
+
+  }; // END: class ExponentialRefractiveIndex
+
+} // namespace corsika
\ No newline at end of file

tests/media/testRefractiveIndex.cpp

@@ -16,6 +16,7 @@
 #include <corsika/media/InhomogeneousMedium.hpp>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/media/UniformRefractiveIndex.hpp>
+#include <corsika/media/ExponentialRefractiveIndex.hpp>
 
 #include <SetupTestTrajectory.hpp>
 
@@ -89,3 +90,82 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
   CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
   CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
 }
+
+
+TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
+
+  logging::set_level(logging::level::info);
+  corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
+
+  // get a CS and a point
+  CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+
+  Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+  // setup our interface types
+  using IModelInterface = IRefractiveIndexModel<IMediumModel>;
+  using AtmModel = ExponentialRefractiveIndex<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});
+
+  // a new refractive index
+  const double n0{2};
+  const InverseLengthType lambda{56 / 1_m};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium(n0, lambda, density, protonComposition);
+  CHECK(n0 == medium.getRefractiveIndex(Point(gCS, 10_m, 14_m, 0_m)));
+
+  // another refractive index
+  const double n0_{1};
+  const InverseLengthType lambda_{1 / 1_km};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium_(n0_, lambda_, density, protonComposition);
+  CHECK(medium_.getRefractiveIndex(Point(gCS, 12_m, 56_m, 1_km)) == Approx(0.3678794));
+
+  // another refractive index
+  const double n0__{4};
+  const InverseLengthType lambda__{2 / 1_m};
+
+  // create the atmospheric model and check refractive index
+  AtmModel medium__(n0__, lambda__, density, protonComposition);
+  CHECK(medium__.getRefractiveIndex(Point(gCS, 0_m, 0_m, 35_km)) == Approx(0));
+
+  // define our axis vector
+  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+  // check the density and nuclear composition
+  REQUIRE(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+  medium.getNuclearComposition();
+  REQUIRE(density == medium_.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+  medium_.getNuclearComposition();
+  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
+  setup::Trajectory const track =
+      setup::testing::make_track<setup::Trajectory>(line, tEnd);
+//  // and the associated trajectory
+//  Trajectory<Line> const trajectory(line, tEnd);
+
+  // and check the integrated grammage
+  REQUIRE((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+  REQUIRE((medium_.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+  REQUIRE((medium__.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
+  REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+}