more complete test

Environment/CMakeLists.txt

@@ -6,7 +6,7 @@ set (
   HomogeneousMedium.h
   InhomogeneousMedium.h
   HomogeneousMedium.h
-  LinearIntegrator.h
+  LinearApproximationIntegrator.h
   DensityFunction.h
   Environment.h
   )

Environment/DensityFunction.h

@@ -11,7 +11,7 @@
 #ifndef _include_environment_DensityFunction_h_
 #define _include_environment_DensityFunction_h_
 
-#include <corsika/environment/LinearIntegrator.h>
+#include <corsika/environment/LinearApproximationIntegrator.h>
 #include <corsika/geometry/Line.h>
 #include <corsika/geometry/Point.h>
 #include <corsika/geometry/Trajectory.h>
@@ -19,10 +19,10 @@
 namespace corsika::environment {
 
   template <class TDerivableRho,
-            template <typename> class TApproximator = LinearApproximator>
+            template <typename> class TIntegrator = LinearApproximationIntegrator>
   class DensityFunction
-      : public TApproximator<DensityFunction<TDerivableRho, TApproximator>> {
-    friend class TApproximator<DensityFunction<TDerivableRho, TApproximator>>;
+      : public TIntegrator<DensityFunction<TDerivableRho, TIntegrator>> {
+    friend class TIntegrator<DensityFunction<TDerivableRho, TIntegrator>>;
 
     TDerivableRho fRho; //!< functor for density
 

Environment/InhomogeneousMedium.h

@@ -20,8 +20,8 @@
 #include <corsika/units/PhysicalUnits.h>
 
 /**
- * An inhomogeneous medium. The mass density distribution TDensityFunction must be a
- * \f$C^1\f$-function.
+ * A general inhomogeneous medium. The mass density distribution TDensityFunction must be
+ * a \f$C^2\f$-function.
  */
 
 namespace corsika::environment {
@@ -46,7 +46,7 @@ namespace corsika::environment {
     corsika::units::si::GrammageType IntegratedGrammage(
         corsika::geometry::Trajectory<corsika::geometry::Line> const& pLine,
         corsika::units::si::LengthType pTo) const override {
-      return fDensityFunction.IntegratedGrammage(pLine, pTo);
+      return fDensityFunction.IntegrateGrammage(pLine, pTo);
     }
 
     corsika::units::si::LengthType ArclengthFromGrammage(

Environment/LinearIntegrator.h → Environment/LinearApproximationIntegrator.h

@@ -8,8 +8,8 @@
  * the license.
  */
 
-#ifndef _include_environment_LinearIntegrator_h_
-#define _include_environment_LinearIntegrator_h_
+#ifndef _include_environment_LinearApproximationIntegrator_h_
+#define _include_environment_LinearApproximationIntegrator_h_
 
 #include <limits>
 
@@ -19,7 +19,7 @@
 
 namespace corsika::environment {
   template <class TDerived>
-  class LinearApproximator {
+  class LinearApproximationIntegrator {
     auto const& GetImplementation() const { return *static_cast<TDerived const*>(this); }
 
   public:
@@ -43,9 +43,9 @@ namespace corsika::environment {
     }
 
     auto MaximumLength(corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
-                       double relError) const {
+                       [[maybe_unused]] double relError) const {
       using namespace corsika::units::si;
-      auto const c1 = GetImplementation().fRho.SecondDerivative(
+      [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
           line.GetPosition(0), line.NormalizedDirection());
 
       // todo: provide a real, working implementation

Environment/testEnvironment.cc

@@ -16,6 +16,7 @@
 #include <corsika/environment/HomogeneousMedium.h>
 #include <corsika/environment/IMediumModel.h>
 #include <corsika/environment/InhomogeneousMedium.h>
+#include <corsika/environment/LinearApproximationIntegrator.h>
 #include <corsika/environment/NuclearComposition.h>
 #include <corsika/environment/VolumeTreeNode.h>
 #include <corsika/geometry/Line.h>
@@ -27,12 +28,12 @@
 
 using namespace corsika::geometry;
 using namespace corsika::environment;
+using namespace corsika::particles;
 using namespace corsika::units::si;
 
 TEST_CASE("HomogeneousMedium") {
-  NuclearComposition const protonComposition(
-      std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
-      std::vector<float>{1.f});
+  NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+                                             std::vector<float>{1.f});
   HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
 }
 
@@ -52,10 +53,15 @@ struct Exponential {
   auto FirstDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
     return Derivative<1>(p, v);
   }
+
+  auto SecondDerivative(Point const& p, Vector<dimensionless_d> const& v) const {
+    return Derivative<2>(p, v);
+  }
 };
 
-TEST_CASE("DensityFunction") {
-  CoordinateSystem& cs = RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+TEST_CASE("InhomogeneousMedium") {
+  CoordinateSystem const& cs =
+      RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
   Point const origin(cs, {0_m, 0_m, 0_m});
 
@@ -68,19 +74,34 @@ TEST_CASE("DensityFunction") {
   Trajectory<Line> const trajectory(line, tEnd);
 
   Exponential const e;
-  REQUIRE(e.Derivative<1>(origin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
-          Approx(1));
+  DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
 
-  DensityFunction const rho(e);
-  REQUIRE(rho.EvaluateAt(origin) == e(origin));
+  SECTION("DensityFunction") {
+    REQUIRE(e.Derivative<1>(origin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+            Approx(1));
+    REQUIRE(rho.EvaluateAt(origin) == e(origin));
+  }
 
   auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
   auto const exactLength = [](auto X) { return 1_m * log(1 + X / (rho0 * 1_m)); };
 
   auto constexpr l = 15_cm;
-  CHECK(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
-        Approx(1).epsilon(1e-2));
-  CHECK(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
-            exactLength(exactGrammage(l)) ==
-        Approx(1).epsilon(1e-2));
+
+  NuclearComposition const composition{{Code::Proton}, {1.f}};
+  InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
+
+  SECTION("Integration") {
+    REQUIRE(rho.IntegrateGrammage(trajectory, l) / exactGrammage(l) ==
+            Approx(1).epsilon(1e-2));
+    REQUIRE(rho.ArclengthFromGrammage(trajectory, exactGrammage(l)) /
+                exactLength(exactGrammage(l)) ==
+            Approx(1).epsilon(1e-2));
+    REQUIRE(rho.MaximumLength(trajectory, 1e-2) >
+            l); // todo: write reasonable test when implementation is working
+
+    REQUIRE(rho.IntegrateGrammage(trajectory, l) ==
+            inhMedium.IntegratedGrammage(trajectory, l));
+    REQUIRE(rho.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
+            inhMedium.ArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
+  }
 }