diff --git a/Environment/LinearIntegrator.h b/Environment/LinearIntegrator.h index 1bea44a296a0c1a4c77190676219c755845681e8..d1ee5d5ef065d5d6302f59ade3f0d0eb7d78eb7d 100644 --- a/Environment/LinearIntegrator.h +++ b/Environment/LinearIntegrator.h @@ -24,10 +24,9 @@ namespace corsika::environment { auto IntegrateGrammage( corsika::geometry::Trajectory<corsika::geometry::Line> const& line, corsika::units::si::LengthType length) const { - auto const c0 = GetImplementation().fRho(line.GetPosition(0)); - auto const c1 = GetImplementation().fRho.Derivative<1>(line.GetPosition(0), - line.NormalizedDirection()); - + auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0)); + auto const c1 = GetImplementation().fRho.FirstDerivative( + line.GetPosition(0), line.NormalizedDirection()); return (c0 + 0.5 * c1 * length) * length; } @@ -35,8 +34,8 @@ namespace corsika::environment { corsika::geometry::Trajectory<corsika::geometry::Line> const& line, corsika::units::si::GrammageType grammage) const { auto const c0 = GetImplementation().fRho(line.GetPosition(0)); - auto const c1 = GetImplementation().fRho.Derivative<1>(line.GetPosition(0), - line.NormalizedDirection()); + auto const c1 = GetImplementation().fRho.FirstDerivative( + line.GetPosition(0), line.NormalizedDirection()); return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0; } diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc index 01edec05fb0844b908f6967ab8cd0c43bd5c6bf0..3358cfa198a696a6e9313528de2cf04ffffacd50 100644 --- a/Environment/testEnvironment.cc +++ b/Environment/testEnvironment.cc @@ -47,6 +47,10 @@ struct Exponential { return v.GetComponents()[0] * (*this)(p) / corsika::units::si::detail::static_pow<N>(1_m); } + + auto FirstDerivative(Point const& p, Vector<dimensionless_d> const& v) const { + return Derivative<1>(p, v); + } }; TEST_CASE("DensityFunction") { @@ -57,7 +61,7 @@ TEST_CASE("DensityFunction") { Vector direction(cs, QuantityVector<dimensionless_d>(1, 0, 0)); Line line(origin, Vector<SpeedType::dimension_type>( - cs, {200_m / second, 0_m / second, 0_m / second})); + cs, {20_m / second, 0_m / second, 0_m / second})); auto const tEnd = 5_s; Trajectory<Line> const trajectory(line, tEnd); @@ -69,7 +73,13 @@ TEST_CASE("DensityFunction") { DensityFunction const rho(e); REQUIRE(rho.EvaluateAt(origin) == e(origin)); - auto const exactGrammage = [](auto x) { return rho0 * exp(x / 1_m); }; - REQUIRE(rho.IntegrateGrammage(trajectory, 5_cm) / exactGrammage(5_cm) == - Approx(1).epsilon(1e-2)); + 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)); }