finished LayeredSphericalAtmosphereBuilder

Environment/CMakeLists.txt

+set (
+  ENVIRONMENT_SOURCES
+  Convenience.cpp
+)
+
 set (
   ENVIRONMENT_HEADERS
   VolumeTreeNode.h
@@ -13,6 +18,7 @@ set (
   BaseExponential.h
   FlatExponential.h
   SlidingPlanarExponential.h
+  Convenience.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/Convenience.cpp

+/*
+ * (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/Convenience.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 * exp(seaLevel_ / c);
+
+  node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>(center_, rho0, -c,
+                                                                   *composition_);
+
+  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;
+
+  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> environment_;
+
+  auto& universe = environment_.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;
+  }
+
+  return environment_;
+}

Environment/Convenience.h

 /*
- * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
  *
  * See file AUTHORS for a list of contributors.
  *
@@ -8,42 +8,47 @@
  * the license.
  */
 
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/IMediumModel.h>
 #include <corsika/environment/NuclearComposition.h>
-#include <corsika/environment/FlatExponential.h>
-#include <corsika/environment/SlidingPlanarExponential.h>
 #include <corsika/environment/VolumeTreeNode.h>
-#include <corsika/particles/ParticleProperties.h>
+#include <corsika/geometry/Point.h>
 #include <corsika/units/PhysicalUnits.h>
 
 #include <memory>
+#include <stack>
 
 namespace corsika::environment {
 
-class LayeredAtmosphereBuilder {
-    std::unique_ptr<NuclearComposition> nuclearComposition_;
+  class LayeredSphericalAtmosphereBuilder {
+    std::unique_ptr<NuclearComposition> composition_;
     geometry::Point center_;
-    units::si::LengthType previousRadius_{LengthType::zero()};
-    static auto constexpr earthRadius_ = 6371_km;
-    
-public:
-    LayeredAtmosphereBuilder(corsika::geometry::Point center) : center_(center) {}
-    
-    void setNuclearComposition(NuclearComposition composition) {
-        nuclearComposition_ = std::make_unique<NuclearComposition>(composition);
-    }
-    
-    void addLayer(units::si::GrammageType a, units::si::GrammageType b, units::si::LengthType c, units::si::LengthType thickness) {
-        auto const radius = previousRadius_ + thickness;
-        
-        std::make_unique<BaseNodeType>(
-          std::make_unique<geometry::Sphere>(center_, radius);
-          
-        previousRadius_ = radiusl
-    }
-    
-    auto const assemble() {
-        // get the outmost VolumeTreeNode
-    }
-};
-
-}
+    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);
+
+    Environment<IMediumModel> assemble();
+  };
+
+} // namespace corsika::environment

Environment/testEnvironment.cc

@@ -8,6 +8,7 @@
  * the license.
  */
 
+#include <corsika/environment/Convenience.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);
+}