Merge branch 'layered_atmo' into 'master'

Layered atmosphere builder

See merge request !171

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 = {

Environment/CMakeLists.txt

+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>
   )

Environment/LayeredSphericalAtmosphereBuilder.cc

+/*
+ * (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;
+  }
+}

Environment/LayeredSphericalAtmosphereBuilder.h

+/*
+ * (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

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);
     }
   };
 

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);
+}