first steps

Environment/BaseExponential.h

+/*
+ * (c) Copyright 2019 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_Environment_BaseExponential_h_
+#define _include_Environment_BaseExponential_h_
+
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <cassert>
+#include <limits>
+
+/**
+ *
+ */
+
+namespace corsika::environment {
+
+  template <class TDerived>
+  class BaseExponential {
+  protected:
+    corsika::units::si::MassDensityType const fRho0;
+    units::si::LengthType const fLambda;
+    units::si::InverseLengthType const fInvLambda;
+    geometry::Point const fP0;
+
+    auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
+
+    corsika::units::si::GrammageType IntegratedGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+        corsika::units::si::LengthType pTo,
+        geometry::Vector<units::si::dimensionless_d> const& axis) const {
+      auto const vDotA = line.NormalizedDirection().dot(axis).magnitude();
+
+      if (vDotA == 0) {
+        return pTo * GetImplementation().GetMassDensity(line.GetR0());
+      } else {
+        return GetImplementation().GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
+               (exp(vDotA * pTo / fLambda) - 1);
+      }
+    }
+
+    corsika::units::si::LengthType ArclengthFromGrammage(
+        corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+        corsika::units::si::GrammageType pGrammage,
+        geometry::Vector<units::si::dimensionless_d> const& axis) const {
+      auto const vDotA = line.NormalizedDirection().dot(axis).magnitude();
+
+      if (vDotA == 0) {
+        return pGrammage / GetImplementation().GetMassDensity(line.GetR0());
+      } else {
+        auto const logArg = pGrammage * vDotA / (fRho0 * fLambda) + 1;
+        if (logArg > 0) {
+          return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+        } else {
+          return std::numeric_limits<typename decltype(
+                     pGrammage)::value_type>::infinity() *
+                 corsika::units::si::meter;
+        }
+      }
+    }
+
+  public:
+    BaseExponential(geometry::Point const& p0, units::si::MassDensityType rho,
+                    units::si::LengthType lambda)
+        : fRho0(rho)
+        , fLambda(lambda)
+        , fInvLambda(1 / lambda)
+        , fP0(p0) {}
+  };
+
+} // namespace corsika::environment
+#endif

Environment/CMakeLists.txt

@@ -9,6 +9,7 @@ set (
   LinearApproximationIntegrator.h
   DensityFunction.h
   Environment.h
+  BaseExponential.h
   FlatExponential.h
   SlidingPlanarExponential.h
   )

Environment/FlatExponential.h

@@ -11,12 +11,12 @@
 #ifndef _include_Environment_FlatExponential_h_
 #define _include_Environment_FlatExponential_h_
 
+#include <corsika/environment/BaseExponential.h>
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/particles/ParticleProperties.h>
-#include <corsika/random/RNGManager.h>
 #include <corsika/units/PhysicalUnits.h>
 
 #include <cassert>
@@ -29,64 +29,39 @@
 namespace corsika::environment {
 
   template <class T>
-  class FlatExponential : public T {
-    corsika::units::si::MassDensityType const fRho0;
-    units::si::LengthType const fLambda;
-    units::si::InverseLengthType const fInvLambda;
-    NuclearComposition const fNuclComp;
+  class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
     geometry::Vector<units::si::dimensionless_d> const fAxis;
-    geometry::Point const fP0;
+    NuclearComposition const fNuclComp;
+
+    using Base = BaseExponential<FlatExponential<T>>;
 
   public:
     FlatExponential(geometry::Point const& p0,
                     geometry::Vector<units::si::dimensionless_d> const& axis,
                     units::si::MassDensityType rho, units::si::LengthType lambda,
-                    NuclearComposition pNuclComp)
-        : fRho0(rho)
-        , fLambda(lambda)
-        , fInvLambda(1 / lambda)
-        , fNuclComp(pNuclComp)
+                    NuclearComposition nuclComp)
+        : Base(p0, rho, lambda)
         , fAxis(axis)
-        , fP0(p0) {}
+        , fNuclComp(nuclComp) {}
 
     corsika::units::si::MassDensityType GetMassDensity(
         corsika::geometry::Point const& p) const override {
-      return fRho0 * exp(fInvLambda * (p - fP0).dot(fAxis));
+      return Base::fRho0 * exp(Base::fInvLambda * (p - Base::fP0).dot(fAxis));
     }
+
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
     corsika::units::si::GrammageType IntegratedGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
         corsika::units::si::LengthType pTo) const override {
-      auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
-
-      if (vDotA == 0) {
-        return pTo * GetMassDensity(line.GetR0());
-      } else {
-        return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
-               (exp(vDotA * pTo / fLambda) - 1);
-      }
+      return Base::IntegratedGrammage(line, pTo, fAxis);
     }
 
     corsika::units::si::LengthType ArclengthFromGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
         corsika::units::si::GrammageType pGrammage) const override {
-      auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
-
-      if (vDotA == 0) {
-        return pGrammage / GetMassDensity(line.GetR0());
-      } else {
-        auto const logArg = pGrammage * vDotA / (fRho0 * fLambda) + 1;
-        if (logArg > 0) {
-          return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
-        } else {
-          return std::numeric_limits<typename decltype(
-                     pGrammage)::value_type>::infinity() *
-                 corsika::units::si::meter;
-        }
-      }
+      return Base::ArclengthFromGrammage(line, pGrammage, fAxis);
     }
   };
-
 } // namespace corsika::environment
 #endif

Environment/SlidingPlanarExponential.h

@@ -11,15 +11,16 @@
 #ifndef _include_Environment_SlidingPlanarExponential_h_
 #define _include_Environment_SlidingPlanarExponential_h_
 
+#include <corsika/environment/FlatExponential.h>
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/geometry/Line.h>
-#include <corsika/environment/FlatExponential.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/random/RNGManager.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <fenv.h>
 #include <cassert>
 #include <limits>
 
@@ -30,30 +31,22 @@
 namespace corsika::environment {
 
   template <class T>
-  class SlidingPlanarExponential : public T {
-    corsika::units::si::MassDensityType const fRho0;
-    units::si::LengthType const fLambda;
-    units::si::InverseLengthType const fInvLambda;
+  class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
+                                   public T {
     NuclearComposition const fNuclComp;
-    geometry::Vector<units::si::dimensionless_d> const fAxis;
-    geometry::Point const fP0;
+
+    using Base = BaseExponential<SlidingPlanarExponential<T>>;
 
   public:
-    SlidingPlanarExponential(geometry::Point const& p0,
-                    geometry::Vector<units::si::dimensionless_d> const& axis,
-                    units::si::MassDensityType rho, units::si::LengthType lambda,
-                    NuclearComposition pNuclComp)
-        : fRho0(rho)
-        , fLambda(lambda)
-        , fInvLambda(1 / lambda)
-        , fNuclComp(pNuclComp)
-        , fAxis(axis)
-        , fP0(p0) {}
+    SlidingPlanarExponential(geometry::Point const& p0, units::si::MassDensityType rho,
+                             units::si::LengthType lambda, NuclearComposition nuclComp)
+        : Base(p0, rho, lambda)
+        , fNuclComp(nuclComp) {}
 
     corsika::units::si::MassDensityType GetMassDensity(
         corsika::geometry::Point const& p) const override {
-      auto const height = (p - fP0).norm();
-      return fRho0 * exp(fInvLambda * height);
+      auto const height = (p - Base::fP0).norm();
+      return Base::fRho0 * exp(Base::fInvLambda * height);
     }
 
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
@@ -61,36 +54,17 @@ namespace corsika::environment {
     corsika::units::si::GrammageType IntegratedGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
         corsika::units::si::LengthType pTo) const override {
-
-      auto const axis = (fP0 - line.GetR0()).normalized();
-      auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
-
-      if (vDotA == 0) {
-        return pTo * GetMassDensity(line.GetR0());
-      } else {
-        return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
-               (exp(vDotA * pTo / fLambda) - 1);
-      }
+      feenableexcept(FE_INVALID);
+      auto const axis = (Base::fP0 - line.GetR0()).normalized();
+      return Base::IntegratedGrammage(line, pTo, axis);
     }
 
     corsika::units::si::LengthType ArclengthFromGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
         corsika::units::si::GrammageType pGrammage) const override {
-      auto const axis = (fP0 - line.GetR0()).normalized();
-      auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
-
-      if (vDotA == 0) {
-        return pGrammage / GetMassDensity(line.GetR0());
-      } else {
-        auto const logArg = pGrammage * vDotA / (fRho0 * fLambda) + 1;
-        if (logArg > 0) {
-          return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
-        } else {
-          return std::numeric_limits<typename decltype(
-                     pGrammage)::value_type>::infinity() *
-                 corsika::units::si::meter;
-        }
-      }
+      feenableexcept(FE_INVALID);
+      auto const axis = (Base::fP0 - line.GetR0()).normalized();
+      return Base::ArclengthFromGrammage(line, pGrammage, axis);
     }
   };
 

Environment/testEnvironment.cc

@@ -19,6 +19,7 @@
 #include <corsika/environment/InhomogeneousMedium.h>
 #include <corsika/environment/LinearApproximationIntegrator.h>
 #include <corsika/environment/NuclearComposition.h>
+#include <corsika/environment/SlidingPlanarExponential.h>
 #include <corsika/environment/VolumeTreeNode.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/RootCoordinateSystem.h>
@@ -100,6 +101,41 @@ TEST_CASE("FlatExponential") {
   }
 }
 
+TEST_CASE("SlidingPlanarExponential") {
+  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+                                             std::vector<float>{1.f});
+
+  LengthType const lambda = 3_m;
+  auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+  auto const tEnd = 5_s;
+
+  SlidingPlanarExponential<IMediumModel> const medium(gOrigin, rho0, lambda,
+                                                      protonComposition);
+
+  SECTION("density") {
+    CHECK(medium.GetMassDensity({gCS, {0_m, 0_m, 3_m}}) /
+              medium.GetMassDensity({gCS, {0_m, 3_m, 0_m}}) ==
+          Approx(1));
+  }
+
+  SECTION("vertical") {
+    Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+    FlatExponential<IMediumModel> const flat(gOrigin, axis, rho0, lambda,
+                                             protonComposition);
+    Line const line({gCS, {0_m, 0_m, 1_m}},
+                    Vector<SpeedType::dimension_type>(
+                        gCS, {0_m / second, 0_m / second, 5_m / second}));
+    Trajectory<Line> const trajectory(line, tEnd);
+
+    CHECK(medium.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
+          flat.GetMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
+    CHECK(medium.IntegratedGrammage(trajectory, 2_m).magnitude() ==
+          flat.IntegratedGrammage(trajectory, 2_m).magnitude());
+    CHECK(medium.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
+          flat.ArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
+  }
+}
+
 auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
 
 struct Exponential {