diff --git a/corsika/detail/media/ExponentialRefractiveIndex.inl b/corsika/detail/media/ExponentialRefractiveIndex.inl
index ac60050c50147ad8a58672ae3459bb1d1192eeb6..17263125073e487400194a0a06ac34199065e4ac 100644
--- a/corsika/detail/media/ExponentialRefractiveIndex.inl
+++ b/corsika/detail/media/ExponentialRefractiveIndex.inl
@@ -15,14 +15,17 @@ namespace corsika {
template <typename T>
template <typename... Args>
ExponentialRefractiveIndex<T>::ExponentialRefractiveIndex(
- double const n0, InverseLengthType const lambda, Args&&... args)
+ double const n0, InverseLengthType const lambda, Point const center,
+ LengthType const radius, Args&&... args)
: T(std::forward<Args>(args)...)
- , n_0(n0)
- , lambda_(lambda) {}
+ , n0_(n0)
+ , lambda_(lambda)
+ , center_(center)
+ , radius_(radius) {}
template <typename T>
double ExponentialRefractiveIndex<T>::getRefractiveIndex(Point const& point) const {
- return n_0 * exp((-lambda_) * point.getCoordinates().getZ());
+ return n0_ * exp((-lambda_) * (distance(point, center_) - radius_));
}
} // namespace corsika
diff --git a/corsika/media/ExponentialRefractiveIndex.hpp b/corsika/media/ExponentialRefractiveIndex.hpp
index 441f7b33fbb14c7ef9568aa5d942845c7591f0cb..02f74dfa88cd21a5bd9c2c4b8b0a6893730d9e3d 100644
--- a/corsika/media/ExponentialRefractiveIndex.hpp
+++ b/corsika/media/ExponentialRefractiveIndex.hpp
@@ -22,8 +22,10 @@ namespace corsika {
template <typename T>
class ExponentialRefractiveIndex : public T {
- double n_0; ///< n0 constant.
+ double n0_; ///< n0 constant.
InverseLengthType lambda_; ///< lambda parameter.
+ LengthType radius_; ///< the planet radius.
+ Point center_; ///< center of the planet.
public:
/**
@@ -37,6 +39,7 @@ namespace corsika {
*/
template <typename... Args>
ExponentialRefractiveIndex(double const n0, InverseLengthType const lambda,
+ Point const center, LengthType const radius,
Args&&... args);
/**
diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index 53c9a694c7c6c3f4f2a05974597856f47bb36140..63d1310f3903671e3e35c50cb5de6b4dc3a04d79 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -10,13 +10,17 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
-#include <corsika/media/FlatExponential.hpp>
+
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMediumModel.hpp>
-#include <corsika/media/InhomogeneousMedium.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/UniformRefractiveIndex.hpp>
#include <corsika/media/ExponentialRefractiveIndex.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <SetupTestTrajectory.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -24,11 +28,13 @@
#include <catch2/catch.hpp>
using namespace corsika;
+template <typename TInterface>
+using MyExtraEnv =
+ ExponentialRefractiveIndex<MediumPropertyModel<UniformMagneticField<TInterface>>>;
TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
logging::set_level(logging::level::info);
- corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
@@ -97,14 +103,13 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
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
+ // setup interface types
using IModelInterface = IRefractiveIndexModel<IMediumModel>;
using AtmModel = ExponentialRefractiveIndex<HomogeneousMedium<IModelInterface>>;
@@ -115,30 +120,33 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
NuclearComposition const protonComposition({Code::Proton}, {1.});
// a new refractive index
- const double n0{2};
- const InverseLengthType lambda{56 / 1_m};
+ const double n0{1};
+ const InverseLengthType lambda{6 / 1_m};
+
+ // the center of the earth
+ Point const center_{gCS, 0_m, 0_m, 0_m};
+ // earth's radius
+ LengthType const radius_{constants::EarthRadius::Mean};
// 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)));
+ AtmModel medium(n0, lambda, center_, constants::EarthRadius::Mean, density,
+ protonComposition);
+ CHECK(n0 - medium.getRefractiveIndex(
+ Point(gCS, 0_m, 0_m, constants::EarthRadius::Mean)) ==
+ Approx(0));
// 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};
+ // distance from the center
+ LengthType const dist_{4_km};
// 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));
+ AtmModel medium_(n0_, lambda_, center_, dist_, density, protonComposition);
+ CHECK(medium_.getRefractiveIndex(Point(gCS, 4_km, 3_km, 0_km)) == Approx(0.3678794412));
- // define our axis vector
+ // define axis vector
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
// check the density and nuclear composition
@@ -146,8 +154,6 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
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();
SpeedType const velocity = 1_m / second;
@@ -169,6 +175,41 @@ TEST_CASE("ExponentialRefractiveIndex w/ Homogeneous medium") {
REQUIRE((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
REQUIRE((medium_.getIntegratedGrammage(track) / (density * length)) == Approx(1));
REQUIRE((medium_.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
- REQUIRE((medium__.getIntegratedGrammage(track) / (density * length)) == Approx(1));
- REQUIRE((medium__.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
+}
+
+TEST_CASE("ExponentialRefractiveIndex 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};
+
+ // another refractive index
+ const double n0{2};
+ const InverseLengthType lambda{1 / 1_km};
+
+ // a reference point to calculate the refractive index there
+ Point const ref_{gCS, 0_m, 0_m, constants::EarthRadius::Mean};
+
+ // setup a 5-layered environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType env;
+
+ create_5layer_atmosphere<EnvironmentInterface, MyExtraEnv>(
+ env, AtmosphereId::LinsleyUSStd, center_, n0, lambda, center_,
+ constants::EarthRadius::Mean, 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));
}