diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 78972df91f591a56a23c0bf6aeb6e8022a293549..ca27eac3d43e6f37cac296af3e7523de6ac7fa6b 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -20,6 +20,7 @@
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/Environment.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/NuclearComposition.h>
@@ -79,26 +80,19 @@ int main() {
using EnvType = Environment<setup::IEnvironmentModel>;
EnvType env;
const CoordinateSystem& rootCS = env.GetCoordinateSystem();
- auto& universe = *(env.GetUniverse());
-
- auto theMedium = EnvType::CreateNode<Sphere>(
- Point{rootCS, 0_m, 0_m, 0_m}, 1_km * std::numeric_limits<double>::infinity());
-
- auto constexpr temperature = 295_K; // AIRES default temperature for isothermal model
- auto constexpr lambda =
- -constants::R * temperature / (constants::g_sub_n * 28.966_g / mole);
- // 1036 g/cm² taken from AIRES code
- auto constexpr rho0 = -1036_g / square(1_cm) / lambda;
- using FlatExp = environment::FlatExponential<environment::IMediumModel>;
- theMedium->SetModelProperties<FlatExp>(
- Point{rootCS, 0_m, 0_m, 0_m}, Vector<dimensionless_d>{rootCS, {0., 0., 1.}}, rho0,
- lambda,
- environment::NuclearComposition(
- std::vector<particles::Code>{particles::Code::Nitrogen,
- particles::Code::Oxygen},
- std::vector<float>{
- 0.7847f,
- 1.f - 0.7847f})); // values taken from AIRES manual, Ar removed for now
+
+ environment::LayeredSphericalAtmosphereBuilder builder(Point{rootCS, 0_m, 0_m, 0_m});
+ builder.setNuclearComposition(
+ {{particles::Code::Nitrogen, particles::Code::Oxygen},
+ {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+
+ 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);
+
+ builder.assemble(env);
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -110,7 +104,9 @@ int main() {
double phi = 0.;
Point const injectionPos(
- rootCS, 0_m, 0_m, 112.8_km); // this is the CORSIKA 7 start of atmosphere/universe
+ rootCS, 0_m, 0_m,
+ 112.8_km * 0.999 +
+ builder.earthRadius); // this is the CORSIKA 7 start of atmosphere/universe
// {
auto elab2plab = [](HEPEnergyType Elab, HEPMassType m) {
@@ -137,17 +133,10 @@ int main() {
Line const line(injectionPos, plab.normalized() * 1_m * 1_Hz);
auto const velocity = line.GetV0().norm();
- auto const observationHeight = 1.425_km;
+ auto const observationHeight = 1.425_km + builder.earthRadius;
setup::Trajectory const showerAxis(line, (112.8_km - observationHeight) / velocity);
- auto const grammage = theMedium->GetModelProperties().IntegratedGrammage(
- showerAxis, (112.8_km - observationHeight));
- std::cout << "Grammage to ground: " << grammage / (1_g / square(1_cm)) << " g/cm²"
- << std::endl;
-
- universe.AddChild(std::move(theMedium));
-
// setup processes, decays and interactions
const std::vector<particles::Code> trackedHadrons = {
diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index 28e5523163014cef5709d39d8d113e738d282889..f4b93c759a2f32128d98406b8e4130dd851b8c45 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -1,3 +1,8 @@
+set (
+ ENVIRONMENT_SOURCES
+ LayeredSphericalAtmosphereBuilder.cc
+)
+
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
@@ -13,6 +18,7 @@ set (
BaseExponential.h
FlatExponential.h
SlidingPlanarExponential.h
+ LayeredSphericalAtmosphereBuilder.h
)
set (
@@ -20,13 +26,20 @@ set (
corsika/environment
)
-add_library (CORSIKAenvironment INTERFACE)
+add_library (CORSIKAenvironment STATIC ${ENVIRONMENT_SOURCES})
CORSIKA_COPY_HEADERS_TO_NAMESPACE (CORSIKAenvironment ${ENVIRONMENT_NAMESPACE} ${ENVIRONMENT_HEADERS})
+set_target_properties (
+ CORSIKAgeometry
+ PROPERTIES
+ VERSION ${PROJECT_VERSION}
+ SOVERSION 1
+ PUBLIC_HEADER "${ENVIRONMENT_HEADERS}"
+ )
+
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
CORSIKAenvironment
- INTERFACE
CORSIKAgeometry
CORSIKAparticles
CORSIKAunits
@@ -35,7 +48,7 @@ target_link_libraries (
target_include_directories (
CORSIKAenvironment
- INTERFACE
+ PUBLIC
$<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
$<INSTALL_INTERFACE:include>
)
diff --git a/Environment/LayeredSphericalAtmosphereBuilder.cc b/Environment/LayeredSphericalAtmosphereBuilder.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4d1631d54e42431cd1b36299af9083c643bce0e7
--- /dev/null
+++ b/Environment/LayeredSphericalAtmosphereBuilder.cc
@@ -0,0 +1,85 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/LayeredSphericalAtmosphereBuilder.h>
+#include <corsika/environment/FlatExponential.h>
+#include <corsika/environment/HomogeneousMedium.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 = seaLevel_ + 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_, seaLevel_);
+
+ layers_.push(std::move(node));
+}
+
+void LayeredSphericalAtmosphereBuilder::addLinearLayer(
+ units::si::LengthType c, units::si::LengthType upperBoundary) {
+ using namespace units::si;
+
+ auto const radius = seaLevel_ + 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
new file mode 100644
index 0000000000000000000000000000000000000000..ab2b8a25d95887d198dc568536fec81c3d40b461
--- /dev/null
+++ b/Environment/LayeredSphericalAtmosphereBuilder.h
@@ -0,0 +1,55 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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/Environment.h>
+#include <corsika/environment/IMediumModel.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/VolumeTreeNode.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <memory>
+#include <stack>
+
+namespace corsika::environment {
+
+ class LayeredSphericalAtmosphereBuilder {
+ std::unique_ptr<NuclearComposition> composition_;
+ geometry::Point center_;
+ units::si::LengthType previousRadius_{units::si::LengthType::zero()};
+ units::si::LengthType seaLevel_;
+
+ std::stack<VolumeTreeNode<environment::IMediumModel>::VTNUPtr>
+ layers_; // innermost layer first
+
+ void checkRadius(units::si::LengthType) const;
+
+ public:
+ static auto constexpr earthRadius = 6'371'000 * units::si::meter;
+
+ LayeredSphericalAtmosphereBuilder(corsika::geometry::Point center,
+ units::si::LengthType seaLevel = earthRadius)
+ : center_(center)
+ , seaLevel_(seaLevel) {}
+
+ void setNuclearComposition(NuclearComposition);
+
+ void addExponentialLayer(units::si::GrammageType, units::si::LengthType,
+ units::si::LengthType);
+
+ auto size() const { return layers_.size(); }
+
+ void addLinearLayer(units::si::LengthType, units::si::LengthType);
+
+ void assemble(Environment<IMediumModel>&);
+ Environment<IMediumModel> assemble();
+ };
+
+} // namespace corsika::environment
diff --git a/Environment/SlidingPlanarExponential.h b/Environment/SlidingPlanarExponential.h
index 71f8e27d917e3e34fc5d8c2fc663d5142d07aa59..fdbde0f4e1d3257ba73768ea1e4d123e6ee6fbbc 100644
--- a/Environment/SlidingPlanarExponential.h
+++ b/Environment/SlidingPlanarExponential.h
@@ -38,36 +38,38 @@ namespace corsika::environment {
template <class T>
class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
public T {
- NuclearComposition const fNuclComp;
+ NuclearComposition const nuclComp_;
+ units::si::LengthType const referenceHeight_;
using Base = BaseExponential<SlidingPlanarExponential<T>>;
public:
- SlidingPlanarExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
- units::si::LengthType vLambda, NuclearComposition vNuclComp)
- : Base(vP0, vRho, vLambda)
- , fNuclComp(vNuclComp) {}
+ SlidingPlanarExponential(geometry::Point const& p0, units::si::MassDensityType rho0,
+ units::si::LengthType lambda, NuclearComposition nuclComp, units::si::LengthType referenceHeight = units::si::LengthType::zero())
+ : Base(p0, rho0, lambda)
+ , nuclComp_(nuclComp),
+ referenceHeight_(referenceHeight) {}
units::si::MassDensityType GetMassDensity(
- geometry::Point const& vP) const override {
- auto const height = (vP - Base::fP0).norm();
+ geometry::Point const& p) const override {
+ auto const height = (p - Base::fP0).norm() - referenceHeight_;
return Base::fRho0 * exp(Base::fInvLambda * height);
}
- NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+ NuclearComposition const& GetNuclearComposition() const override { return nuclComp_; }
units::si::GrammageType IntegratedGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
- units::si::LengthType vL) const override {
- auto const axis = (vLine.GetR0() - Base::fP0).normalized();
- return Base::IntegratedGrammage(vLine, vL, axis);
+ geometry::Trajectory<geometry::Line> const& line,
+ units::si::LengthType l) const override {
+ auto const axis = (line.GetR0() - Base::fP0).normalized();
+ return Base::IntegratedGrammage(line, l, axis);
}
units::si::LengthType ArclengthFromGrammage(
- geometry::Trajectory<geometry::Line> const& vLine,
- units::si::GrammageType vGrammage) const override {
- auto const axis = (vLine.GetR0() - Base::fP0).normalized();
- return Base::ArclengthFromGrammage(vLine, vGrammage, axis);
+ geometry::Trajectory<geometry::Line> const& line,
+ units::si::GrammageType grammage) const override {
+ auto const axis = (line.GetR0() - Base::fP0).normalized();
+ return Base::ArclengthFromGrammage(line, grammage, axis);
}
};
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index cae91a86000b027c3e2687d9515df977c23d0eb9..5c314d1d07b5b9a5faad15ea2345114a5fb8f3f1 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -8,6 +8,7 @@
* the license.
*/
+#include <corsika/environment/LayeredSphericalAtmosphereBuilder.h>
#include <corsika/environment/DensityFunction.h>
#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
@@ -196,3 +197,35 @@ TEST_CASE("InhomogeneousMedium") {
inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
}
}
+
+TEST_CASE("LayeredSphericalAtmosphereBuilder") {
+ LayeredSphericalAtmosphereBuilder builder(gOrigin);
+ builder.setNuclearComposition(
+ {{{particles::Code::Nitrogen, particles::Code::Oxygen}}, {{.6, .4}}});
+
+ builder.addLinearLayer(1_km, 10_km);
+ builder.addLinearLayer(2_km, 20_km);
+ builder.addLinearLayer(3_km, 30_km);
+
+ REQUIRE(builder.size() == 3);
+
+ auto const builtEnv = builder.assemble();
+ auto const& univ = builtEnv.GetUniverse();
+
+ REQUIRE(builder.size() == 0);
+
+ auto constexpr R = LayeredSphericalAtmosphereBuilder::earthRadius;
+
+ 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);
+}