Add earthRadius as runtime constructor arg.

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

@@ -20,11 +20,21 @@
 
 namespace corsika::environment {
 
+  /**
+   * A namespace containing various Earth radii.
+   */
+  namespace EarthRadius {
+    static constexpr auto Mean{6'371'000 * units::si::meter};
+    static constexpr auto Eqautorial{6'378'137 * units::si::meter};
+    static constexpr auto Polar{6'356'752 * units::si::meter};
+    static constexpr auto PolarCurvature{6'399'593 * units::si::meter};
+  } // namespace EarthRadius
+
   class LayeredSphericalAtmosphereBuilder {
     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 +42,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 = EarthRadius::Mean)
         : center_(center)
-        , seaLevel_(seaLevel) {}
+        , earthRadius_(earthRadius) {}
 
     void setNuclearComposition(NuclearComposition);
 
@@ -50,6 +59,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);