diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index b04a58bdce307668b346ab1957b8e8b6325216a1..ef1d24035859fc6668d902ae89ee7aebc994a95b 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -6,7 +6,7 @@ set (
HomogeneousMedium.h
InhomogeneousMedium.h
HomogeneousMedium.h
- LinearIntegrator.h
+ LinearApproximationIntegrator.h
DensityFunction.h
Environment.h
)
diff --git a/Environment/DensityFunction.h b/Environment/DensityFunction.h
index 943410afc2a40226de94f9efac61c00af86ab225..32ee30adedfcbe9fccf1ebea2f09c8b9877f6015 100644
--- a/Environment/DensityFunction.h
+++ b/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
diff --git a/Environment/InhomogeneousMedium.h b/Environment/InhomogeneousMedium.h
index c506fb176d2d10370330c7b4d4d16eb45b304f95..feb6497cfcd82a28b8eafa183c068d106e187520 100644
--- a/Environment/InhomogeneousMedium.h
+++ b/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(
diff --git a/Environment/LinearIntegrator.h b/Environment/LinearApproximationIntegrator.h
similarity index 85%
rename from Environment/LinearIntegrator.h
rename to Environment/LinearApproximationIntegrator.h
index 00f5904655cd1cd4cfe6fe66add749178716ea66..910d2397189e7aa9b1456e00df94f72cfd49ff4f 100644
--- a/Environment/LinearIntegrator.h
+++ b/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
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 7f39ca1043bcb6df61800d7b0b041df95844335f..5bb18bfe09e48d4865f5aafcd8ad4aa01f4afeb8 100644
--- a/Environment/testEnvironment.cc
+++ b/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)));
+ }
}