diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl index ea0fe5b0e8bd96e0c9a68a0680bec2b1b0f7590c..1c5a98569ad4197dcff0c8423d85b5f82426e7eb 100644 --- a/corsika/detail/media/FlatExponential.inl +++ b/corsika/detail/media/FlatExponential.inl @@ -29,7 +29,7 @@ namespace corsika { template <typename T> inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const { return BaseExponential<FlatExponential<T>>::getMassDensity( - (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm()); + (point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).dot(axis_)); } template <typename T> diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp index 81c1ac75701baaf8c670ddb7525d3d41dfc2d1f5..50dfa4ecd36a28ab4b0d968d9b5fb8cd0e27a0ec 100644 --- a/tests/media/testEnvironment.cpp +++ b/tests/media/testEnvironment.cpp @@ -101,11 +101,19 @@ TEST_CASE("FlatExponential") { LengthType const length = 2_m; TimeType const tEnd = length / speed; + CHECK((medium.getMassDensity(gOrigin)) == rho0); CHECK(medium.getNuclearComposition().getFractions() == std::vector<double>{1.}); CHECK(medium.getNuclearComposition().getComponents() == std::vector<Code>{Code::Proton}); SECTION("horizontal") { + // Check that not moving along axis does not change density + CHECK(medium.getMassDensity(Point(gCS, 1_m, 0_m, 0_m)) == rho0); + CHECK(medium.getMassDensity(Point(gCS, 1_m, 0_m, 0_m)) == rho0); + CHECK(medium.getMassDensity(Point(gCS, -1_m, 0_m, 0_m)) == rho0); + CHECK(medium.getMassDensity(Point(gCS, 0_m, 1_m, 0_m)) == rho0); + CHECK(medium.getMassDensity(Point(gCS, 0_m, -1_m, 0_m)) == rho0); + Line const line(gOrigin, Vector<SpeedType::dimension_type>( gCS, {speed, 0_m / second, 0_m / second})); setup::Trajectory const trajectory = @@ -117,6 +125,10 @@ TEST_CASE("FlatExponential") { } SECTION("vertical") { + // Moving along axis does change density + CHECK(medium.getMassDensity(Point(gCS, 0_m, 0_m, 1_m)) > rho0); + CHECK(medium.getMassDensity(Point(gCS, 0_m, 0_m, -1_m)) < rho0); + Line const line(gOrigin, Vector<SpeedType::dimension_type>( gCS, {0_m / second, 0_m / second, speed})); setup::Trajectory const trajectory =