only density parameters as arugments to LayeredSphericalAtmosphereBuilder::addXYZLayer

Documentation/Examples/vertical_EAS.cc

@@ -78,6 +78,13 @@ void registerRandomStreams(const int seed) {
 template <typename T>
 using MEnv = environment::MediumPropertyModel<environment::UniformMagneticField<T>>;
 
+template <typename... TArgs>
+auto make_builder(geometry::Point const& center, TArgs... args) {
+  return environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv,
+                                                        TArgs...>{
+      args..., center, units::constants::EarthRadius::Mean};
+}
+
 int main(int argc, char** argv) {
 
   logging::SetLevel(logging::level::info);
@@ -101,26 +108,17 @@ int main(int argc, char** argv) {
   EnvType env;
   const CoordinateSystem& rootCS = env.GetCoordinateSystem();
   Point const center{rootCS, 0_m, 0_m, 0_m};
-  environment::LayeredSphericalAtmosphereBuilder<setup::EnvironmentInterface, MEnv>
-      builder{center};
+  auto builder = make_builder(center, environment::Medium::AirDry1Atm,
+                              geometry::Vector{rootCS, 0_T, 0_T, 1_T});
   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,
-                              environment::Medium::AirDry1Atm,
-                              geometry::Vector(rootCS, 0_T, 0_T, 1_T));
-  builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km,
-                              environment::Medium::AirDry1Atm,
-                              geometry::Vector(rootCS, 0_T, 0_T, 1_T));
-  builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km,
-                              environment::Medium::AirDry1Atm,
-                              geometry::Vector(rootCS, 0_T, 0_T, 1_T));
-  builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km,
-                              environment::Medium::AirDry1Atm,
-                              geometry::Vector(rootCS, 0_T, 0_T, 1_T));
-  builder.addLinearLayer(1e9_cm, 112.8_km, environment::Medium::AirDry1Atm,
-                         geometry::Vector(rootCS, 0_T, 0_T, 1_T));
+  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

Environment/LayeredSphericalAtmosphereBuilder.h

@@ -18,8 +18,10 @@
 #include <corsika/units/PhysicalConstants.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <functional>
 #include <memory>
 #include <stack>
+#include <tuple>
 #include <type_traits>
 
 namespace corsika::environment {
@@ -52,12 +54,14 @@ namespace corsika::environment {
   } // namespace detail
 
   template <typename TMediumInterface = environment::IMediumModel,
-            template <typename> typename TMediumModelExtra = detail::NoExtraModel>
+            template <typename> typename TMediumModelExtra = detail::NoExtraModel,
+            typename... TModelArgs>
   class LayeredSphericalAtmosphereBuilder {
     std::unique_ptr<NuclearComposition> composition_;
     geometry::Point center_;
     units::si::LengthType previousRadius_{units::si::LengthType::zero()};
     units::si::LengthType earthRadius_;
+    std::tuple<TModelArgs...> const additionalModelArgs_;
 
     std::stack<typename VolumeTreeNode<TMediumInterface>::VTNUPtr>
         layers_; // innermost layer first
@@ -76,18 +80,18 @@ namespace corsika::environment {
 
   public:
     LayeredSphericalAtmosphereBuilder(
-        corsika::geometry::Point center,
+        TModelArgs... args, corsika::geometry::Point center,
         units::si::LengthType earthRadius = units::constants::EarthRadius::Mean)
         : center_(center)
-        , earthRadius_(earthRadius) {}
+        , earthRadius_(earthRadius)
+        , additionalModelArgs_{args...} {}
 
     void setNuclearComposition(NuclearComposition composition) {
       composition_ = std::make_unique<NuclearComposition>(composition);
     }
 
-    template <typename... TArgs>
     void addExponentialLayer(units::si::GrammageType b, units::si::LengthType c,
-                             units::si::LengthType upperBoundary, TArgs&&... args) {
+                             units::si::LengthType upperBoundary) {
       using namespace units::si;
 
       auto const radius = earthRadius_ + upperBoundary;
@@ -100,43 +104,56 @@ namespace corsika::environment {
       auto const rho0 = b / c;
       std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl;
 
-      if constexpr (detail::has_extra_models<TMediumModelExtra>::value)
-        node->template SetModelProperties<
-            TMediumModelExtra<environment::SlidingPlanarExponential<TMediumInterface>>>(
-            args..., center_, rho0, -c, *composition_, earthRadius_);
-      else
-        node->template SetModelProperties<
-            environment::SlidingPlanarExponential<TMediumInterface>>(
+      if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+        // helper lambda in which the last 5 arguments to make_shared<...> are bound
+        auto lastBound = [&](auto... argPack) {
+          return std::make_shared<
+              TMediumModelExtra<environment::SlidingPlanarExponential<TMediumInterface>>>(
+              argPack..., center_, rho0, -c, *composition_, earthRadius_);
+        };
+
+        // now unpack the additional arguments
+        auto model = std::apply(lastBound, additionalModelArgs_);
+        node->SetModelProperties(std::move(model));
+      } else {
+        node->template SetModelProperties<SlidingPlanarExponential<TMediumInterface>>(
             center_, rho0, -c, *composition_, earthRadius_);
+      }
 
       layers_.push(std::move(node));
     }
 
-    template <typename... TArgs>
-    void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary,
-                        TArgs&&... args) {
+    void addLinearLayer(units::si::LengthType c, units::si::LengthType upperBoundary) {
       using namespace units::si;
 
       auto const radius = earthRadius_ + upperBoundary;
       checkRadius(radius);
       previousRadius_ = radius;
 
+      auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
+          std::make_unique<geometry::Sphere>(center_, 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));
+      if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
+        // helper lambda in which the last 2 arguments to make_shared<...> are bound
+        auto lastBound = [&](auto... argPack) {
+          return std::make_shared<
+              TMediumModelExtra<environment::HomogeneousMedium<TMediumInterface>>>(
+              argPack..., rho0, *composition_);
+        };
 
-      if constexpr (detail::has_extra_models<TMediumModelExtra>::value)
-        node->template SetModelProperties<
-            TMediumModelExtra<environment::HomogeneousMedium<TMediumInterface>>>(
-            args..., rho0, *composition_);
-      else
+        // now unpack the additional arguments
+        auto model = std::apply(lastBound, additionalModelArgs_);
+
+        node->SetModelProperties(std::move(model));
+      } else {
         node->template SetModelProperties<
-            environment::HomogeneousMedium<TMediumInterface>>(args..., rho0,
-                                                              *composition_);
+            environment::HomogeneousMedium<TMediumInterface>>(rho0, *composition_);
+      }
 
       layers_.push(std::move(node));
     }