diff --git a/corsika/detail/media/CORSIKA7Atmospheres.inl b/corsika/detail/media/CORSIKA7Atmospheres.inl index 1e18dffb35912ff2c0e472accf82f562dd169c82..3e5c86e76e37fb8c1a90211c2522a9d94b721ef6 100644 --- a/corsika/detail/media/CORSIKA7Atmospheres.inl +++ b/corsika/detail/media/CORSIKA7Atmospheres.inl @@ -20,13 +20,7 @@ namespace corsika { 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}}); + builder.setNuclearComposition(standardAirComposition); // add the standard atmosphere layers auto const params = atmosphereParameterList[static_cast<uint8_t>(atmId)]; diff --git a/corsika/detail/media/NuclearComposition.inl b/corsika/detail/media/NuclearComposition.inl index c499a593d98f9ba064b0656ead63956fc3a78244..9d1241aabfd5ddff1812cf8b0163e5ae7cd13d70 100644 --- a/corsika/detail/media/NuclearComposition.inl +++ b/corsika/detail/media/NuclearComposition.inl @@ -46,7 +46,7 @@ namespace corsika { } template <typename TFunction> - inline auto NuclearComposition::getWeighted(TFunction const& func) const { + inline auto NuclearComposition::getWeighted(TFunction func) const { using ResultQuantity = decltype(func(std::declval<Code>())); auto const product = [&](auto const compID, auto const fraction) { return func(compID) * fraction; @@ -66,7 +66,7 @@ namespace corsika { } // namespace corsika template <typename TFunction> - inline auto NuclearComposition::getWeightedSum(TFunction const& func) const + inline auto NuclearComposition::getWeightedSum(TFunction func) const -> decltype(func(std::declval<Code>())) { using ResultQuantity = decltype(func(std::declval<Code>())); diff --git a/corsika/detail/modules/conex/CONEXhybrid.inl b/corsika/detail/modules/conex/CONEXhybrid.inl index e42ca29fcaf6aa90ddeb4bd88026e68b083de3a9..d2ae2a2b45e628ba4e06d2586a85f8df3bcc97ee 100644 --- a/corsika/detail/modules/conex/CONEXhybrid.inl +++ b/corsika/detail/modules/conex/CONEXhybrid.inl @@ -7,6 +7,7 @@ */ #include <corsika/framework/core/Logging.hpp> +#include <corsika/media/CORSIKA7Atmospheres.hpp> #include <corsika/modules/conex/CONEXhybrid.hpp> #include <corsika/modules/conex/CONEX_f.hpp> #include <corsika/framework/random/RNGManager.hpp> @@ -17,7 +18,7 @@ #include <algorithm> #include <fstream> -#include <iomanip> +#include <numeric> #include <utility> namespace corsika { @@ -83,6 +84,28 @@ namespace corsika { CORSIKA_LOG_DEBUG("showerCore (C8): {}", showerCore_.getCoordinates(center.getCoordinateSystem())); + auto const& components = ::corsika::standardAirComposition.getComponents(); + auto const& fractions = ::corsika::standardAirComposition.getFractions(); + if (::corsika::standardAirComposition.getSize() != 3) { + throw std::runtime_error{"CONEXhybrid only usable with standard 3-component air"}; + } + + std::transform(components.cbegin(), components.cend(), ::conex::cxair_.aira.begin(), + get_nucleus_A); + std::transform(components.cbegin(), components.cend(), ::conex::cxair_.airz.begin(), + get_nucleus_Z); + std::copy(fractions.cbegin(), fractions.cend(), ::conex::cxair_.airw.begin()); + + ::conex::cxair_.airava = + std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(), + ::conex::cxair_.aira.cbegin(), 0.); + ::conex::cxair_.airavz = + std::inner_product(::conex::cxair_.airw.cbegin(), ::conex::cxair_.airw.cend(), + ::conex::cxair_.airz.cbegin(), 0.); + + // this is the CONEX default but actually unused there + ::conex::cxair_.airi = {82.0e-09, 95.0e-09, 188.e-09}; + int randomSeeds[3] = {1234, 0, 0}; // SEEDS ARE NOT USED. All random numbers are obtained from // the CORSIKA 8 stream "conex" and "epos"! diff --git a/corsika/media/CORSIKA7Atmospheres.hpp b/corsika/media/CORSIKA7Atmospheres.hpp index 0eb0b6a7a38d3dc8b4a9de238fa9c47f6555c5d8..cba81771f84ae43f7d4910f50150acbafa72f172 100644 --- a/corsika/media/CORSIKA7Atmospheres.hpp +++ b/corsika/media/CORSIKA7Atmospheres.hpp @@ -11,6 +11,7 @@ #include <corsika/media/IRefractiveIndexModel.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> #include <corsika/framework/utility/ImplementsMixin.hpp> +#include <corsika/media/NuclearComposition.hpp> // for detail namespace, NoExtraModelInner, NoExtraModel and traits #include <corsika/detail/media/LayeredSphericalAtmosphereBuilder.hpp> @@ -204,6 +205,10 @@ namespace corsika { void create_5layer_atmosphere(TEnvironment& env, AtmosphereId const atmId, Point const& center, TArgs... args); + //! The standard/default air composition with fraction values based on CORSIKA 7 + static inline NuclearComposition const standardAirComposition{ + {Code::Nitrogen, Code::Oxygen, Code::Argon}, {0.78479, .21052, 0.00469}}; + } // namespace corsika -#include <corsika/detail/media/CORSIKA7Atmospheres.inl> \ No newline at end of file +#include <corsika/detail/media/CORSIKA7Atmospheres.inl> diff --git a/corsika/media/NuclearComposition.hpp b/corsika/media/NuclearComposition.hpp index 184853ef2cdc620c821b516c4760144d503ffc63..402969edb42b36705815548953546a49c8381e1d 100644 --- a/corsika/media/NuclearComposition.hpp +++ b/corsika/media/NuclearComposition.hpp @@ -51,7 +51,7 @@ namespace corsika { * @retval returns the vector with weighted return types of func. */ template <typename TFunction> - auto getWeighted(TFunction const& func) const; + auto getWeighted(TFunction func) const; /** * Sum all all relative composition weighted by func(element) @@ -65,8 +65,7 @@ namespace corsika { * @retval returns the weighted sum with the type defined by the return type of func. */ template <typename TFunction> - auto getWeightedSum(TFunction const& func) const - -> decltype(func(std::declval<Code>())); + auto getWeightedSum(TFunction func) const -> decltype(func(std::declval<Code>())); /** * Number of elements in the composition array diff --git a/corsika/modules/conex/CONEX_f.hpp b/corsika/modules/conex/CONEX_f.hpp index db80137ad9537ca7e1f01ce7240cd4c8473955d1..3339903ce96f8642a5ea7852869f48c56c0b1fbc 100644 --- a/corsika/modules/conex/CONEX_f.hpp +++ b/corsika/modules/conex/CONEX_f.hpp @@ -36,13 +36,18 @@ namespace conex { extern double double_rndm_interface(); - extern "C" {} - // the CONEX fortran interface extern "C" { extern struct { std::array<double, 16> dptl; } cxoptl_; + //! common block for atmosphere composition + extern struct { + std::array<double, 3> airz, aira, airw; //!< nuclear Z, A, composition fraction + double airavz, airava; //!< average Z, A + std::array<double, 3> airi; //!< ionization potential, not used in cxroot + } cxair_; + void cegs4_(int&, int&); void initconex_(int&, int*, int&, int&, diff --git a/modules/conex b/modules/conex index 344c875bed534590eb86aa5e12b0382a9868fc64..c26ca0fada6ef182147d014bb7d413814f9cf202 160000 --- a/modules/conex +++ b/modules/conex @@ -1 +1 @@ -Subproject commit 344c875bed534590eb86aa5e12b0382a9868fc64 +Subproject commit c26ca0fada6ef182147d014bb7d413814f9cf202