From 7ac853657a178d0fc2455ea631000f3c468f73aa Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Thu, 21 Oct 2021 20:09:26 +0200
Subject: [PATCH] added all c7 atmospheres
---
.../LayeredSphericalAtmosphereBuilder.inl | 12 +-
corsika/framework/utility/ImplementsMixin.hpp | 2 +-
.../LayeredSphericalAtmosphereBuilder.hpp | 10 +-
corsika/media/USStandardAtmosphere.hpp | 181 +++++++++++++-----
examples/vertical_EAS.cpp | 22 +--
tests/media/testEnvironment.cpp | 10 +-
tests/modules/testCONEX.cpp | 2 +-
7 files changed, 161 insertions(+), 78 deletions(-)
diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index 6ccd877f6..3659e8873 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
@@ -1,4 +1,4 @@
-/*
+w/*
* (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* This software is distributed under the terms of the GNU General Public
@@ -19,8 +19,8 @@ namespace corsika {
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
- inline void LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
- TModelArgs...>::checkRadius(LengthType r)
+ inline void LayeredSphericalAtmosphereBuilder<
+ TMediumInterface, TMediumModelExtra, TModelArgs...>::checkRadius(LengthType const r)
const {
if (r <= previousRadius_) {
throw std::runtime_error("radius must be greater than previous");
@@ -40,7 +40,8 @@ namespace corsika {
inline typename LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra,
TModelArgs...>::volume_tree_node*
LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
- addExponentialLayer(GrammageType const b, LengthType const scaleHeight, LengthType const upperBoundary) {
+ addExponentialLayer(GrammageType const b, LengthType const scaleHeight,
+ LengthType const upperBoundary) {
// outer radius
auto const radius = planetRadius_ + upperBoundary;
@@ -76,7 +77,7 @@ namespace corsika {
typename... TModelArgs>
inline void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
- TModelArgs...>::addLinearLayer(LengthType const scaleHeight,
+ TModelArgs...>::addLinearLayer(GrammageType const b, LengthType const scaleHeight,
LengthType const upperBoundary) {
// outer radius
auto const radius = planetRadius_ + upperBoundary;
@@ -86,7 +87,6 @@ namespace corsika {
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<Sphere>(center_, radius));
- GrammageType constexpr b = 1_g / (1_cm * 1_cm);
auto const rho0 = b / scaleHeight;
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
diff --git a/corsika/framework/utility/ImplementsMixin.hpp b/corsika/framework/utility/ImplementsMixin.hpp
index cd6857b96..9082a9bcd 100644
--- a/corsika/framework/utility/ImplementsMixin.hpp
+++ b/corsika/framework/utility/ImplementsMixin.hpp
@@ -37,7 +37,7 @@ namespace corsika {
};
template <template <typename> typename Mixin, typename T>
- typedef implements_mixin<Mixin, T>::value implements_mixin_v;
+ bool constexpr implements_mixin_v = implements_mixin<Mixin, T>::value;
} // namespace detail
} // namespace corsika
diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
index a5c1ab418..3b794a34f 100644
--- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+++ b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
@@ -68,7 +68,7 @@ namespace corsika {
protected:
LayeredSphericalAtmosphereBuilder(TModelArgs... args, Point const& center,
- LengthType planetRadius)
+ LengthType const planetRadius)
: center_(center)
, planetRadius_(planetRadius)
, additionalModelArgs_{args...} {}
@@ -78,9 +78,11 @@ namespace corsika {
typedef typename VolumeTreeNode<TMediumInterface>::VTNUPtr volume_tree_node_uptr;
void setNuclearComposition(NuclearComposition const& composition);
- volume_tree_node* addExponentialLayer(GrammageType const b, LengthType const scaleHeight,
+ volume_tree_node* addExponentialLayer(GrammageType const b,
+ LengthType const scaleHeight,
LengthType const upperBoundary);
- void addLinearLayer(LengthType const c, LengthType const upperBoundary);
+ void addLinearLayer(GrammageType const b, LengthType const scaleHeight,
+ LengthType const upperBoundary);
void addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
unsigned int const nBins, LengthType const deltaHeight,
@@ -97,7 +99,7 @@ namespace corsika {
LengthType getPlanetRadius() const { return planetRadius_; }
private:
- void checkRadius(LengthType r) const;
+ void checkRadius(LengthType const r) const;
std::unique_ptr<NuclearComposition> composition_;
Point center_;
diff --git a/corsika/media/USStandardAtmosphere.hpp b/corsika/media/USStandardAtmosphere.hpp
index 5b228598f..bc353c795 100644
--- a/corsika/media/USStandardAtmosphere.hpp
+++ b/corsika/media/USStandardAtmosphere.hpp
@@ -24,48 +24,147 @@ namespace corsika {
*/
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,*/
+ 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;
- };
+ namespace detail {
+ struct AtmosphereLayerParameters {
+ LengthType altitude;
+ GrammageType offset;
+ LengthType scaleHeight;
+ };
+
+ typedef std::array<AtmosphereLayerParameters, 5> AtmosphereParameters;
- typedef std::array<AtmosphereLayerParameters, 5> AtmosphereParameters;
+ auto constexpr grammage(double const v) { return v * 1_g / (1_cm * 1_cm); }
- 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}}};
+ std::array<AtmosphereParameters,
+ static_cast<uint8_t>(
+ AtmosphereId::LastAtmosphere)> constexpr atmosphereParameterList{
+ {{{{4_km, grammage(1222.6562), 994186.38_cm},
+ {10_km, grammage(1144.9069), 878153.55_cm},
+ {40_km, grammage(1305.5948), 636143.04_cm},
+ {100_km, grammage(540.1778), 772170.16_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1173.9861), 919546_cm},
+ {10_km, grammage(1205.7625), 963267.92_cm},
+ {40_km, grammage(1386.7807), 614315_cm},
+ {100_km, grammage(555.8935), 739059.6_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1240.48), 933697_cm},
+ {10_km, grammage(1117.85), 765229_cm},
+ {40_km, grammage(1210.9), 636790_cm},
+ {100_km, grammage(608.2128), 733793.8_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1285.2782), 1088310_cm},
+ {10_km, grammage(1173.1616), 935485_cm},
+ {40_km, grammage(1320.4561), 635137_cm},
+ {100_km, grammage(680.6803), 727312.6_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1251.474), 1032310_cm},
+ {10_km, grammage(1173.321), 925528_cm},
+ {40_km, grammage(1307.826), 645330_cm},
+ {100_km, grammage(763.1139), 720851.4_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(113.3362), 923077_cm},
+ {10_km, grammage(1226.5761), 1109960_cm},
+ {40_km, grammage(1382.6933), 630217_cm},
+ {100_km, grammage(685.6073), 726901.3_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1262.7013), 1059360_cm},
+ {10_km, grammage(1139.0249), 888814_cm},
+ {10_km, grammage(1270.2886), 639902_cm},
+ {40_km, grammage(681.4061), 727251.8_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1210.4), 970276_cm},
+ {10_km, grammage(1103.8629), 820946_cm},
+ {40_km, grammage(1215.3545), 639074_cm},
+ {100_km, grammage(629.7611), 731776.5_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{4_km, grammage(1130.74), 867358_cm},
+ {10_km, grammage(1052.05), 741208_cm},
+ {40_km, grammage(1137.21), 633846_cm},
+ {100_km, grammage(442.512), 759850_cm},
+ {112.8_km, grammage(1), 5.4303203e9_cm}}},
+ {{{4_km, grammage(1183.70), 875221_cm},
+ {10_km, grammage(1108.06), 753212_cm},
+ {40_km, grammage(1424.02), 545846_cm},
+ {100_km, grammage(207.595), 793043_cm},
+ {112.8_km, grammage(1), 5.9787908e9_cm}}},
+ {{{4_km, grammage(1177.19), 861745_cm},
+ {10_km, grammage(1125.11), 765925_cm},
+ {40_km, grammage(1304.77), 581351_cm},
+ {100_km, grammage(433.823), 775155_cm},
+ {112.8_km, grammage(1), 7.4095699e9_cm}}},
+ {{{4_km, grammage(1139.99), 861913_cm},
+ {10_km, grammage(1073.82), 744955_cm},
+ {40_km, grammage(1052.96), 675928_cm},
+ {100_km, grammage(492.503), 829627_cm},
+ {112.8_km, grammage(1), 5.8587010e9_cm}}},
+ {{{2.67_km, grammage(1133.10), 861730_cm},
+ {5.33_km, grammage(1101.20), 826340_cm},
+ {8_km, grammage(1085.00), 790950_cm},
+ {100_km, grammage(1098.00), 682800_cm},
+ {112.8_km, grammage(1), 26798156e9_cm}}},
+ {{{6.67_km, grammage(1079.00), 764170_cm},
+ {13.33_km, grammage(1071.90), 699910_cm},
+ {20_km, grammage(1182.00), 635650_cm},
+ {100_km, grammage(1647.10), 551010_cm},
+ {112.8_km, grammage(1), 59.329575e9_cm}}},
+ {{{8_km, grammage(1198.5972), 945766.30_cm},
+ {18.1_km, grammage(1198.8796), 681780.12_cm},
+ {34.5_km, grammage(1419.4152), 620224.52_cm},
+ {100_km, grammage(730.6380), 728157.92_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{8.3_km, grammage(1179.5010), 939228.66_cm},
+ {12.9_km, grammage(1172.4883), 787969.34_cm},
+ {34_km, grammage(1437.4911), 620008.53_cm},
+ {100_km, grammage(761.3281), 724585.33_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{5.9_km, grammage(1202.8804), 977139.52_cm},
+ {12.0_km, grammage(1148.6275), 858087.01_cm},
+ {34.5_km, grammage(1432.0312), 614451.60_cm},
+ {100_km, grammage(696.42788), 730875.73_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{9_km, grammage(1175.3347), 986169.72_cm},
+ {14.6_km, grammage(1180.3694), 793171.45_cm},
+ {33_km, grammage(1614.5404), 600120.95_cm},
+ {100_km, grammage(755.56438), 725247.87_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{8_km, grammage(1196.9290), 985241.1_cm},
+ {13_km, grammage(1173.2537), 819245.00_cm},
+ {33.5_km, grammage(1502.1837), 611220.86_cm},
+ {100_km, grammage(750.89704), 725797.06_cm},
+ {112.8_km, grammage(1), 1e9_cm}}},
+ {{{7_km, grammage(1183.6071), 954248.34_cm},
+ {11.4_km, grammage(1143.0425), 800005.34_cm},
+ {37_km, grammage(1322.9748), 629568.93_cm},
+ {100_km, grammage(655.67307), 737521.77_cm},
+ {112.8_km, grammage(1), 1e9_cm}}}}};
+ } // namespace detail
template <typename TEnvironmentInterface, template <typename> typename TExtraEnv,
typename TEnvironment, typename... TArgs>
- auto create_5layer_atmosphere(AtmosphereId const atmId, TEnvironment& env,
+ auto create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
Point const& center, TArgs... args) {
// construct the atmosphere builder
@@ -79,22 +178,16 @@ namespace corsika {
{{Code::Nitrogen, Code::Oxygen}, {0.7847, 1. - 0.7847}});
// add the standard atmosphere layers
- auto params = AtmosphereParameterList[static_cast<uint8_t>(atmId)];
+ auto const params = detail::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.addExponentialLayer(params[i].offset, params[i].scaleHeight,
+ params[i].altitude);
}
- 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);*/
+ builder.addLinearLayer(params[4].offset, params[4].scaleHeight, params[4].altitude);
// check if we want to also add the US standard refractivity
if constexpr (detail::implements_mixin_v<IRefractiveIndexModel,
TEnvironmentInterface>) {
-
// TODO: Add US Standard refractivity
}
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 5bc1d5fcd..925824a8e 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -123,23 +123,11 @@ int main(int argc, char** argv) {
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
- auto builder = make_layered_spherical_atmosphere_builder<
- setup::EnvironmentInterface, MyExtraEnv>::create(center,
- constants::EarthRadius::Mean,
- Medium::AirDry1Atm,
- MagneticFieldVector{rootCS, 0_T,
- 50_uT, 0_T});
- builder.setNuclearComposition(
- {{Code::Nitrogen, Code::Oxygen},
- {0.7847, 1. - 0.7847}}); // values taken from AIRES manual, Ar removed for now
-
- 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);
- 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.assemble(env);
+
+ // build a Linsley US Standard atmosphere into `env`
+ create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>(
+ env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
+ MagneticFieldVector{rootCS, 0_T, 50_uT, 0_T});
// pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1]));
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 372cfdb9a..c9d6a673d 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -555,11 +555,11 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
builder.setNuclearComposition({{{Code::Nitrogen, Code::Oxygen}}, {{.6, .4}}});
- builder.addLinearLayer(1_km, 10_km);
- builder.addLinearLayer(2_km, 20_km);
+ builder.addLinearLayer(1_g / (1_cm * 1_cm), 1_km, 10_km);
+ builder.addLinearLayer(1_g / (1_cm * 1_cm), 2_km, 20_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 30_km);
- CHECK_THROWS(builder.addLinearLayer(0.5_km, 5_km));
+ CHECK_THROWS(builder.addLinearLayer(1_g / (1_cm * 1_cm), 0.5_km, 5_km));
CHECK(builder.getSize() == 3);
@@ -603,7 +603,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
B0);
builder.setNuclearComposition({{{Code::Nitrogen, Code::Oxygen}}, {{.6, .4}}});
- builder.addLinearLayer(1_km, 10_km);
+ builder.addLinearLayer(1_g / (1_cm * 1_cm), 1_km, 10_km);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 20_km);
CHECK(builder.getSize() == 2);
@@ -644,7 +644,7 @@ TEST_CASE("media", "LayeredSphericalAtmosphereBuilder USStd") {
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.addLinearLayer(1_g / (1_cm * 1_cm), 1e9_cm, 112.8_km);
Environment<IMediumModel> env;
builder.assemble(env);
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 985043d4d..16539ba56 100644
--- a/tests/modules/testCONEX.cpp
+++ b/tests/modules/testCONEX.cpp
@@ -78,7 +78,7 @@ TEST_CASE("CONEXSourceCut") {
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.addLinearLayer(1_g / (1_cm * 1_cm), 1e9_cm, 112.8_km);
builder.assemble(env);
--
GitLab