diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 64ed627e03c8c7687a2b75330dd50d46e1168627..b8551c87ac0e730014c65439825252b6c527866f 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -45,7 +45,6 @@ using namespace corsika::process;
using namespace corsika::units;
using namespace corsika::particles;
using namespace corsika::random;
-using namespace corsika::setup;
using namespace corsika::geometry;
using namespace corsika::environment;
@@ -68,27 +67,30 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("cascade");
// setup environment, geometry
- using EnvType = environment::Environment<setup::IEnvironmentModel>;
- EnvType env;
+ setup::Environment env;
auto& universe = *(env.GetUniverse());
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
- auto outerMedium = EnvType::CreateNode<Sphere>(
+ auto outerMedium = setup::EnvironmentType::CreateNode<Sphere>(
Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
+
+ using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
+
// fraction of oxygen
const float fox = 0.20946;
auto const props =
outerMedium
- ->SetModelProperties<environment::HomogeneousMedium<setup::IEnvironmentModel>>(
- 1_kg / (1_m * 1_m * 1_m),
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Nitrogen,
- particles::Code::Oxygen},
- std::vector<float>{1.f - fox, fox}));
-
- auto innerMedium = EnvType::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
+ ->SetModelProperties<EnvironmentModel>(environment::EMediumType::eAir,
+ Vector(rootCS, 0_T, 0_T, 0_T),
+ 1_kg / (1_m * 1_m * 1_m),
+ environment::NuclearComposition(
+ std::vector<particles::Code>{particles::Code::Nitrogen,
+ particles::Code::Oxygen},
+ std::vector<float>{1.f - fox, fox}));
+
+ auto innerMedium = setup::Environment::CreateNode<Sphere>(Point{rootCS, 0_m, 0_m, 0_m}, 5000_m);
innerMedium->SetModelProperties(props);
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index f0adcf74a8fcbd55abf8eaec6b75a42700630f6d..9edd10b54c3a6daa25a00920c58a0fe6a6c2ad5f 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -25,6 +25,9 @@ set (
UniformMagneticField.h
IRefractiveIndexModel.h
UniformRefractiveIndex.h
+ MediumTypes.h
+ IMediumTypeModel.h
+ UniformMediumType.h
)
set (
diff --git a/Environment/Environment.h b/Environment/Environment.h
index ef86670e8e35c1d772ecc88e80b6bfe3c1b84358..774b7492761469affbfc2d6bf326cac5501dbbf2 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -38,8 +38,6 @@ namespace corsika::environment {
, fUniverse(std::make_unique<BaseNodeType>(
std::make_unique<Universe>(fCoordinateSystem))) {}
- // using IEnvironmentModel = corsika::setup::IEnvironmentModel;
-
auto& GetUniverse() { return fUniverse; }
auto const& GetUniverse() const { return fUniverse; }
@@ -61,7 +59,4 @@ namespace corsika::environment {
typename BaseNodeType::VTNUPtr fUniverse;
};
- // using SetupBaseNodeType = VolumeTreeNode<corsika::setup::IEnvironmentModel>;
- // using SetupEnvironment = Environment<corsika::setup::IEnvironmentModel>;
-
} // namespace corsika::environment
diff --git a/Environment/IMediumType.h b/Environment/IMediumType.h
new file mode 100644
index 0000000000000000000000000000000000000000..abf64479d9c24a267ca84be255da5f3f0c869cfe
--- /dev/null
+++ b/Environment/IMediumType.h
@@ -0,0 +1,41 @@
+/*
+ * (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 <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::environment {
+
+ /**
+ * An interface for type of media, needed e.g. to determine energy losses
+ *
+ * This is the base interface for media types.
+ *
+ */
+ template <typename Model>
+ class IMediumType : public Model {
+
+ public:
+ /**
+ * Evaluate the medium type at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The media type
+ */
+ virtual double medium_type(corsika::geometry::Point const&) const = 0;
+
+ /**
+ * A virtual default destructor.
+ */
+ virtual ~IMediumType() = default;
+
+ }; // END: class IMediumType
+
+} // namespace corsika::environment
diff --git a/Environment/IMediumTypeModel.h b/Environment/IMediumTypeModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..3e8d0471618a9463c08ebcb79e44ad26076b0f8d
--- /dev/null
+++ b/Environment/IMediumTypeModel.h
@@ -0,0 +1,43 @@
+/*
+ * (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 <corsika/environment/MediumTypes.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::environment {
+
+ /**
+ * An interface for type of media, needed e.g. to determine energy losses
+ *
+ * This is the base interface for media types.
+ *
+ */
+ template <typename Model>
+ class IMediumTypeModel : public Model {
+
+ public:
+ /**
+ * Evaluate the medium type at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The media type
+ */
+ virtual EMediumType medium_type(corsika::geometry::Point const&) const = 0;
+
+ /**
+ * A virtual default destructor.
+ */
+ virtual ~IMediumTypeModel() = default;
+
+ }; // END: class IMediumTypeModel
+
+} // namespace corsika::environment
diff --git a/Environment/MediumTypes.h b/Environment/MediumTypes.h
new file mode 100644
index 0000000000000000000000000000000000000000..d07e68aeccc979077ca241ce3167cdd8e6df3eaa
--- /dev/null
+++ b/Environment/MediumTypes.h
@@ -0,0 +1,15 @@
+/*
+ * (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
+
+namespace corsika::environment {
+
+ enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
+
+}
diff --git a/Environment/UniformMediaType.h b/Environment/UniformMediaType.h
new file mode 100644
index 0000000000000000000000000000000000000000..fcd7e35f0357c0e085add466625342823a571768
--- /dev/null
+++ b/Environment/UniformMediaType.h
@@ -0,0 +1,61 @@
+/*
+ * (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 <corsika/environment/IMediaTypeModel.h>
+
+namespace corsika::environment {
+
+ /**
+ * A uniform refractive index.
+ *
+ * This class returns the same refractive index
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class UniformMediaType : public T {
+
+ EMediaType double n_; ///< The constant refractive index that we use.
+
+ public:
+ /**
+ * Construct a UniformMediaType.
+ *
+ * This is initialized with a fixed refractive index
+ * and returns this refractive index at all locations.
+ *
+ * @param field The refractive index to return everywhere.
+ */
+ template <typename... Args>
+ UniformMediaType(double const n, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , n_(n) {}
+
+ /**
+ * Evaluate the refractive index at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this point.
+ */
+ double GetMediaType(corsika::geometry::Point const&) const final override {
+ return n_;
+ }
+
+ /**
+ * Set the refractive index returned by this instance.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this location.
+ */
+ void SetMediaType(double const& n) { n_ = n; }
+
+ }; // END: class MediaType
+
+} // namespace corsika::environment
diff --git a/Environment/UniformMediumType.h b/Environment/UniformMediumType.h
new file mode 100644
index 0000000000000000000000000000000000000000..5529b8a50a8b0c39d68c4b5859bc30d4929e949b
--- /dev/null
+++ b/Environment/UniformMediumType.h
@@ -0,0 +1,62 @@
+/*
+ * (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 <corsika/environment/IMediumTypeModel.h>
+#include <corsika/environment/MediumTypes.h>
+
+namespace corsika::environment {
+
+ /**
+ * A uniform refractive index.
+ *
+ * This class returns the same refractive index
+ * for all evaluated locations.
+ *
+ */
+ template <typename T>
+ class UniformMediumType : public T {
+
+ EMediumType type_; ///< The constant medium type
+
+ public:
+ /**
+ * Construct a UniformMediumType.
+ *
+ * This is initialized with a fixed medium type
+ * and returns this at all locations.
+ *
+ * @param field The medium type to return everywhere.
+ */
+ template <typename... Args>
+ UniformMediumType(EMediumType const type, Args&&... args)
+ : T(std::forward<Args>(args)...)
+ , type_(type) {}
+
+ /**
+ * Evaluate the refractive index at a given location.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this point.
+ */
+ EMediumType medium_type(corsika::geometry::Point const&) const final override {
+ return type_;
+ }
+
+ /**
+ * Set the refractive index returned by this instance.
+ *
+ * @param point The location to evaluate at.
+ * @returns The refractive index at this location.
+ */
+ void set_medium_type(EMediumType const& type) { type_ = type; }
+
+ }; // END: class MediumType
+
+} // namespace corsika::environment
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 088dacb1c281befe76301145f1e5de4992e7d085..7394abf42b66e03f6489e8ca9e3885a599c92b65 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -11,6 +11,7 @@
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/IMediumTypeModel.h>
#include <corsika/environment/IRefractiveIndexModel.h>
#include <corsika/environment/InhomogeneousMedium.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
@@ -18,6 +19,7 @@
#include <corsika/environment/NuclearComposition.h>
#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/UniformMagneticField.h>
+#include <corsika/environment/UniformMediumType.h>
#include <corsika/environment/UniformRefractiveIndex.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Line.h>
@@ -61,8 +63,8 @@ TEST_CASE("FlatExponential") {
gCS, {20_cm / second, 0_m / second, 0_m / second}));
Trajectory<Line> const trajectory(line, tEnd);
- REQUIRE((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
- REQUIRE((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
+ CHECK((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
+ CHECK((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
}
SECTION("vertical") {
@@ -72,8 +74,8 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2 * lambda;
GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
- REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
- REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
+ CHECK((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
+ CHECK((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
SECTION("escape grammage") {
@@ -83,9 +85,9 @@ TEST_CASE("FlatExponential") {
GrammageType const escapeGrammage = rho0 * lambda;
- REQUIRE(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
- REQUIRE(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
- std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
+ CHECK(trajectory.NormalizedDirection().dot(axis).magnitude() < 0);
+ CHECK(medium.ArclengthFromGrammage(trajectory, 1.2 * escapeGrammage) ==
+ std::numeric_limits<typename GrammageType::value_type>::infinity() * 1_m);
}
SECTION("inclined") {
@@ -96,8 +98,8 @@ TEST_CASE("FlatExponential") {
LengthType const length = 2 * lambda;
GrammageType const exact =
rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
- REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
- REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
+ CHECK((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
+ CHECK((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
}
@@ -170,9 +172,9 @@ TEST_CASE("InhomogeneousMedium") {
DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
SECTION("DensityFunction") {
- REQUIRE(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
- Approx(1));
- REQUIRE(rho.EvaluateAt(gOrigin) == e(gOrigin));
+ CHECK(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+ Approx(1));
+ CHECK(rho.EvaluateAt(gOrigin) == e(gOrigin));
}
auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
@@ -184,18 +186,18 @@ TEST_CASE("InhomogeneousMedium") {
InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
SECTION("Integration") {
- REQUIRE(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
- Approx(1).epsilon(1e-2));
- REQUIRE(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
- exactLength(exactGrammage(l)) ==
- Approx(1).epsilon(1e-2));
- REQUIRE(rho.MaximumLength(trajectory, 1e-2) >
- l); // todo: write reasonable test when implementation is working
-
- REQUIRE(rho.IntegrateGrammage(trajectory, l) ==
- inhMedium.IntegratedGrammage(trajectory, l));
- REQUIRE(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
- inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
+ CHECK(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
+ Approx(1).epsilon(1e-2));
+ CHECK(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
+ exactLength(exactGrammage(l)) ==
+ Approx(1).epsilon(1e-2));
+ CHECK(rho.MaximumLength(trajectory, 1e-2) >
+ l); // todo: write reasonable test when implementation is working
+
+ CHECK(rho.IntegrateGrammage(trajectory, l) ==
+ inhMedium.IntegratedGrammage(trajectory, l));
+ CHECK(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
+ inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
}
}
@@ -208,27 +210,27 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
builder.addLinearLayer(2_km, 20_km);
builder.addLinearLayer(3_km, 30_km);
- REQUIRE(builder.size() == 3);
+ CHECK(builder.size() == 3);
auto const builtEnv = builder.assemble();
auto const& univ = builtEnv.GetUniverse();
- REQUIRE(builder.size() == 0);
+ CHECK(builder.size() == 0);
auto const R = builder.getEarthRadius();
- REQUIRE(univ->GetChildNodes().size() == 1);
-
- REQUIRE(univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 35_km)) == univ.get());
- REQUIRE(dynamic_cast<Sphere const&>(
- univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 8_km))->GetVolume())
- .GetRadius() == R + 10_km);
- REQUIRE(dynamic_cast<Sphere const&>(
- univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 12_km))->GetVolume())
- .GetRadius() == R + 20_km);
- REQUIRE(dynamic_cast<Sphere const&>(
- univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 24_km))->GetVolume())
- .GetRadius() == R + 30_km);
+ CHECK(univ->GetChildNodes().size() == 1);
+
+ CHECK(univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 35_km)) == univ.get());
+ CHECK(dynamic_cast<Sphere const&>(
+ univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 8_km))->GetVolume())
+ .GetRadius() == R + 10_km);
+ CHECK(dynamic_cast<Sphere const&>(
+ univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 12_km))->GetVolume())
+ .GetRadius() == R + 20_km);
+ CHECK(dynamic_cast<Sphere const&>(
+ univ->GetContainingNode(Point(gCS, 0_m, 0_m, R + 24_km))->GetVolume())
+ .GetRadius() == R + 30_km);
}
TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
@@ -251,13 +253,13 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
AtmModel medium(B0, density, protonComposition);
// and test at several locations
- REQUIRE(B0.GetComponents(gCS) ==
- medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
- REQUIRE(
+ CHECK(B0.GetComponents(gCS) ==
+ medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
+ CHECK(
B0.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).GetComponents(gCS));
- REQUIRE(B0.GetComponents(gCS) ==
- medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
+ CHECK(B0.GetComponents(gCS) ==
+ medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
// create a new magnetic field vector
Vector B1(gCS, 23_T, 57_T, -4_T);
@@ -266,16 +268,16 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
medium.SetMagneticField(B1);
// and test at several locations
- REQUIRE(B1.GetComponents(gCS) ==
- medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
- REQUIRE(
+ CHECK(B1.GetComponents(gCS) ==
+ medium.GetMagneticField(Point(gCS, -10_m, 4_m, 35_km)).GetComponents(gCS));
+ CHECK(
B1.GetComponents(gCS) ==
medium.GetMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).GetComponents(gCS));
- REQUIRE(B1.GetComponents(gCS) ==
- medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
+ CHECK(B1.GetComponents(gCS) ==
+ medium.GetMagneticField(Point(gCS, 0_m, 0_m, 0_m)).GetComponents(gCS));
// check the density and nuclear composition
- REQUIRE(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
+ CHECK(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.GetNuclearComposition();
// create a line of length 1 m
@@ -289,8 +291,8 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
Trajectory<Line> const trajectory(line, tEnd);
// and check the integrated grammage
- REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
- REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+ CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
@@ -313,10 +315,10 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
AtmModel medium(n, density, protonComposition);
// and require that it is constant
- REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
- REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
- REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
- REQUIRE(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+ CHECK(n == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ CHECK(n == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ CHECK(n == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ CHECK(n == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
// a new refractive index
const double n2{2.3472123};
@@ -325,16 +327,75 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
medium.SetRefractiveIndex(n2);
// check that the returned refractive index is correct
- REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
- REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
- REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
- REQUIRE(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+ CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, -10_m, 4_m, 35_km)));
+ CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, +210_m, 0_m, 7_km)));
+ CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, 0_m, 0_m, 0_km)));
+ CHECK(n2 == medium.GetRefractiveIndex(Point(gCS, 100_km, 400_km, 350_km)));
+
+ // define our axis vector
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+
+ // check the density and nuclear composition
+ CHECK(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
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ // and check the integrated grammage
+ CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+}
+
+TEST_CASE("UniformMediumType w/ Homogeneous") {
+
+ // setup our interface types
+ using IModelInterface = IMediumTypeModel<IMediumModel>;
+ using AtmModel = UniformMediumType<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});
+
+ // the refrative index that we use
+ const EMediumType type = EMediumType::eAir;
+
+ // create the atmospheric model
+ AtmModel medium(type, density, protonComposition);
+
+ // and require that it is constant
+ CHECK(type == medium.medium_type(Point(gCS, -10_m, 4_m, 35_km)));
+ CHECK(type == medium.medium_type(Point(gCS, +210_m, 0_m, 7_km)));
+ CHECK(type == medium.medium_type(Point(gCS, 0_m, 0_m, 0_km)));
+ CHECK(type == medium.medium_type(Point(gCS, 100_km, 400_km, 350_km)));
+
+ // a new refractive index
+ const EMediumType type2 = EMediumType::eRock;
+
+ // update the refractive index of this atmospheric model
+ medium.set_medium_type(type2);
+
+ // check that the returned refractive index is correct
+ CHECK(type2 == medium.medium_type(Point(gCS, -10_m, 4_m, 35_km)));
+ CHECK(type2 == medium.medium_type(Point(gCS, +210_m, 0_m, 7_km)));
+ CHECK(type2 == medium.medium_type(Point(gCS, 0_m, 0_m, 0_km)));
+ CHECK(type2 == medium.medium_type(Point(gCS, 100_km, 400_km, 350_km)));
// 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)));
+ CHECK(density == medium.GetMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.GetNuclearComposition();
// create a line of length 1 m
@@ -348,6 +409,6 @@ TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
Trajectory<Line> const trajectory(line, tEnd);
// and check the integrated grammage
- REQUIRE((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
- REQUIRE((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
+ CHECK((medium.IntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
+ CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
diff --git a/Processes/CONEXSourceCut/testCONEXSourceCut.cc b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
index 82d208521e4b2fd86aca6349305c5dbaec609e58..f70d54789319c8546d5bacd6da6b5e315c85e862 100644
--- a/Processes/CONEXSourceCut/testCONEXSourceCut.cc
+++ b/Processes/CONEXSourceCut/testCONEXSourceCut.cc
@@ -6,18 +6,28 @@
* the license.
*/
+#include <corsika/setup/SetupEnvironment.h>
+
#include <corsika/environment/Environment.h>
#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/environment/UniformMediumType.h>
+#include <corsika/environment/UniformMagneticField.h>
+
#include <corsika/geometry/Point.h>
#include <corsika/geometry/RootCoordinateSystem.h>
#include <corsika/geometry/Vector.h>
+
#include <corsika/particles/ParticleProperties.h>
+
#include <corsika/process/conex_source_cut/CONEXSourceCut.h>
#include <corsika/process/sibyll/Interaction.h>
#include <corsika/process/sibyll/NuclearInteraction.h>
+
#include <corsika/random/RNGManager.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
+
#include <catch2/catch.hpp>
using namespace corsika;
@@ -32,11 +42,14 @@ TEST_CASE("CONEXSourceCut") {
feenableexcept(FE_INVALID);
// setup environment, geometry
- using EnvType = Environment<setup::IEnvironmentModel>;
- EnvType env;
+ setup::Environment env;
+ auto& universe = *(env.GetUniverse());
+ using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::InhomogeneousMedium<setup::IEnvironment>>>;
+
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
environment::LayeredSphericalAtmosphereBuilder builder{center, conex::earthRadius};
+
builder.setNuclearComposition(
{{particles::Code::Nitrogen, particles::Code::Oxygen},
{0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
diff --git a/Processes/Proposal/CMakeLists_PROPOSAL.txt b/Processes/Proposal/CMakeLists_PROPOSAL.txt
index ad1b75b76d839ad1ac2173f4dbc1a6917d5a41c5..7045ae65db2b0a81b82bfd38ba643d1a9056590e 100644
--- a/Processes/Proposal/CMakeLists_PROPOSAL.txt
+++ b/Processes/Proposal/CMakeLists_PROPOSAL.txt
@@ -55,8 +55,8 @@ else (IN_CORSIKA8)
endif (IN_CORSIKA8)
add_library(PROPOSAL::PROPOSAL ALIAS PROPOSAL)
-target_compile_features(PROPOSAL PUBLIC cxx_std_11)
-set_target_properties(PROPOSAL PROPERTIES CXX_EXTENSIONS OFF)
+#target_compile_features(PROPOSAL PUBLIC cxx_std_11)
+#set_target_properties(PROPOSAL PROPERTIES CXX_EXTENSIONS OFF)
target_include_directories(
PROPOSAL PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
diff --git a/Processes/Pythia/testPythia8.cc b/Processes/Pythia/testPythia8.cc
index 6a8ae279badbbb7336357fefb208f926d55b4748..4fd221bd66cdb45378d0f21bd162336d97c82084 100644
--- a/Processes/Pythia/testPythia8.cc
+++ b/Processes/Pythia/testPythia8.cc
@@ -15,9 +15,11 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/geometry/Point.h>
+
#include <corsika/units/PhysicalUnits.h>
#include <corsika/utl/CorsikaFenv.h>
+
#include <catch2/catch.hpp>
TEST_CASE("Pythia", "[processes]") {
@@ -77,10 +79,13 @@ TEST_CASE("Pythia", "[processes]") {
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/UniformMediumType.h>
+#include <corsika/environment/UniformMagneticField.h>
using namespace corsika;
using namespace corsika::units::si;
@@ -97,24 +102,28 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
TEST_CASE("pythia process") {
// setup environment, geometry
- environment::Environment<environment::IMediumModel> env;
-
- geometry::CoordinateSystem const& cs = env.GetCoordinateSystem();
+ setup::Environment env;
+ auto& universe = *(env.GetUniverse());
+ using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
auto theMedium =
- environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
- geometry::Point{cs, 0_m, 0_m, 0_m},
+ setup::Environment::CreateNode<geometry::Sphere>(
+ geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
+ theMedium->SetModelProperties<EnvironmentModel>(
+ environment::EMediumType::eAir,
+ geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Hydrogen},
- std::vector<float>{1.}));
+ std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+
+ auto const* nodePtr = theMedium.get();
+ universe.AddChild(std::move(theMedium));
- auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
+ const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
SECTION("pythia decay") {
feenableexcept(FE_INVALID);
setup::Stack stack;
diff --git a/Processes/QGSJetII/testQGSJetII.cc b/Processes/QGSJetII/testQGSJetII.cc
index eafe725dbcd676afcf402ea562a589d228f468dc..d1e9c7783fd1a57105e01d6d3837b907a36b59ad 100644
--- a/Processes/QGSJetII/testQGSJetII.cc
+++ b/Processes/QGSJetII/testQGSJetII.cc
@@ -111,12 +111,14 @@ TEST_CASE("QgsjetII", "[processes]") {
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/process/qgsjetII/qgsjet-II-04.h>
+#include <corsika/environment/UniformMediumType.h>
+#include <corsika/environment/UniformMagneticField.h>
using namespace corsika::units::si;
using namespace corsika::units;
@@ -124,16 +126,18 @@ using namespace corsika::units;
TEST_CASE("QgsjetIIInterface", "[processes]") {
// setup environment, geometry
- environment::Environment<environment::IMediumModel> env;
+ setup::Environment env;
auto& universe = *(env.GetUniverse());
+ using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
auto theMedium =
- environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ setup::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
+ theMedium->SetModelProperties<EnvironmentModel>(
+ environment::EMediumType::eAir,
+ geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc
index de77d5d33a32b1cb4461abb88552ee14bea415a0..b87d02e5a9fb34d06524b6fa1a522bbac4c015a3 100644
--- a/Processes/Sibyll/NuclearInteraction.cc
+++ b/Processes/Sibyll/NuclearInteraction.cc
@@ -43,7 +43,7 @@ namespace corsika::process::sibyll {
}
template <>
- void NuclearInteraction<SetupEnvironment>::PrintCrossSectionTable(
+ void NuclearInteraction<setup::Environment>::PrintCrossSectionTable(
corsika::particles::Code pCode) {
using namespace corsika::particles;
const int k = targetComponentsIndex_.at(pCode);
@@ -68,7 +68,7 @@ namespace corsika::process::sibyll {
}
template <>
- void NuclearInteraction<SetupEnvironment>::InitializeNuclearCrossSections() {
+ void NuclearInteraction<setup::Environment>::InitializeNuclearCrossSections() {
using namespace corsika::particles;
using namespace units::si;
@@ -134,7 +134,7 @@ namespace corsika::process::sibyll {
}
template <>
- units::si::CrossSectionType NuclearInteraction<SetupEnvironment>::ReadCrossSectionTable(
+ units::si::CrossSectionType NuclearInteraction<setup::Environment>::ReadCrossSectionTable(
const int ia, particles::Code pTarget, units::si::HEPEnergyType elabnuc) {
using namespace corsika::particles;
using namespace units::si;
@@ -154,7 +154,7 @@ namespace corsika::process::sibyll {
template <>
template <>
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
- NuclearInteraction<SetupEnvironment>::GetCrossSection(Particle const& vP,
+ NuclearInteraction<setup::Environment>::GetCrossSection(Particle const& vP,
const particles::Code TargetId) {
using namespace units::si;
if (vP.GetPID() != particles::Code::Nucleus)
@@ -195,7 +195,7 @@ namespace corsika::process::sibyll {
template <>
template <>
- units::si::GrammageType NuclearInteraction<SetupEnvironment>::GetInteractionLength(
+ units::si::GrammageType NuclearInteraction<setup::Environment>::GetInteractionLength(
Particle const& vP) {
using namespace units;
@@ -615,8 +615,8 @@ namespace corsika::process::sibyll {
}
template <>
- NuclearInteraction<SetupEnvironment>::NuclearInteraction(
- process::sibyll::Interaction& hadint, SetupEnvironment const& env)
+ NuclearInteraction<setup::Environment>::NuclearInteraction(
+ process::sibyll::Interaction& hadint, setup::Environment const& env)
: environment_(env)
, hadronicInteraction_(hadint) {
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index dae9741fcb3cba2641d6bd5a9ff885cb3752ed1c..755237efa5d346f133bd0b9f999322e35ba4c5fd 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -75,12 +75,14 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupEnvironment.h>
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
-#include <corsika/process/sibyll/sibyll2.3d.h>
+#include <corsika/environment/UniformMediumType.h>
+#include <corsika/environment/UniformMagneticField.h>
using namespace corsika::units::si;
using namespace corsika::units;
@@ -95,17 +97,19 @@ auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS)
TEST_CASE("SibyllInterface", "[processes]") {
// setup environment, geometry
- environment::Environment<environment::IMediumModel> env;
+ setup::Environment env;
auto& universe = *(env.GetUniverse());
+ using EnvironmentModel = environment::UniformMediumType<environment::UniformMagneticField<environment::HomogeneousMedium<setup::IEnvironment>>>;
auto theMedium =
- environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+ setup::Environment::CreateNode<geometry::Sphere>(
geometry::Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
1_km * std::numeric_limits<double>::infinity());
- using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
- theMedium->SetModelProperties<MyHomogeneousModel>(
- 1_kg / (1_m * 1_m * 1_m),
+ theMedium->SetModelProperties<EnvironmentModel>(
+ environment::EMediumType::eAir,
+ geometry::Vector(env.GetCoordinateSystem(), 0_T, 0_T, 0_T),
+ 1_kg / (1_m * 1_m * 1_m),
environment::NuclearComposition(
std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
diff --git a/Setup/SetupEnvironment.h b/Setup/SetupEnvironment.h
index 94b6ed91f0357fd82f6b217f63f1f6351fce1176..89f5147326dee4c95d038c87590f31e3e0c78e04 100644
--- a/Setup/SetupEnvironment.h
+++ b/Setup/SetupEnvironment.h
@@ -9,10 +9,20 @@
#pragma once
#include <corsika/environment/Environment.h>
+#include <corsika/environment/IMagneticFieldModel.h>
#include <corsika/environment/IMediumModel.h>
-#include <corsika/environment/NameModel.h>
+#include <corsika/environment/IMediumTypeModel.h>
+#include <corsika/environment/IRefractiveIndexModel.h>
+#include <corsika/environment/InhomogeneousMedium.h>
namespace corsika::setup {
- using IEnvironmentModel = environment::IMediumModel;
- using SetupEnvironment = environment::Environment<IEnvironmentModel>;
+
+ /**
+ Definition of the default environemnt model interface. Each model
+ interface provides properties of the environment in a position
+ bdependent way.
+ */
+
+ using IEnvironment = environment::IMediumTypeModel<environment::IMagneticFieldModel<environment::IMediumModel>>;
+ using Environment = environment::Environment<IEnvironment>;
} // namespace corsika::setup