diff --git a/tests/media/testRefractiveIndex.cpp b/tests/media/testRefractiveIndex.cpp
index 53c9a694c7c6c3f4f2a05974597856f47bb36140..aace8c2c0c05249f78f9797b027917f4f54c15e9 100644
--- a/tests/media/testRefractiveIndex.cpp
+++ b/tests/media/testRefractiveIndex.cpp
@@ -10,10 +10,12 @@
#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>
@@ -24,11 +26,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 +101,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 +118,30 @@ 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 +149,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 +170,55 @@ 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/ 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 realistic environment
+ using EnvironmentInterface =
+ IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+ using EnvType = Environment<EnvironmentInterface>;
+ EnvType env;
+
+ auto builder = make_layered_spherical_atmosphere_builder<
+ EnvironmentInterface, MyExtraEnv>::create(center_,
+ constants::EarthRadius::Mean, n0, lambda,
+ center_, constants::EarthRadius::Mean,
+ Medium::AirDry1Atm,
+ MagneticFieldVector{gCS, 10_uT,
+ 0_T, 0_T});
+ builder.setNuclearComposition(
+ {{Code::Nitrogen, Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}});
+
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
+ builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+ builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
+ builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
+ builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
+ builder.assemble(env);
+
+ // 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));
}