Merge branch '151-exponential-flat-atmosphere' into 'master'

Resolve "Exponential flat atmosphere"

Closes #151

See merge request !92

Environment/CMakeLists.txt

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

Environment/FlatExponential.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.
+ */
+
+#ifndef _include_Environment_FlatExponential_h_
+#define _include_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>
+
+#include <cassert>
+
+/**
+ *
+ */
+
+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;
+    geometry::Vector<units::si::dimensionless_d> const fAxis;
+    geometry::Point const fP0;
+
+  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) {}
+
+    corsika::units::si::MassDensityType GetMassDensity(
+        corsika::geometry::Point const& p) const override {
+      return fRho0 * exp(fInvLambda * (p - 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);
+      }
+    }
+
+    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 {
+        return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+      }
+    }
+  };
+
+} // namespace corsika::environment
+#endif

Environment/testEnvironment.cc

@@ -13,6 +13,7 @@
                           // cpp file
 
 #include <corsika/environment/DensityFunction.h>
+#include <corsika/environment/FlatExponential.h>
 #include <corsika/environment/HomogeneousMedium.h>
 #include <corsika/environment/IMediumModel.h>
 #include <corsika/environment/InhomogeneousMedium.h>
@@ -30,6 +31,12 @@ using namespace corsika::geometry;
 using namespace corsika::environment;
 using namespace corsika::particles;
 using namespace corsika::units::si;
+using namespace corsika;
+
+CoordinateSystem const& gCS =
+    RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+Point const gOrigin(gCS, {0_m, 0_m, 0_m});
 
 TEST_CASE("HomogeneousMedium") {
   NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
@@ -37,6 +44,51 @@ TEST_CASE("HomogeneousMedium") {
   HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
 }
 
+TEST_CASE("FlatExponential") {
+  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+                                             std::vector<float>{1.f});
+
+  Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+  LengthType const lambda = 3_m;
+  auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+  FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
+                                             protonComposition);
+  auto const tEnd = 5_s;
+
+  SECTION("horizontal") {
+    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+                                 gCS, {20_cm / second, 0_m / second, 0_m / second}));
+    Trajectory<Line> const trajectory(line, tEnd);
+
+    REQUIRE((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
+    REQUIRE((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
+  }
+
+  SECTION("vertical") {
+    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+                                 gCS, {0_m / second, 0_m / second, 5_m / second}));
+    Trajectory<Line> const trajectory(line, tEnd);
+    LengthType const length = 2 * lambda;
+    GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
+
+    REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
+    REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
+  }
+
+  SECTION("vertical") {
+    Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+                                 gCS, {0_m / second, 5_m / second, 5_m / second}));
+    Trajectory<Line> const trajectory(line, tEnd);
+    double const cosTheta = M_SQRT1_2;
+    LengthType const length = 2 * lambda;
+    GrammageType const exact =
+        rho0 * lambda * (exp(cosTheta * length / lambda) - 1) / cosTheta;
+
+    REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
+    REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
+  }
+}
+
 auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
 
 struct Exponential {
@@ -60,15 +112,10 @@ struct Exponential {
 };
 
 TEST_CASE("InhomogeneousMedium") {
-  CoordinateSystem const& cs =
-      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  Point const origin(cs, {0_m, 0_m, 0_m});
-
-  Vector direction(cs, QuantityVector<dimensionless_d>(1, 0, 0));
+  Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
 
-  Line line(origin, Vector<SpeedType::dimension_type>(
-                        cs, {20_m / second, 0_m / second, 0_m / second}));
+  Line line(gOrigin, Vector<SpeedType::dimension_type>(
+                         gCS, {20_m / second, 0_m / second, 0_m / second}));
 
   auto const tEnd = 5_s;
   Trajectory<Line> const trajectory(line, tEnd);
@@ -77,9 +124,9 @@ TEST_CASE("InhomogeneousMedium") {
   DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
 
   SECTION("DensityFunction") {
-    REQUIRE(e.Derivative<1>(origin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+    REQUIRE(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
             Approx(1));
-    REQUIRE(rho.EvaluateAt(origin) == e(origin));
+    REQUIRE(rho.EvaluateAt(gOrigin) == e(gOrigin));
   }
 
   auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };

Framework/Geometry/Vector.h

@@ -125,17 +125,9 @@ namespace corsika::geometry {
 
     template <typename ScalarDim>
     auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
-      using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
+      using ProdDim = phys::units::detail::product_d<dim, ScalarDim>;
 
-      if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
-                                                               // not a "Quantity" anymore
-      {
-        return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs,
-                                                    BaseVector<dim>::qVector * p);
-      } else {
-        return Vector<typename ProdQuantity::dimension_type>(
-            *BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
-      }
+      return Vector<ProdDim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
     }
 
     template <typename ScalarDim>
@@ -173,16 +165,8 @@ namespace corsika::geometry {
       auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
       auto const bareResult = c1.cross(c2);
 
-      using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
-
-      if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
-                                                               // not a "Quantity" anymore
-      {
-        return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
-      } else {
-        return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs,
-                                                             bareResult);
-      }
+      using ProdDim = phys::units::detail::product_d<dim, dim2>;
+      return Vector<ProdDim>(*BaseVector<dim>::cs, bareResult);
     }
 
     template <typename dim2>
@@ -191,9 +175,10 @@ namespace corsika::geometry {
       auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
       auto const bareResult = c1.dot(c2);
 
-      using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+      using ProdDim = phys::units::detail::product_d<dim, dim2>;
 
-      return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
+      return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
+                                                    bareResult);
     }
   };
 

ThirdParty/phys/units/quantity.hpp

@@ -144,9 +144,18 @@ namespace phys {
           d8 = DX::dim8 + DY::dim8
         };
 
-        typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
       };
 
+      /**
+       * convenience type for product dimension
+       */
+
+      template <typename DX, typename DY>
+      using product_d = typename product<DX, DY, Rep>::dim;
+
       template <typename DX, typename DY, typename X, typename Y>
       using Product = typename product<DX, DY, PromoteMul<X, Y>>::type;
 
@@ -166,9 +175,18 @@ namespace phys {
           d8 = DX::dim8 - DY::dim8
         };
 
-        typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
       };
 
+      /**
+       * convenience type for quotient dimension
+       */
+
+      template <typename DX, typename DY>
+      using quotient_d = typename quotient<DX, DY, Rep>::dim;
+
       template <typename DX, typename DY, typename X, typename Y>
       using Quotient = typename quotient<DX, DY, PromoteMul<X, Y>>::type;
 
@@ -188,9 +206,18 @@ namespace phys {
           d8 = -D::dim8
         };
 
-        typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
       };
 
+      /**
+       * convenience type for reciprocal dimension
+       */
+
+      template <typename DX, typename DY>
+      using reciprocal_d = typename reciprocal<DX, Rep>::dim;
+
       template <typename D, typename X, typename Y>
       using Reciprocal = typename reciprocal<D, PromoteMul<X, Y>>::type;
 
@@ -210,11 +237,20 @@ namespace phys {
           d8 = N * D::dim8
         };
 
-        typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
       };
 
+      /**
+       * convenience type for power dimension
+       */
+
+      template <typename DX, int N>
+      using power_d = typename power<DX, N, Rep>::dim;
+
       template <typename D, int N, typename T>
-      using Power = typename detail::power<D, N, T>::type;
+      using Power = typename power<D, N, T>::type;
 
       /**
        * root type generator.
@@ -238,9 +274,17 @@ namespace phys {
           d8 = D::dim8 / N
         };
 
-        typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+        typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+        typedef Collapse<dim, T> type;
       };
 
+      /**
+       * convenience type for root dimension
+       */
+      template <typename D, int N>
+      using root_d = typename root<D, N, Rep>::dim;
+
       template <typename D, int N, typename T>
       using Root = typename detail::root<D, N, T>::type;