diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 01d745f4530738ae7d60e9edc858dffd8761a58a..8fd61fcde6eec20269e0a489a1ef2e6adc588137 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/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
diff --git a/Environment/LayeredSphericalAtmosphereBuilder.h b/Environment/LayeredSphericalAtmosphereBuilder.h
index 9c81a5a6ff71edb6b8814b2f08e055ab6c4b5c33..2a297005f8a726865f19a16531eefbfbdb991144 100644
--- a/Environment/LayeredSphericalAtmosphereBuilder.h
+++ b/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));
}