diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl
index dbd4180410675b3c61c37438737ccb732db30768..edb849fd90862c608be6d8f5f43dc1e03bbe6dd1 100644
--- a/corsika/detail/media/BaseExponential.inl
+++ b/corsika/detail/media/BaseExponential.inl
@@ -17,8 +17,8 @@ namespace corsika {
template <typename TDerived>
inline BaseExponential<TDerived>::BaseExponential(Point const& point,
LengthType const referenceHeight,
- MassDensityType rho0,
- LengthType lambda)
+ MassDensityType const rho0,
+ LengthType const lambda)
: rho0_(rho0)
, lambda_(lambda)
, invLambda_(1 / lambda)
@@ -43,7 +43,7 @@ namespace corsika {
if (length == LengthType::zero()) { return GrammageType::zero(); }
// this corresponds to height:
- double const uDotA = traj.getDirection(0).dot(axis).magnitude();
+ double const uDotA = traj.getDirection(0).dot(axis);
MassDensityType const rhoStart =
getImplementation().getMassDensity(traj.getPosition(0));
@@ -56,10 +56,10 @@ namespace corsika {
template <typename TDerived>
inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
- BaseTrajectory const& traj, GrammageType grammage,
+ BaseTrajectory const& traj, GrammageType const grammage,
DirectionVector const& axis) const {
// this corresponds to height:
- double const uDotA = traj.getDirection(0).dot(axis).magnitude();
+ double const uDotA = traj.getDirection(0).dot(axis);
MassDensityType const rhoStart =
getImplementation().getMassDensity(traj.getPosition(0));
diff --git a/corsika/detail/media/CORSIKA7Atmospheres.inl b/corsika/detail/media/CORSIKA7Atmospheres.inl
new file mode 100644
index 0000000000000000000000000000000000000000..1e18dffb35912ff2c0e472accf82f562dd169c82
--- /dev/null
+++ b/corsika/detail/media/CORSIKA7Atmospheres.inl
@@ -0,0 +1,43 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+ template <typename TEnvironmentInterface, template <typename> typename TExtraEnv,
+ typename TEnvironment, typename... TArgs>
+ void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
+ Point const& center, TArgs... args) {
+
+ // construct the atmosphere builder
+ auto builder = make_layered_spherical_atmosphere_builder<
+ TEnvironmentInterface, TExtraEnv>::create(center, constants::EarthRadius::Mean,
+ std::forward<TArgs>(args)...);
+
+ // composition values from AIRES manual
+ builder.setNuclearComposition({{
+ Code::Nitrogen,
+ Code::Argon,
+ Code::Oxygen,
+ },
+ {0.7847, 0.0047, 1. - 0.7847 - 0.0047}});
+
+ // add the standard atmosphere layers
+ auto const params = atmosphereParameterList[static_cast<uint8_t>(atmId)];
+ for (int i = 0; i < 4; ++i) {
+ builder.addExponentialLayer(params[i].offset, params[i].scaleHeight,
+ params[i].altitude);
+ }
+ builder.addLinearLayer(params[4].offset, params[4].scaleHeight, params[4].altitude);
+
+ // and assemble the environment
+ builder.assemble(env);
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl
index f665dda992148f252261327c5a7f931cf5a9416c..ea0fe5b0e8bd96e0c9a68a0680bec2b1b0f7590c 100644
--- a/corsika/detail/media/FlatExponential.inl
+++ b/corsika/detail/media/FlatExponential.inl
@@ -19,7 +19,8 @@ namespace corsika {
template <typename T>
inline FlatExponential<T>::FlatExponential(Point const& point,
DirectionVector const& axis,
- MassDensityType rho, LengthType lambda,
+ MassDensityType const rho,
+ LengthType const lambda,
NuclearComposition const& nuclComp)
: BaseExponential<FlatExponential<T>>(point, 0_m, rho, lambda)
, axis_(axis)
@@ -44,7 +45,7 @@ namespace corsika {
template <typename T>
inline LengthType FlatExponential<T>::getArclengthFromGrammage(
- BaseTrajectory const& line, GrammageType grammage) const {
+ BaseTrajectory const& line, GrammageType const grammage) const {
return BaseExponential<FlatExponential<T>>::getArclengthFromGrammage(line, grammage,
axis_);
}
diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index 0f17cf7c8ddf03f5a9e036038240e8a786ed0f9c..5e697b51c7d379895b8a879759a4f4e2da9ea5cb 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
@@ -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,8 +40,10 @@ 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 +51,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 +66,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 +77,9 @@ namespace corsika {
typename... TModelArgs>
inline void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
- TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
+ TModelArgs...>::addLinearLayer(GrammageType const b, LengthType const scaleHeight,
+ LengthType const upperBoundary) {
+ // outer radius
auto const radius = planetRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
@@ -83,8 +87,7 @@ 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;
+ 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
@@ -167,7 +170,8 @@ namespace corsika {
template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
struct make_layered_spherical_atmosphere_builder {
template <typename... TArgs>
- static auto create(Point const& center, LengthType planetRadius, TArgs... args) {
+ static auto create(Point const& center, LengthType const planetRadius,
+ TArgs... args) {
return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
TArgs...>{std::forward<TArgs>(args)...,
center, planetRadius};
diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl
index 9b921872811e157c9d0d232286b3714444983727..b9c1728484ae622f117cbd628edbdf0f4fbc6e9a 100644
--- a/corsika/detail/media/SlidingPlanarExponential.inl
+++ b/corsika/detail/media/SlidingPlanarExponential.inl
@@ -12,8 +12,8 @@ namespace corsika {
template <typename TDerived>
inline SlidingPlanarExponential<TDerived>::SlidingPlanarExponential(
- Point const& p0, MassDensityType rho0, LengthType lambda,
- NuclearComposition const& nuclComp, LengthType referenceHeight)
+ Point const& p0, MassDensityType const rho0, LengthType const lambda,
+ NuclearComposition const& nuclComp, LengthType const referenceHeight)
: BaseExponential<SlidingPlanarExponential<TDerived>>(p0, referenceHeight, rho0,
lambda)
, nuclComp_(nuclComp) {}
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index ed42473ec865d339e540e7b7e1904df763136e09..07f390b7385a6be73b155b8fe8dd8a64b4308c7a 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -127,7 +127,10 @@ namespace corsika::sibyll {
si::CrossSectionType weightedProdCrossSection = mediumComposition.getWeightedSum(
[=](corsika::Code targetID) -> si::CrossSectionType {
- return std::get<0>(this->getCrossSection(corsikaBeamId, targetID, ECoM));
+ // Argon needs special handling ....
+ return targetID == Code::Argon ? CrossSectionType::zero()
+ : std::get<0>(this->getCrossSection(
+ corsikaBeamId, targetID, ECoM));
});
CORSIKA_LOG_DEBUG(
@@ -243,6 +246,7 @@ namespace corsika::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
+ if (targetId == Code::Argon) continue; // skip Argon ....
const auto [sigProd, sigEla] = getCrossSection(corsikaBeamId, targetId, Ecm);
[[maybe_unused]] const auto& dummy_sigEla = sigEla;
cross_section_of_components[i] = sigProd;
diff --git a/corsika/detail/modules/sibyll/NuclearInteraction.inl b/corsika/detail/modules/sibyll/NuclearInteraction.inl
index db7788a117dd5a6cb5cb75422a10e37ebee0cffb..1460b1133616c1347cdd2eaa156bbcd06ac307b1 100644
--- a/corsika/detail/modules/sibyll/NuclearInteraction.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteraction.inl
@@ -51,6 +51,10 @@ namespace corsika::sibyll {
template <typename TEnvironment>
inline void NuclearInteraction<TEnvironment>::printCrossSectionTable(Code pCode) {
+ if (pCode == Code::Argon) {
+ CORSIKA_LOG_WARN("SIBYLL cannot handle Argon as target!");
+ return;
+ }
const int k = targetComponentsIndex_.at(pCode);
Code pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
Code::Neon, Code::Argon, Code::Iron};
@@ -97,6 +101,7 @@ namespace corsika::sibyll {
// loop over target components, at most 4!!
int k = -1;
for (auto& ptarg : allElementsInUniverse) {
+ if (ptarg == Code::Argon) continue; // NEED TO IGNORE Argon ....
++k;
CORSIKA_LOG_DEBUG("NuclearInteraction: init target component: {}", ptarg);
const int ib = get_nucleus_A(ptarg);
@@ -265,6 +270,7 @@ namespace corsika::sibyll {
const auto& w = mediumComposition.getFractions();
// loop over components in medium
for (auto const targetId : mediumComposition.getComponents()) {
+ if (targetId == Code::Argon) continue; // NEED TO IGNORE Argon ....
i++;
CORSIKA_LOG_DEBUG("NuclearInteraction: get interaction length for target: {}",
get_name(targetId));
@@ -422,6 +428,7 @@ namespace corsika::sibyll {
for (size_t i = 0; i < compVec.size(); ++i) {
auto const targetId = compVec[i];
+ if (targetId == Code::Argon) continue; // NEED TO IGNORE Argon ....
CORSIKA_LOG_DEBUG("target component: {}", get_name(targetId));
CORSIKA_LOG_DEBUG("beam id: {}", get_name(beamId));
const auto [sigProd, sigEla] =
diff --git a/corsika/framework/utility/ImplementsMixin.hpp b/corsika/framework/utility/ImplementsMixin.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9082a9bcd1ca44b9415d5c304d0efcc26608f13b
--- /dev/null
+++ b/corsika/framework/utility/ImplementsMixin.hpp
@@ -0,0 +1,43 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+#pragma once
+
+/**
+ * @file ImplementsMixin.hpp
+ */
+
+#include <type_traits>
+
+namespace corsika {
+
+ namespace detail {
+
+ /**
+ Helper traits class (partial) for static compile time checking.
+
+ This method checks whether a given class implements a particular
+ mixin - this is particularly useful in the media heirarchy that utilizes
+ mixins for the interface definition.
+ */
+ template <template <typename> typename Mixin, typename T>
+ struct implements_mixin {
+
+ template <typename U>
+ static std::true_type test(Mixin<U>&);
+ static std::false_type test(...);
+
+ public:
+ using type = decltype(test(std::declval<T&>()));
+ static constexpr bool value = type::value;
+ };
+
+ template <template <typename> typename Mixin, typename T>
+ bool constexpr implements_mixin_v = implements_mixin<Mixin, T>::value;
+
+ } // namespace detail
+} // namespace corsika
diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp
index 507ed640fa130eda007a5e71a3dec3437a530d9b..7e8aae9f7526d314a247275a53fc7a009d31da18 100644
--- a/corsika/media/BaseExponential.hpp
+++ b/corsika/media/BaseExponential.hpp
@@ -15,6 +15,10 @@
#include <corsika/framework/geometry/BaseTrajectory.hpp>
#include <limits>
+/**
+ * @file corsika/media/BaseExponential.hpp
+ */
+
namespace corsika {
/**
@@ -26,7 +30,7 @@ namespace corsika {
public:
BaseExponential(Point const& point, LengthType const referenceHeight,
- MassDensityType rho0, LengthType lambda);
+ MassDensityType const rho0, LengthType const lambda);
Point const& getAnchorPoint() const { return point_; }
@@ -38,14 +42,15 @@ namespace corsika {
// clang-format off
/**
* For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
- * direction \f$ \vec{u} \f$ is given by
+ * direction \f$ \vec{u} \f$ is given by:
* \f[
- * X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right)
- * \f], where \f$ \varrho_0 \f$ is the density at the starting point.
+ * X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right) \quad \text{,}
+ * \f]
+ * where \f$ \varrho_0 \f$ is the density at the starting point.
*
* If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[
- * X = \varrho_0 l;
+ * X = \varrho_0 l
* \f]
*/
// clang-format on
@@ -55,22 +60,23 @@ namespace corsika {
// clang-format off
/**
* For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized)
- * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by
+ * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by:
* \f[
* l = \begin{cases}
- * \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 +
- * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda}
- * \infty & \text{else,}
+ * \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} & Y := 1 +
+ * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda} > 0 \\
+ * \infty & \text{else} & \text{,}
* \end{cases}
* \f] where \f$ \varrho_0 \f$ is the density at the starting point.
*
- * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
+ * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like for a homogeneous density:
* \f[
* l = \frac{X}{\varrho_0}
* \f]
*/
// clang-format on
- LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage,
+ LengthType getArclengthFromGrammage(BaseTrajectory const& line,
+ GrammageType const grammage,
DirectionVector const& axis) const;
private:
diff --git a/corsika/media/CORSIKA7Atmospheres.hpp b/corsika/media/CORSIKA7Atmospheres.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0eb0b6a7a38d3dc8b4a9de238fa9c47f6555c5d8
--- /dev/null
+++ b/corsika/media/CORSIKA7Atmospheres.hpp
@@ -0,0 +1,209 @@
+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/media/IRefractiveIndexModel.hpp>
+#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
+#include <corsika/framework/utility/ImplementsMixin.hpp>
+
+// for detail namespace, NoExtraModelInner, NoExtraModel and traits
+#include <corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp>
+
+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
+ };
+
+ namespace {
+ /**
+ * Struct to hold parameters of one layer of a CORSIKA7 atmosphere.
+ *
+ * The definition of each layer is according to a BaseExponential:
+ * @f[
+ * \varrho = offset/scaleHeight \cdot
+ * \exp\left(-(height-altitude)/scaleHeight\right)
+ * @f],
+ * where @f$ altitude @f$ is the height where the atmosphere layer starts, @f$
+ * offset/scaleHeight
+ * @f$ is the density at this height.
+ */
+ struct AtmosphereLayerParameters {
+ LengthType altitude;
+ GrammageType offset;
+ LengthType scaleHeight;
+ };
+
+ /**
+ * All the 5 layers of a CORSIKA7 atmosphere.
+ */
+ typedef std::array<AtmosphereLayerParameters, 5> AtmosphereParameters;
+
+ /**
+ * Local, internal helper function to provide "Grammage" type.
+ *
+ * @param v
+ * @return value v with g/cm2 units
+ */
+ auto constexpr grammage(double const v) { return v * 1_g / (1_cm * 1_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
+
+ /**
+ * Function to create a CORSIKA 7 5-layer atmosphere.
+ *
+ * @tparam TEnvironmentInterface
+ * @tparam TExtraEnv
+ * @tparam TEnvironment
+ * @tparam TArgs
+ * @param env
+ * @param atmId
+ * @param center
+ * @param args
+ */
+ template <typename TEnvironmentInterface,
+ template <typename> typename TExtraEnv = detail::NoExtraModel,
+ typename TEnvironment, typename... TArgs>
+ void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId,
+ Point const& center, TArgs... args);
+
+} // namespace corsika
+
+#include <corsika/detail/media/CORSIKA7Atmospheres.inl>
\ No newline at end of file
diff --git a/corsika/media/Environment.hpp b/corsika/media/Environment.hpp
index a52aeaebe3c33548f52db8aeddcc096940ff9307..5d5d67c94d36b7c7a7bc3f2884437d9f9871ce6d 100644
--- a/corsika/media/Environment.hpp
+++ b/corsika/media/Environment.hpp
@@ -20,9 +20,11 @@
namespace corsika {
- /** Base Evnironment class
- * Describes the Environment in which the shower is propagated
- **/
+ /**
+ * Base Environment class.
+ *
+ * Describes the Environment in which the shower is propagated.
+ */
template <typename IEnvironmentModel>
class Environment {
public:
@@ -30,30 +32,34 @@ namespace corsika {
Environment();
- /** Getters for the universe stored in the Environment
+ /**
+ * Getters for the universe stored in the Environment.
*
- * @retval Retuns reference to a Universe object with infinite size
- **/
+ * @retval Retuns reference to a Universe object with infinite size.
+ */
///@{
- //* Get non const universe */
+ //! Get non const universe
typename BaseNodeType::VTNUPtr& getUniverse();
- //* Get const universe */
+ //! Get const universe
typename BaseNodeType::VTNUPtr const& getUniverse() const;
///@}
- /** Getter for the CoordinateSystem used in the Environment
+ /**
+ * Getter for the CoordinateSystem used in the Environment.
*
- * @retval Retuns a const reference to the CoordinateSystem used
- **/
+ * @retval Retuns a const reference to the CoordinateSystem used.
+ */
CoordinateSystemPtr const& getCoordinateSystem() const;
- /** Factory method for creation of VolumeTreeNodes
+ /**
+ * Factory method for creation of VolumeTreeNodes.
+ *
* @tparam TVolumeType Type of volume to be created
* @tparam TVolumeArgs Types to forward to the constructor
* @param args Parameter forwarded to the constructor of TVolumeType
- * @retval Retuns unique pointer to a VolumeTreeNode with the same EnvitonmentModel as
- *this class
- **/
+ * @retval Returns unique pointer to a VolumeTreeNode with the same EnvitonmentModel
+ * as this class.
+ */
template <class TVolumeType, typename... TVolumeArgs>
static std::unique_ptr<BaseNodeType> createNode(TVolumeArgs&&... args);
diff --git a/corsika/media/FlatExponential.hpp b/corsika/media/FlatExponential.hpp
index 6f64295e83efcc36a86ce4e43cc94a37e6593d78..042273f43d1fe13f0e1d081b1e1c3bdbcb1c5da8 100644
--- a/corsika/media/FlatExponential.hpp
+++ b/corsika/media/FlatExponential.hpp
@@ -21,11 +21,14 @@ namespace corsika {
/**
* flat exponential density distribution with
* \f[
- * \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
+ * \varrho(\vec{r}) = \varrho_0 \exp\left( \frac{1}{\lambda} (\vec{r} - \vec{p}) \cdot
* \vec{a} \right).
* \f]
* \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
- * with the scale parameter \f$ \lambda \f$.
+ * with the scale parameter \f$ \lambda \f$, \f$ \vec{r} \f$ is the location of
+ * the evaluation, \f$ \vec{p} \f$ is the anchor point at which \f$ \varrho_0 \f$
+ * is given. Thus, the unit vector \f$ \vec{a} \f$ specifies the direction of
+ * <em>decreasing</em> <b>height/altitude</b>.
*/
// clang-format on
template <typename T>
@@ -34,7 +37,7 @@ namespace corsika {
public:
FlatExponential(Point const& point, Vector<dimensionless_d> const& axis,
- MassDensityType rho, LengthType lambda,
+ MassDensityType const rho, LengthType const lambda,
NuclearComposition const& nuclComp);
MassDensityType getMassDensity(Point const& point) const override;
@@ -44,7 +47,7 @@ namespace corsika {
GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& line,
- GrammageType grammage) const override;
+ GrammageType const grammage) const override;
private:
DirectionVector const axis_;
diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
index 09bcabc06472717a79ed14f62778843aac042f4f..dfbb1504390440042a181260461b3b65903dd586 100644
--- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+++ b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
@@ -24,10 +24,14 @@
#include <type_traits>
#include <functional>
+/**
+ * @file LayeredSphericalAtmosphereBuilder.hpp
+ */
+
namespace corsika {
/**
- * \class make_layered_spherical_atmosphere_builder
+ * make_layered_spherical_atmosphere_builder.
*
* Helper class to create LayeredSphericalAtmosphereBuilder, the
* extra environment models have to be passed as template-template
@@ -35,7 +39,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 +53,6 @@ namespace corsika {
*
* Each layer by definition has a density profile and a (constant)
* nuclear composition model.
- *
*/
template <typename TMediumInterface = IMediumModel,
@@ -69,7 +72,7 @@ namespace corsika {
protected:
LayeredSphericalAtmosphereBuilder(TModelArgs... args, Point const& center,
- LengthType planetRadius)
+ LengthType const planetRadius)
: center_(center)
, planetRadius_(planetRadius)
, additionalModelArgs_{args...} {}
@@ -79,9 +82,11 @@ 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(GrammageType const b, LengthType const scaleHeight,
+ LengthType const upperBoundary);
void addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
unsigned int const nBins, LengthType const deltaHeight,
@@ -98,7 +103,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/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp
index 4ef91bb4e1ca338241c61d805f1fa3dd68c47ba4..3e1b74836fd44df3d7c02becc3fd6054890f524d 100644
--- a/corsika/media/SlidingPlanarExponential.hpp
+++ b/corsika/media/SlidingPlanarExponential.hpp
@@ -20,7 +20,7 @@ namespace corsika {
// clang-format off
/**
- * The SlidingPlanarExponential models mass density as
+ * The SlidingPlanarExponential models mass density as:
* \f[
* \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right).
* \f]
@@ -37,9 +37,9 @@ namespace corsika {
using Base = BaseExponential<SlidingPlanarExponential<T>>;
public:
- SlidingPlanarExponential(Point const& p0, MassDensityType rho0, LengthType lambda,
- NuclearComposition const& nuclComp,
- LengthType referenceHeight = LengthType::zero());
+ SlidingPlanarExponential(Point const& p0, MassDensityType const rho0,
+ LengthType const lambda, NuclearComposition const& nuclComp,
+ LengthType const referenceHeight = LengthType::zero());
MassDensityType getMassDensity(Point const& point) const override;
@@ -48,7 +48,7 @@ namespace corsika {
GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& line,
- GrammageType grammage) const override;
+ GrammageType const grammage) const override;
private:
NuclearComposition const nuclComp_;
diff --git a/documentation/Doxyfile.in b/documentation/Doxyfile.in
index f3768a1d22014397db7d836339de10dd020695ec..f942b9ee098840314ca2a4aab5dc81927f982c9b 100644
--- a/documentation/Doxyfile.in
+++ b/documentation/Doxyfile.in
@@ -1,5 +1,5 @@
-PROJECT_NAME = CORSIKA 8
-PROJECT_NUMBER = 0.0.0
+PROJECT_NAME = CORSIKA
+PROJECT_NUMBER = @c8_version@
PROJECT_BRIEF = "The framework to simulate particle cascades for astroparticle physics"
GENERATE_HTML = YES
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index f57ec5f133749f34122636710de35268eea8ce2f..a3db1ccfa0c4d35899c25d9603e9655c69a7443d 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -37,7 +37,7 @@
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
@@ -158,24 +158,29 @@ int main(int argc, char** argv) {
->default_str("info")
->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
->group("Misc.");
+ app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.")
+ ->default_val("info")
+ ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
+ ->group("Misc.");
// parse the command line options into the variables
CLI11_PARSE(app, argc, argv);
- string const loglevel =
- (app.count("--verbosity") ? app["--verbosity"]->as<string>() : "info");
- if (loglevel == "warn") {
- logging::set_level(logging::level::warn);
- } else if (loglevel == "info") {
- logging::set_level(logging::level::info);
- } else if (loglevel == "debug") {
- logging::set_level(logging::level::debug);
- } else if (loglevel == "trace") {
+ if (app.count("--verbosity")) {
+ string const loglevel = app["verbosity"]->as<string>();
+ if (loglevel == "warn") {
+ logging::set_level(logging::level::warn);
+ } else if (loglevel == "info") {
+ logging::set_level(logging::level::info);
+ } else if (loglevel == "debug") {
+ logging::set_level(logging::level::debug);
+ } else if (loglevel == "trace") {
#ifndef DEBUG
- CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
- return 1;
+ CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
+ return 1;
#endif
- logging::set_level(logging::level::trace);
+ logging::set_level(logging::level::trace);
+ }
}
// check that we got either PDG or A/Z
@@ -197,28 +202,17 @@ 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);
- 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});
+
/* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
ofstream atmout("earth.dat");
- for (LengthType h = 0_m; h < 110_km; h += 10_m) {
- Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h};
+ for (LengthType h = 0_m; h < 110_km; h += 100_m) {
+ Point const ptest{rootCS, 0_m, 0_m, constants::EarthRadius::Mean + h};
auto rho =
env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity(
ptest);
@@ -263,8 +257,8 @@ int main(int argc, char** argv) {
/* === END: CONSTRUCT PRIMARY PARTICLE === */
/* === START: CONSTRUCT GEOMETRY === */
- auto const observationHeight = 0_km + builder.getPlanetRadius();
- auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
+ auto const observationHeight = 0_km + constants::EarthRadius::Mean;
+ auto const injectionHeight = 111.75_km + constants::EarthRadius::Mean;
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 4730561718a80dd962729670a726cf71407ae387..62aa91e564010d12a0fbcb3ab48062917ca42da4 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -26,6 +26,7 @@
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
#include <corsika/modules/ObservationPlane.hpp>
@@ -88,21 +89,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,
- Vector{rootCS, 0_T, 50_uT, 0_T});
- builder.setNuclearComposition(
- {{Code::Nitrogen, 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);
- 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);
+
+ // 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});
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -128,8 +119,8 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
- auto const observationHeight = 1.4_km + builder.getPlanetRadius();
- auto const injectionHeight = 112.75_km + builder.getPlanetRadius();
+ auto const observationHeight = 1.4_km + constants::EarthRadius::Mean;
+ auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index dabed46e9f4e598b4146eebd673cc9ad07721fd0..4a85351b0e9e5ebdb4c24910a9fa77655c14215a 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -34,6 +34,7 @@
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
@@ -138,21 +139,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,
- Vector{rootCS, 0_T, 50_uT, 0_T});
- builder.setNuclearComposition(
- {{Code::Nitrogen, 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);
- 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);
+
+ // 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});
// setup particle stack, and add primary particle
setup::Stack stack;
@@ -180,8 +171,8 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
- auto const observationHeight = 0_km + builder.getPlanetRadius();
- auto const injectionHeight = 112.75_km + builder.getPlanetRadius();
+ auto const observationHeight = 0_km + constants::EarthRadius::Mean;
+ auto const injectionHeight = 112.75_km + constants::EarthRadius::Mean;
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/examples/mars.cpp b/examples/mars.cpp
index e1f525020a0712b7aca32dee4449e409794be082..14eab4096b2c07aa34d65ffa7fafdc3e12d3f1e7 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -241,7 +241,7 @@ int main(int argc, char** argv) {
builder.addTabularLayer(layer1, 100, 100_m, 7_km);
builder.addTabularLayer(layer2, 300, 500_m, 100_km);
- builder.addLinearLayer(1e9_cm, 112.8_km);
+ builder.addLinearLayer(1_g / square(1_cm), 1e9_cm, 112.8_km);
builder.assemble(env);
/* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 09e229bca9df09409e1146e34fcf505a11c47120..911f45c528f6f440c45b7a94048c2e3478b5a8a3 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -34,12 +34,11 @@
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
-#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/MediumPropertyModel.hpp>
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
-#include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
@@ -123,23 +122,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]));
@@ -180,8 +167,8 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
- auto const observationHeight = 0_km + builder.getPlanetRadius();
- auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
+ auto const observationHeight = 0_km + constants::EarthRadius::Mean;
+ auto const injectionHeight = 111.75_km + constants::EarthRadius::Mean;
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/tests/media/CMakeLists.txt b/tests/media/CMakeLists.txt
index 292b5793e02fc47efa55b9089f9ebf1efe77dbc4..ffcbe11d185219c6c35e34d91e4a4dcf48b75794 100644
--- a/tests/media/CMakeLists.txt
+++ b/tests/media/CMakeLists.txt
@@ -5,6 +5,7 @@ set (test_media_sources
testMedium.cpp
testRefractiveIndex.cpp
testMagneticField.cpp
+ testCORSIKA7Atmospheres.cpp
)
CORSIKA_ADD_TEST (testMedia SOURCES ${test_media_sources})
diff --git a/tests/media/testCORSIKA7Atmospheres.cpp b/tests/media/testCORSIKA7Atmospheres.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..55b2081675c0385912fa4a4725fd6986929595c4
--- /dev/null
+++ b/tests/media/testCORSIKA7Atmospheres.cpp
@@ -0,0 +1,38 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+
+#include <corsika/media/CORSIKA7Atmospheres.hpp>
+
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+
+TEST_CASE("CORSIKA7Atmospheres") {
+
+ logging::set_level(logging::level::info);
+
+ CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
+ Point const gOrigin(gCS, {0_m, 0_m, 0_m});
+
+ Environment<IMediumModel> env;
+
+ // build a Linsley US Standard atmosphere into `env`
+ create_5layer_atmosphere<IMediumModel>(env, AtmosphereId::LinsleyUSStd, gOrigin);
+
+ typedef typename Environment<IMediumModel>::BaseNodeType::VTN_type node_type;
+ node_type const* universe = env.getUniverse().get();
+
+ Point const p(gCS, {constants::EarthRadius::Mean, 0_m, 0_m});
+ CHECK(universe->getContainingNode(p)->getModelProperties().getMassDensity(p) / 1_g *
+ static_pow<3>(1_cm) ==
+ Approx(1222.6562 / 994186.38));
+}
\ No newline at end of file
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 372cfdb9a578e321354517d3772c812920022fb5..43476644e38aa44b987a0386c2d123a7e630cf67 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -79,8 +79,7 @@ TEST_CASE("HomogeneousMedium") {
logging::set_level(logging::level::info);
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, {1.});
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
CHECK(protonComposition.getFractions() == std::vector<float>{1.});
@@ -94,8 +93,7 @@ TEST_CASE("FlatExponential") {
logging::set_level(logging::level::info);
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
LengthType const lambda = 3_m;
@@ -162,8 +160,7 @@ TEST_CASE("SlidingPlanarExponential") {
logging::set_level(logging::level::info);
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, {1.});
LengthType const lambda = 3_m;
auto const rho0 = 1_g / static_pow<3>(1_cm);
@@ -245,8 +242,7 @@ TEST_CASE("SlidingPlanarTabular") {
logging::set_level(logging::level::info);
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, {1.});
RhoFuncConst rhoFunc;
SlidingPlanarTabular<IMediumModel> const medium(gOrigin, rhoFunc, 1000, 10_m,
@@ -555,11 +551,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);
@@ -592,8 +588,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
using ModelInterface = IMagneticFieldModel<IMediumModel>;
// the composition we use for the homogenous medium
- NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
- std::vector<float>{1.f});
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton}, {1.});
// create magnetic field vectors
Vector B0(gCS, 0_T, 0_T, 1_T);
@@ -603,7 +598,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 +639,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);
@@ -652,7 +647,7 @@ TEST_CASE("media", "LayeredSphericalAtmosphereBuilder USStd") {
typedef typename Environment<IMediumModel>::BaseNodeType::VTN_type node_type;
node_type const* universe = env.getUniverse().get();
- // far out ther is the universe
+ // far out there is the universe
CHECK(universe->getContainingNode(Point(gCS, {10000_km, 0_m, 0_m})) == universe);
CHECK(universe->getContainingNode(Point(gCS, {0_m, 10000_km, 0_m})) == universe);
diff --git a/tests/modules/testCONEX.cpp b/tests/modules/testCONEX.cpp
index 985043d4dca501fab27359ff8be5eb557615ca8f..16539ba56ee241c6afed50bd3748a7ab941314c8 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);