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..ba148f83dd6681a4ebac84b9a8fb8cec39c9c172 --- /dev/null +++ b/corsika/detail/media/CORSIKA7Atmospheres.inl @@ -0,0 +1,49 @@ +/* + * (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); + + // check if we want to also add the US standard refractivity + if constexpr (detail::implements_mixin_v<IRefractiveIndexModel, + TEnvironmentInterface>) { + // TODO: Add US Standard refractivity + } + + // 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/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/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/USStandardAtmosphere.hpp b/corsika/media/CORSIKA7Atmospheres.hpp similarity index 83% rename from corsika/media/USStandardAtmosphere.hpp rename to corsika/media/CORSIKA7Atmospheres.hpp index bc353c7952b184751ab3fd01d035896d5ef88d0d..0eb0b6a7a38d3dc8b4a9de238fa9c47f6555c5d8 100644 --- a/corsika/media/USStandardAtmosphere.hpp +++ b/corsika/media/CORSIKA7Atmospheres.hpp @@ -12,6 +12,9 @@ #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 { /** @@ -46,15 +49,36 @@ namespace corsika { LastAtmosphere }; - namespace detail { + 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, @@ -160,39 +184,26 @@ namespace corsika { {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 + } // namespace - template <typename TEnvironmentInterface, template <typename> typename TExtraEnv, + /** + * 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> - auto 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)...); - - // 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.7847, 1. - 0.7847}}); - - // add the standard atmosphere layers - auto const params = detail::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); - - // check if we want to also add the US standard refractivity - if constexpr (detail::implements_mixin_v<IRefractiveIndexModel, - TEnvironmentInterface>) { - // TODO: Add US Standard refractivity - } - - // and assemble the environment - builder.assemble(env); - }; + 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 3b794a34fd3ace8ed596f8743659a16933daf938..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 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/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/testEnvironment.cpp b/tests/media/testEnvironment.cpp index c9d6a673d3c01ca2d58da1ba3a05f9e11aea27df..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, @@ -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); @@ -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);