Merge branch 'sliding_planar_exp_atmo'

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 <limits>
+
+namespace corsika::environment {
+
+  /**
+   * This class provides the grammage/length conversion functionality for
+   * (locally) flat exponential atmospheres.
+   */
+  template <class TDerived>
+  class BaseExponential {
+  protected:
+    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); }
+
+    // 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
+     * \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.
+     * 
+     * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
+     * \f[
+     *   X = \varrho_0 l;
+     * \f]
+     */
+    // clang-format on
+    units::si::GrammageType IntegratedGrammage(
+        geometry::Trajectory<geometry::Line> const& vLine, units::si::LengthType vL,
+        geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
+      auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
+      auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+
+      if (uDotA == 0) {
+        return vL * rhoStart;
+      } else {
+        return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1);
+      }
+    }
+
+    // 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
+     * \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,}
+     *   \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:
+     * \f[
+     *   l =  \frac{X}{\varrho_0}
+     * \f]
+     */
+    // clang-format on
+    units::si::LengthType ArclengthFromGrammage(
+        geometry::Trajectory<geometry::Line> const& vLine,
+        units::si::GrammageType vGrammage,
+        geometry::Vector<units::si::dimensionless_d> const& vAxis) const {
+      auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
+      auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+
+      if (uDotA == 0) {
+        return vGrammage / rhoStart;
+      } else {
+        auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1;
+        if (logArg > 0) {
+          return fLambda / uDotA * log(logArg);
+        } else {
+          return std::numeric_limits<typename decltype(
+                     vGrammage)::value_type>::infinity() *
+                 units::si::meter;
+        }
+      }
+    }
+
+  public:
+    BaseExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
+                    units::si::LengthType vLambda)
+        : fRho0(vRho)
+        , fLambda(vLambda)
+        , fInvLambda(1 / vLambda)
+        , fP0(vP0) {}
+  };
+
+} // namespace corsika::environment
+#endif

Environment/CMakeLists.txt

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

Environment/FlatAtmosphere/FlatAtmosphere.h

-
-/*
- * (c) Copyright 2018 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.
- */

Environment/FlatExponential.h

@@ -12,82 +12,60 @@
 #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>
-#include <limits>
-
-/**
- *
- */
-
 namespace corsika::environment {
 
+  //clang-format off
+  /**
+   * flat exponential density distribution with
+   * \f[
+   *  \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - 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$.
+   */
+  //clang-format on
   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)
-        , fAxis(axis)
-        , fP0(p0) {}
+    FlatExponential(geometry::Point const& vP0,
+                    geometry::Vector<units::si::dimensionless_d> const& vAxis,
+                    units::si::MassDensityType vRho, units::si::LengthType vLambda,
+                    NuclearComposition vNuclComp)
+        : Base(vP0, vRho, vLambda)
+        , fAxis(vAxis)
+        , fNuclComp(vNuclComp) {}
 
-    corsika::units::si::MassDensityType GetMassDensity(
-        corsika::geometry::Point const& p) const override {
-      return fRho0 * exp(fInvLambda * (p - fP0).dot(fAxis));
+    units::si::MassDensityType GetMassDensity(geometry::Point const& vP) const override {
+      return Base::fRho0 * exp(Base::fInvLambda * (vP - 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();
+    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
 
-      if (vDotA == 0) {
-        return pTo * GetMassDensity(line.GetR0());
-      } else {
-        return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
-               (exp(vDotA * pTo / fLambda) - 1);
-      }
+    units::si::GrammageType IntegratedGrammage(
+        geometry::Trajectory<geometry::Line> const& vLine,
+        units::si::LengthType vTo) const override {
+      return Base::IntegratedGrammage(vLine, vTo, 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;
-        }
-      }
+    units::si::LengthType ArclengthFromGrammage(
+        geometry::Trajectory<geometry::Line> const& vLine,
+        units::si::GrammageType vGrammage) const override {
+      return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
     }
   };
-
 } // namespace corsika::environment
 #endif

Environment/SlidingPlanarExponential.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_SlidingPlanarExponential_h_
+#define _include_Environment_SlidingPlanarExponential_h_
+
+#include <corsika/environment/FlatExponential.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>
+
+namespace corsika::environment {
+
+  // clang-format off
+  /**
+   * The SlidingPlanarExponential models mass density as
+   * \f[
+   *   \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right).
+   * \f]
+   * For grammage/length conversion, the density distribution is approximated as
+   * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the axis pointing
+   * from \f$ p_0 \f$ to \f$ r_0 \f$.
+   */
+   //clang-format on
+   
+  template <class T>
+  class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
+                                   public T {
+    NuclearComposition const fNuclComp;
+
+    using Base = BaseExponential<SlidingPlanarExponential<T>>;
+
+  public:
+    SlidingPlanarExponential(geometry::Point const& vP0, units::si::MassDensityType vRho,
+                             units::si::LengthType vLambda, NuclearComposition vNuclComp)
+        : Base(vP0, vRho, vLambda)
+        , fNuclComp(vNuclComp) {}
+
+    units::si::MassDensityType GetMassDensity(
+        geometry::Point const& vP) const override {
+      auto const height = (vP - Base::fP0).norm();
+      return Base::fRho0 * exp(Base::fInvLambda * height);
+    }
+
+    NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+
+    units::si::GrammageType IntegratedGrammage(
+        geometry::Trajectory<geometry::Line> const& vLine,
+        units::si::LengthType vL) const override {
+      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
+      return Base::IntegratedGrammage(vLine, vL, axis);
+    }
+
+    units::si::LengthType ArclengthFromGrammage(
+        geometry::Trajectory<geometry::Line> const& vLine,
+        units::si::GrammageType vGrammage) const override {
+      auto const axis = (vLine.GetR0() - Base::fP0).normalized();
+      return Base::ArclengthFromGrammage(vLine, vGrammage, axis);
+    }
+  };
+
+} // namespace corsika::environment
+#endif

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 {