Merge branch 'rprechelt-earth-radius' into 'master'

Add support for different Earth radii approximations into LayeredSphericalAtmosphereBuilder

See merge request !209

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));

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;
 

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

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);
 

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