diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index efa57f50f27d68a871b4af6c5c74e29ae5c73709..6ccd877f6b55e45ed1f50d3eb1d94f26ef54281e 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/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
diff --git a/corsika/framework/utility/ImplementsMixin.hpp b/corsika/framework/utility/ImplementsMixin.hpp
index 45eec946ed475b755a6670d989664939cbaee7e1..cd6857b96d07a5b46dddd701fea9066d60c56ee6 100644
--- a/corsika/framework/utility/ImplementsMixin.hpp
+++ b/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
diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
index 09bcabc06472717a79ed14f62778843aac042f4f..a5c1ab418953637931b337f8e9ea8ad833818aed 100644
--- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+++ b/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,
diff --git a/corsika/media/USStandardAtmosphere.hpp b/corsika/media/USStandardAtmosphere.hpp
index 7194d79f187f602a5dfcbd750e8f29b7cd9fea61..5b228598f80f6fe462739e80bd8c84b95e0d3785 100644
--- a/corsika/media/USStandardAtmosphere.hpp
+++ b/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
}