diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 9edd10b54c3a6daa25a00920c58a0fe6a6c2ad5f..3a420038a2784f6e16ebe979c3f5c3fa16877d46 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -1,6 +1,6 @@
set (
ENVIRONMENT_SOURCES
- LayeredSphericalAtmosphereBuilder.cc
+ # LayeredSphericalAtmosphereBuilder.cc
ShowerAxis.cc
)
diff --git a/Environment/LayeredSphericalAtmosphereBuilder.cc b/Environment/LayeredSphericalAtmosphereBuilder.cc
deleted file mode 100644
index 4f5b33857d19718e1b9502b11751e90aec5d1f7c..0000000000000000000000000000000000000000
--- a/Environment/LayeredSphericalAtmosphereBuilder.cc
+++ /dev/null
@@ -1,83 +0,0 @@
-/*
- * (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.
- */
-
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/HomogeneousMedium.h>
-#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
-#include <corsika/environment/SlidingPlanarExponential.h>
-
-using namespace corsika::environment;
-
-void LayeredSphericalAtmosphereBuilder::checkRadius(units::si::LengthType r) const {
- if (r <= previousRadius_) {
- throw std::runtime_error("radius must be greater than previous");
- }
-}
-
-void LayeredSphericalAtmosphereBuilder::setNuclearComposition(
- NuclearComposition composition) {
- composition_ = std::make_unique<NuclearComposition>(composition);
-}
-
-void LayeredSphericalAtmosphereBuilder::addExponentialLayer(
- units::si::GrammageType b, units::si::LengthType c,
- units::si::LengthType upperBoundary) {
- auto const radius = earthRadius_ + upperBoundary;
- checkRadius(radius);
- previousRadius_ = radius;
-
- auto node = std::make_unique<VolumeTreeNode<IMediumModel>>(
- std::make_unique<geometry::Sphere>(center_, radius));
-
- auto const rho0 = b / c;
- std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
-
- node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(
- center_, rho0, -c, *composition_, earthRadius_);
-
- layers_.push(std::move(node));
-}
-
-void LayeredSphericalAtmosphereBuilder::addLinearLayer(
- units::si::LengthType c, units::si::LengthType upperBoundary) {
- using namespace units::si;
-
- auto const radius = earthRadius_ + upperBoundary;
- checkRadius(radius);
- previousRadius_ = radius;
-
- units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
- auto const rho0 = b / c;
-
- std::cout << "rho0 = " << rho0;
-
- auto node = std::make_unique<VolumeTreeNode<environment::IMediumModel>>(
- std::make_unique<geometry::Sphere>(center_, radius));
- node->SetModelProperties<HomogeneousMedium<IMediumModel>>(rho0, *composition_);
-
- layers_.push(std::move(node));
-}
-
-Environment<IMediumModel> LayeredSphericalAtmosphereBuilder::assemble() {
- Environment<IMediumModel> env;
- assemble(env);
- return env;
-}
-
-void LayeredSphericalAtmosphereBuilder::assemble(Environment<IMediumModel>& env) {
- auto& universe = env.GetUniverse();
- auto* outmost = universe.get();
-
- while (!layers_.empty()) {
- auto l = std::move(layers_.top());
- auto* tmp = l.get();
- outmost->AddChild(std::move(l));
- layers_.pop();
- outmost = tmp;
- }
-}
diff --git a/Environment/LayeredSphericalAtmosphereBuilder.h b/Environment/LayeredSphericalAtmosphereBuilder.h
index ad945cb8461a378145e1d3606461f30711dc69cf..3c383c7b7b4107ff3a3639d4832ae4a478f5b672 100644
--- a/Environment/LayeredSphericalAtmosphereBuilder.h
+++ b/Environment/LayeredSphericalAtmosphereBuilder.h
@@ -9,6 +9,7 @@
#include <corsika/environment/Environment.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/SlidingPlanarExponential.h>
#include <corsika/environment/VolumeTreeNode.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalConstants.h>
@@ -19,16 +20,33 @@
namespace corsika::environment {
+ /**
+ * Helper class to setup concentric spheres of layered atmosphere
+ * with spcified density profiles (exponential, linear, ...).
+ *
+ * This can be used most importantly to replicate CORSIKA7
+ * atmospheres.
+ *
+ * Each layer by definition has a density profile and a (constant)
+ * nuclear composition model.
+ *
+ */
+
+ template <typename TMediumInterface = environment::IMediumModel>
class LayeredSphericalAtmosphereBuilder {
std::unique_ptr<NuclearComposition> composition_;
geometry::Point center_;
units::si::LengthType previousRadius_{units::si::LengthType::zero()};
units::si::LengthType earthRadius_;
- std::stack<VolumeTreeNode<environment::IMediumModel>::VTNUPtr>
+ std::stack<typename VolumeTreeNode<TMediumInterface>::VTNUPtr>
layers_; // innermost layer first
- void checkRadius(units::si::LengthType) const;
+ void checkRadius(units::si::LengthType r) const {
+ if (r <= previousRadius_) {
+ throw std::runtime_error("radius must be greater than previous");
+ }
+ }
public:
LayeredSphericalAtmosphereBuilder(
@@ -37,22 +55,119 @@ namespace corsika::environment {
: center_(center)
, earthRadius_(earthRadius) {}
- void setNuclearComposition(NuclearComposition);
+ void setNuclearComposition(NuclearComposition composition) {
+ composition_ = std::make_unique<NuclearComposition>(composition);
+ }
+
+ void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
+ units::si::LengthType upperBoundary) {
+ auto const radius = earthRadius_ + upperBoundary;
+ checkRadius(radius);
+ previousRadius_ = radius;
+
+ auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+ std::make_unique<geometry::Sphere>(center_, radius));
+
+ auto const rho0 = b / c;
+ std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
+
+ node->template SetModelProperties<
+ environment::SlidingPlanarExponential<TMediumInterface>>(
+ center_, rho0, -c, *composition_, earthRadius_);
+
+ layers_.push(std::move(node));
+ }
+
+ template <template <typename> typename MExtraModels, typename... TArgs>
+ void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
+ units::si::LengthType upperBoundary, TArgs&&... args) {
+ auto const radius = earthRadius_ + upperBoundary;
+ checkRadius(radius);
+ previousRadius_ = radius;
+
+ auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+ std::make_unique<geometry::Sphere>(center_, radius));
+
+ auto const rho0 = b / c;
+ std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
+
+ node->template SetModelProperties<
+ MExtraModels<SlidingPlanarExponential<TMediumInterface>>>(
+ args..., center_, rho0, -c, *composition_, earthRadius_);
+
+ layers_.push(std::move(node));
+ }
+
+ int size() const { return layers_.size(); }
+
+ template <template <typename> typename MExtraModels, typename... TArgs>
+ void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary,
+ TArgs&&... args) {
+ using namespace units::si;
- void addExponentialLayer(units::si::GrammageType, units::si::LengthType,
- units::si::LengthType);
+ auto const radius = earthRadius_ + upperBoundary;
+ checkRadius(radius);
+ previousRadius_ = radius;
- auto size() const { return layers_.size(); }
+ units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
+ auto const rho0 = b / c;
- void addLinearLayer(units::si::LengthType, units::si::LengthType);
+ std::cout << "rho0 = " << rho0;
- void assemble(Environment<IMediumModel>&);
- Environment<IMediumModel> assemble();
+ auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+ std::make_unique<geometry::Sphere>(center_, radius));
+
+ node->template SetModelProperties<
+ MExtraModels<HomogeneousMedium<TMediumInterface>>>(args..., rho0,
+ *composition_);
+
+ layers_.push(std::move(node));
+ }
+
+ void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary) {
+ using namespace units::si;
+
+ auto const radius = earthRadius_ + upperBoundary;
+ checkRadius(radius);
+ previousRadius_ = radius;
+
+ units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
+ auto const rho0 = b / c;
+
+ std::cout << "rho0 = " << rho0;
+
+ auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+ std::make_unique<geometry::Sphere>(center_, radius));
+ node->template SetModelProperties<HomogeneousMedium<TMediumInterface>>(
+ rho0, *composition_);
+
+ layers_.push(std::move(node));
+ }
+
+ void assemble(Environment<TMediumInterface>& env) {
+ auto& universe = env.GetUniverse();
+ auto* outmost = universe.get();
+
+ while (!layers_.empty()) {
+ auto l = std::move(layers_.top());
+ auto* tmp = l.get();
+ outmost->AddChild(std::move(l));
+ layers_.pop();
+ outmost = tmp;
+ }
+ }
+
+ Environment<TMediumInterface> assemble() {
+ Environment<TMediumInterface> env;
+ assemble(env);
+ return env;
+ }
/**
* Get the current Earth radius.
*/
- units::si::LengthType getEarthRadius() const { return earthRadius_; };
- };
+ units::si::LengthType getEarthRadius() const { return earthRadius_; }
+
+ }; // end class LayeredSphericalAtmosphereBuilder
} // namespace corsika::environment
diff --git a/Environment/MediumTypes.h b/Environment/MediumTypes.h
index d07e68aeccc979077ca241ce3167cdd8e6df3eaa..b2f7a0aaa4a5fce02d7ad3150a1f04828a42e1ba 100644
--- a/Environment/MediumTypes.h
+++ b/Environment/MediumTypes.h
@@ -10,6 +10,13 @@
namespace corsika::environment {
+ /**
+ * Medium types are useful most importantly for effective models
+ * like energy losses. a particular medium (mixture of components)
+ * may have specif properties not reflected by its mixture of
+ * components.
+ */
+
enum class EMediumType { eUnknown, eAir, eWater, eIce, eRock };
}
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 523eac4ecfc8861f33a4e0e8d3361b2224085d52..3acb988e85dbf2f097cf2c3517202af3cf65166f 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -107,7 +107,7 @@ namespace corsika::environment {
template <typename ModelProperties, typename... Args>
auto SetModelProperties(Args&&... args) {
static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
- "unusable type provided");
+ "unusable model properties type provided");
fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
return fModelProperties;
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 7394abf42b66e03f6489e8ca9e3885a599c92b65..99aec33f04e4d7aa0f0489c66397680217674907 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -208,7 +208,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
builder.addLinearLayer(1_km, 10_km);
builder.addLinearLayer(2_km, 20_km);
- builder.addLinearLayer(3_km, 30_km);
+ builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 30_km);
CHECK(builder.size() == 3);
@@ -295,6 +295,50 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
CHECK((medium.ArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
+TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
+
+ // setup our interface types
+ using IModelInterface = IMagneticFieldModel<IMediumModel>;
+
+ // the composition we use for the homogenous medium
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ // create magnetic field vectors
+ Vector B0(gCS, 0_T, 0_T, 1_T);
+ Vector B1(gCS, 1_T, 1_T, 0_T);
+
+ // LayeredSphericalAtmosphereBuilder<IMediumModel> builder(gOrigin);
+ LayeredSphericalAtmosphereBuilder<IModelInterface> builder(gOrigin);
+ builder.setNuclearComposition(
+ {{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
+
+ builder.addLinearLayer<UniformMagneticField>(1_km, 10_km, B0);
+ builder.addExponentialLayer<UniformMagneticField>(1222.6562_g / (1_cm * 1_cm),
+ 994186.38_cm, 20_km, B1);
+
+ CHECK(builder.size() == 2);
+
+ auto const builtEnv = builder.assemble();
+ auto const& univ = builtEnv.GetUniverse();
+
+ CHECK(builder.size() == 0);
+ CHECK(univ->GetChildNodes().size() == 1);
+ auto const R = builder.getEarthRadius();
+
+ // check magnetic field at several locations
+ const Point pTest(gCS, -10_m, 4_m, R+35_m);
+ CHECK(B0.GetComponents(gCS) == univ->GetContainingNode(pTest)
+ ->GetModelProperties()
+ .GetMagneticField(pTest)
+ .GetComponents(gCS));
+ const Point pTest2(gCS, 10_m, -4_m, R+15_km);
+ CHECK(B1.GetComponents(gCS) == univ->GetContainingNode(pTest2)
+ ->GetModelProperties()
+ .GetMagneticField(pTest2)
+ .GetComponents(gCS));
+}
+
TEST_CASE("UniformRefractiveIndex w/ Homogeneous") {
// setup our interface types