us -> 5layer

corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl

@@ -40,8 +40,9 @@ namespace corsika {
   inline typename LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
                                                     TModelArgs...>::volume_tree_node*
   LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
-      addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary) {
+      addExponentialLayer(GrammageType const b, LengthType const scaleHeight, LengthType const upperBoundary) {
 
+    // outer radius
     auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
@@ -49,14 +50,14 @@ namespace corsika {
     auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
         std::make_unique<Sphere>(center_, radius));
 
-    auto const rho0 = b / c;
+    auto const rho0 = b / scaleHeight;
 
     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<SlidingPlanarExponential<TMediumInterface>>>(
-            argPack..., center_, rho0, -c, *composition_, planetRadius_);
+            argPack..., center_, rho0, -scaleHeight, *composition_, planetRadius_);
       };
 
       // now unpack the additional arguments
@@ -64,7 +65,7 @@ namespace corsika {
       node->setModelProperties(std::move(model));
     } else {
       node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
-          center_, rho0, -c, *composition_, planetRadius_);
+          center_, rho0, -scaleHeight, *composition_, planetRadius_);
     }
 
     layers_.push(std::move(node));
@@ -75,7 +76,9 @@ namespace corsika {
             typename... TModelArgs>
   inline void LayeredSphericalAtmosphereBuilder<
       TMediumInterface, TMediumModelExtra,
-      TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
+      TModelArgs...>::addLinearLayer(LengthType const scaleHeight,
+                                     LengthType const upperBoundary) {
+    // outer radius
     auto const radius = planetRadius_ + upperBoundary;
     checkRadius(radius);
     previousRadius_ = radius;
@@ -83,8 +86,8 @@ namespace corsika {
     auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
         std::make_unique<Sphere>(center_, radius));
 
-    units::si::GrammageType constexpr b = 1_g / (1_cm * 1_cm);
-    auto const rho0 = b / c;
+    GrammageType constexpr b = 1_g / (1_cm * 1_cm);
+    auto const rho0 = b / scaleHeight;
 
     if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
       // helper lambda in which the last 2 arguments to make_shared<...> are bound

corsika/framework/utility/ImplementsMixin.hpp

@@ -36,5 +36,8 @@ namespace corsika {
       static constexpr bool value = type::value;
     };
 
+    template <template <typename> typename Mixin, typename T>
+    typedef implements_mixin<Mixin, T>::value implements_mixin_v;
+
   } // namespace detail
 } // namespace corsika

corsika/media/LayeredSphericalAtmosphereBuilder.hpp

@@ -35,7 +35,7 @@ namespace corsika {
    * function `create` does then take an unspecified number of extra
    * parameters to internalize those models for all layers later
    * produced.
-   **/
+   */
   template <typename TMediumInterface = IMediumModel,
             template <typename> typename MExtraEnvirnoment = detail::NoExtraModel>
   struct make_layered_spherical_atmosphere_builder;
@@ -49,7 +49,6 @@ namespace corsika {
    *
    * Each layer by definition has a density profile and a (constant)
    * nuclear composition model.
-   *
    */
 
   template <typename TMediumInterface = IMediumModel,
@@ -79,9 +78,9 @@ namespace corsika {
     typedef typename VolumeTreeNode<TMediumInterface>::VTNUPtr volume_tree_node_uptr;
 
     void setNuclearComposition(NuclearComposition const& composition);
-    volume_tree_node* addExponentialLayer(GrammageType b, LengthType c,
-                                          LengthType upperBoundary);
-    void addLinearLayer(LengthType c, LengthType upperBoundary);
+    volume_tree_node* addExponentialLayer(GrammageType const b, LengthType const scaleHeight,
+                                          LengthType const upperBoundary);
+    void addLinearLayer(LengthType const c, LengthType const upperBoundary);
 
     void addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
                          unsigned int const nBins, LengthType const deltaHeight,

corsika/media/USStandardAtmosphere.hpp

@@ -14,32 +14,86 @@
 
 namespace corsika {
 
+  /**
+   * Atmosphere Ids following the CORSIKA 7 5-layered atmosphere models.
+   *
+   * Each model corresponds to a standard 5-layered atmosphere model. 4 Layers are
+   * exponential, the outer layer is with constant density.
+   *
+   * All atmospheres are valid for heights (above Earth sea level) up to 112.8 km.
+   */
+  enum class AtmosphereId : uint8_t {
+    LinsleyUSStd = 0,
+    /*    MiddleEuropeJan,
+        MiddleEuropeFeb,
+        MiddleEuropeMay,
+        MiddleEuropeJun,
+        MiddleEuropeAug,
+        MiddleEuropeOct,
+        MiddleEuropeDec,
+        SouthPoleMar,
+        SouthPoleJul,
+        SouthPoleOct,
+        SouthPoleDec,
+        SouthPoleJan,
+        SouthPoleAug,
+        MalargueWinterI,
+        MalargueWinterII,
+        MalargueSpring,
+        MalargueSummer,
+        MalargueAutumn,
+        USStdBK,*/
+    LastAtmosphere
+  };
+
+  struct AtmosphereLayerParameters {
+    LengthType altitude;
+    GrammageType offset;
+    LengthType scaleHeight;
+  };
+
+  typedef std::array<AtmosphereLayerParameters, 5> AtmosphereParameters;
+
+  std::arrary<AtmosphereParameters,
+              static_cast<uint8_t>(
+                  AtmosphereId::LastAtmosphere)> constexpr AtmosphereParameterList{
+      {{4_km, 1222.6562_g / cube(1_cm), 994186.38_cm},
+       {10_km, 1144.9069_g / cube(1_cm), 878153.55_cm},
+       {40_km, 1305.5948_g / cube(1_cm), 636143.04_cm},
+       {100_km, 540.1778_g / cube(1_cm), 772170.16_cm, 100_km},
+       {112.8_km, 1_g / cube(1_cm), 1e9_cm}}};
+
   template <typename TEnvironmentInterface, template <typename> typename TExtraEnv,
             typename TEnvironment, typename... TArgs>
-  auto create_us_standard_atmosphere(TEnvironment& env, Point const& center,
-                                     LengthType const& earthRadius, TArgs... args) {
+  auto create_5layer_atmosphere(AtmosphereId const atmId, TEnvironment& env,
+                                Point const& center, TArgs... args) {
 
     // construct the atmosphere builder
     auto builder = make_layered_spherical_atmosphere_builder<
-        TEnvironmentInterface, TExtraEnv>::create(center, earthRadius,
+        TEnvironmentInterface, TExtraEnv>::create(center, constants::EarthRadius::Mean,
                                                   std::forward<TArgs>(args)...);
 
     // as per the vertical_EAS reference, we do not include Ar for now.
     // TODO: This is not a US standard atmosphere
     builder.setNuclearComposition(
-        {{Code::Nitrogen, Code::Oxygen}, {0.7847f, 1.f - 0.7847f}});
+        {{Code::Nitrogen, Code::Oxygen}, {0.7847, 1. - 0.7847}});
 
     // add the standard atmosphere layers
-    builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
-    builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
+    auto params = AtmosphereParameterList[static_cast<uint8_t>(atmId)];
+    for (int i = 0; i < 4; ++i) {
+      builder.addExponentialLayer(params[i][1], params[i][2], params[i][0]);
+    }
+    builder.addLinearLayer(params[4][2], params[4][0]);
+    /*
+      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 + constants::EarthRadius::Mean);
+    builder.addLinearLayer(1e9_cm, 112.8_km);*/
 
     // check if we want to also add the US standard refractivity
-    if constexpr (detail::implements_mixin<IRefractiveIndexModel,
-                                           TEnvironmentInterface>::value) {
+    if constexpr (detail::implements_mixin_v<IRefractiveIndexModel,
+                                             TEnvironmentInterface>) {
 
       // TODO: Add US Standard refractivity
     }