diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc index 36fa8d964dbfe6ceac77b373abc632c0299b5c98..a81414dc1e373aa21cea0fb4fa1e44a94b9b8336 100644 --- a/Documentation/Examples/vertical_EAS.cc +++ b/Documentation/Examples/vertical_EAS.cc @@ -117,8 +117,8 @@ int main(int argc, char** argv) { cout << "input momentum: " << plab.GetComponents() / 1_GeV << ", norm = " << plab.norm() << endl; - auto const observationHeight = 1.4_km + builder.earthRadius; - auto const injectionHeight = 112.75_km + builder.earthRadius; + auto const observationHeight = 1.4_km + builder.getEarthRadius(); + auto const injectionHeight = 112.75_km + builder.getEarthRadius(); auto const t = -observationHeight * cos(thetaRad) + sqrt(-si::detail::static_pow<2>(sin(thetaRad) * observationHeight) + si::detail::static_pow<2>(injectionHeight)); diff --git a/Environment/LayeredSphericalAtmosphereBuilder.cc b/Environment/LayeredSphericalAtmosphereBuilder.cc index 2366266e278606a684f935acf89f33ee73d49252..86ee4dfc450796da863b7b7a72c53e45878b04dd 100644 --- a/Environment/LayeredSphericalAtmosphereBuilder.cc +++ b/Environment/LayeredSphericalAtmosphereBuilder.cc @@ -29,7 +29,7 @@ void LayeredSphericalAtmosphereBuilder::setNuclearComposition( void LayeredSphericalAtmosphereBuilder::addExponentialLayer( units::si::GrammageType b, units::si::LengthType c, units::si::LengthType upperBoundary) { - auto const radius = seaLevel_ + upperBoundary; + auto const radius = earthRadius_ + upperBoundary; checkRadius(radius); previousRadius_ = radius; @@ -40,7 +40,7 @@ void LayeredSphericalAtmosphereBuilder::addExponentialLayer( std::cout << "rho0 = " << rho0 << ", c = " << c << std::endl; node->SetModelProperties<SlidingPlanarExponential<IMediumModel>>( - center_, rho0, -c, *composition_, seaLevel_); + center_, rho0, -c, *composition_, earthRadius_); layers_.push(std::move(node)); } @@ -49,7 +49,7 @@ void LayeredSphericalAtmosphereBuilder::addLinearLayer( units::si::LengthType c, units::si::LengthType upperBoundary) { using namespace units::si; - auto const radius = seaLevel_ + upperBoundary; + auto const radius = earthRadius_ + upperBoundary; checkRadius(radius); previousRadius_ = radius; diff --git a/Environment/LayeredSphericalAtmosphereBuilder.h b/Environment/LayeredSphericalAtmosphereBuilder.h index ab2b8a25d95887d198dc568536fec81c3d40b461..f051485ca7de45460003628be03052a988995466 100644 --- a/Environment/LayeredSphericalAtmosphereBuilder.h +++ b/Environment/LayeredSphericalAtmosphereBuilder.h @@ -13,6 +13,7 @@ #include <corsika/environment/NuclearComposition.h> #include <corsika/environment/VolumeTreeNode.h> #include <corsika/geometry/Point.h> +#include <corsika/units/PhysicalConstants.h> #include <corsika/units/PhysicalUnits.h> #include <memory> @@ -24,7 +25,7 @@ namespace corsika::environment { std::unique_ptr<NuclearComposition> composition_; geometry::Point center_; units::si::LengthType previousRadius_{units::si::LengthType::zero()}; - units::si::LengthType seaLevel_; + units::si::LengthType earthRadius_; std::stack<VolumeTreeNode<environment::IMediumModel>::VTNUPtr> layers_; // innermost layer first @@ -32,12 +33,11 @@ namespace corsika::environment { void checkRadius(units::si::LengthType) const; public: - static auto constexpr earthRadius = 6'371'000 * units::si::meter; - - LayeredSphericalAtmosphereBuilder(corsika::geometry::Point center, - units::si::LengthType seaLevel = earthRadius) + LayeredSphericalAtmosphereBuilder( + corsika::geometry::Point center, + units::si::LengthType earthRadius = units::constants::EarthRadius::Mean) : center_(center) - , seaLevel_(seaLevel) {} + , earthRadius_(earthRadius) {} void setNuclearComposition(NuclearComposition); @@ -50,6 +50,11 @@ namespace corsika::environment { void assemble(Environment<IMediumModel>&); Environment<IMediumModel> assemble(); + + /** + * Get the current Earth radius. + */ + units::si::LengthType getEarthRadius() const { return earthRadius_; }; }; } // namespace corsika::environment diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc index 31955e15ac8618e431315fb56685429f6a07f046..ed2bf9a1a7b232e95c1e04096abbcd3062c3839c 100644 --- a/Environment/testEnvironment.cc +++ b/Environment/testEnvironment.cc @@ -218,7 +218,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") { REQUIRE(builder.size() == 0); - auto constexpr R = LayeredSphericalAtmosphereBuilder::earthRadius; + auto const R = builder.getEarthRadius(); REQUIRE(univ->GetChildNodes().size() == 1); diff --git a/Framework/Units/PhysicalConstants.h b/Framework/Units/PhysicalConstants.h index eacccfe165b779f810156cdf3c4d30df35d99dba..da86bbd1edc3564badee1dedfd0bf7ead7d2ad4f 100644 --- a/Framework/Units/PhysicalConstants.h +++ b/Framework/Units/PhysicalConstants.h @@ -53,6 +53,16 @@ namespace corsika::units::constants { // molar gas constant auto constexpr R = Rep(8.314'459'8) * joule / (mole * kelvin); + /** + * A namespace containing various Earth radii. + */ + namespace EarthRadius { + static constexpr auto Mean{6'371'000 * meter}; + static constexpr auto Eqautorial{6'378'137 * meter}; + static constexpr auto Polar{6'356'752 * meter}; + static constexpr auto PolarCurvature{6'399'593 * meter}; + } // namespace EarthRadius + // etc. } // namespace corsika::units::constants