IAP GITLAB

Skip to content
Snippets Groups Projects
Commit b7deecfa authored by ralfulrich's avatar ralfulrich Committed by Ralf Ulrich
Browse files

some density modelling refactoring

parent 27604f72
No related branches found
No related tags found
1 merge request!375Resolve "Need martian atmosphere model for ICRC"
Hide whitespace changes
Inline Side-by-side
Showing
with 176 additions and 102 deletions
corsika/detail/media/BaseExponential.inl
+ 30
− 17
View file @ b7deecfa
......@@ -14,23 +14,43 @@
namespace corsika {
template <typename TDerived>
inline BaseExponential<TDerived>::BaseExponential(Point const& point,
LengthType const referenceHeight,
MassDensityType rho0,
LengthType lambda)
: rho0_(rho0)
, lambda_(lambda)
, invLambda_(1 / lambda)
, point_(point)
, referenceHeight_(referenceHeight) {}
template <typename TDerived>
inline auto const& BaseExponential<TDerived>::getImplementation() const {
return *static_cast<TDerived const*>(this);
}
template <typename TDerived>
inline MassDensityType BaseExponential<TDerived>::getMassDensity(
LengthType const height) const {
return rho0_ * exp(invLambda_ * (height - referenceHeight_));
}
template <typename TDerived>
inline GrammageType BaseExponential<TDerived>::getIntegratedGrammage(
BaseTrajectory const& traj, LengthType vL, DirectionVector const& axis) const {
if (vL == LengthType::zero()) { return GrammageType::zero(); }
BaseTrajectory const& traj, DirectionVector const& axis) const {
LengthType const length = traj.getLength();
if (length == LengthType::zero()) { return GrammageType::zero(); }
auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
// this corresponds to height:
double const uDotA = traj.getDirection(0).dot(axis).magnitude();
MassDensityType const rhoStart =
getImplementation().getMassDensity(traj.getPosition(0));
if (uDotA == 0) {
return vL * rhoStart;
return length * rhoStart;
} else {
return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1);
return rhoStart * (lambda_ / uDotA) * (exp(uDotA * length * invLambda_) - 1);
}
}
......@@ -38,8 +58,10 @@ namespace corsika {
inline LengthType BaseExponential<TDerived>::getArclengthFromGrammage(
BaseTrajectory const& traj, GrammageType grammage,
DirectionVector const& axis) const {
auto const uDotA = traj.getDirection(0).dot(axis).magnitude();
auto const rhoStart = getImplementation().getMassDensity(traj.getPosition(0));
// this corresponds to height:
double const uDotA = traj.getDirection(0).dot(axis).magnitude();
MassDensityType const rhoStart =
getImplementation().getMassDensity(traj.getPosition(0));
if (uDotA == 0) {
return grammage / rhoStart;
......@@ -54,13 +76,4 @@ namespace corsika {
}
}
template <typename TDerived>
inline BaseExponential<TDerived>::BaseExponential(Point const& point,
MassDensityType rho0,
LengthType lambda)
: rho0_(rho0)
, lambda_(lambda)
, invLambda_(1 / lambda)
, point_(point) {}
} // namespace corsika
corsika/detail/media/FlatExponential.inl
+ 6
− 8
View file @ b7deecfa
......@@ -18,19 +18,17 @@ namespace corsika {
template <typename T>
inline FlatExponential<T>::FlatExponential(Point const& point,
Vector<dimensionless_d> const& axis,
DirectionVector const& axis,
MassDensityType rho, LengthType lambda,
NuclearComposition const& nuclComp)
: BaseExponential<FlatExponential<T>>(point, rho, lambda)
: BaseExponential<FlatExponential<T>>(point, 0_m, rho, lambda)
, axis_(axis)
, nuclComp_(nuclComp) {}
template <typename T>
inline MassDensityType FlatExponential<T>::getMassDensity(Point const& point) const {
return BaseExponential<FlatExponential<T>>::getRho0() *
exp(BaseExponential<FlatExponential<T>>::getInvLambda() *
(point - BaseExponential<FlatExponential<T>>::getAnchorPoint())
.dot(axis_));
return BaseExponential<FlatExponential<T>>::getMassDensity(
(point - BaseExponential<FlatExponential<T>>::getAnchorPoint()).getNorm());
}
template <typename T>
......@@ -40,8 +38,8 @@ namespace corsika {
template <typename T>
inline GrammageType FlatExponential<T>::getIntegratedGrammage(
BaseTrajectory const& line, LengthType to) const {
return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, to, axis_);
BaseTrajectory const& line) const {
return BaseExponential<FlatExponential<T>>::getIntegratedGrammage(line, axis_);
}
template <typename T>
......
corsika/detail/media/HomogeneousMedium.inl
+ 3
− 3
View file @ b7deecfa
......@@ -32,9 +32,9 @@ namespace corsika {
}
template <typename T>
inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(BaseTrajectory const&,
LengthType to) const {
return to * density_;
inline GrammageType HomogeneousMedium<T>::getIntegratedGrammage(
BaseTrajectory const& track) const {
return track.getLength() * density_;
}
template <typename T>
......
corsika/detail/media/InhomogeneousMedium.inl
+ 2
− 2
View file @ b7deecfa
......@@ -36,8 +36,8 @@ namespace corsika {
template <typename T, typename TDensityFunction>
inline GrammageType InhomogeneousMedium<T, TDensityFunction>::getIntegratedGrammage(
BaseTrajectory const& line, LengthType to) const {
return densityFunction_.getIntegrateGrammage(line, to);
BaseTrajectory const& line) const {
return densityFunction_.getIntegrateGrammage(line);
}
template <typename T, typename TDensityFunction>
......
corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+ 35
− 0
View file @ b7deecfa
......@@ -13,6 +13,7 @@
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
//#include <corsika/media/SlidingPlanarTabular.hpp>
namespace corsika {
......@@ -103,6 +104,40 @@ namespace corsika {
layers_.push(std::move(node));
}
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
inline void
LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
addTabularLayer(GrammageType b,
std::function<MassDensityType(LengthType)> const& funcRho,
LengthType upperBoundary) {
auto const radius = earthRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<Sphere>(center_, radius));
/*
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
// helper lambda in which the last 5 arguments to make_shared<...> are bound
auto lastBound = [&](auto... argPack) {
return std::make_shared<
TMediumModelExtra<SlidingPlanarTabular<TMediumInterface>>>(
argPack..., center_, rho0, -c, *composition_, earthRadius_);
};
// now unpack the additional arguments
auto model = std::apply(lastBound, additionalModelArgs_);
node->setModelProperties(std::move(model));
} else {
node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
center_, rho0, -c, *composition_, earthRadius_);
}
*/
layers_.push(std::move(node));
}
template <typename TMediumInterface, template <typename> typename TMediumModelExtra,
typename... TModelArgs>
inline Environment<TMediumInterface> LayeredSphericalAtmosphereBuilder<
......
corsika/detail/media/LinearApproximationIntegrator.inl
+ 5
− 1
View file @ b7deecfa
......@@ -19,10 +19,14 @@ namespace corsika {
template <typename TDerived>
inline auto LinearApproximationIntegrator<TDerived>::getIntegrateGrammage(
BaseTrajectory const& line, LengthType length) const {
BaseTrajectory const& line) const {
LengthType const length = line.getLength();
auto const c0 = getImplementation().evaluateAt(line.getPosition(0));
auto const c1 = getImplementation().rho_.getFirstDerivative(line.getPosition(0),
line.getDirection(0));
CORSIKA_LOG_INFO("length={} c0={} c1={} pos={} dir={} return={}", length, c0, c1,
line.getPosition(0), line.getDirection(0),
(c0 + 0.5 * c1 * length) * length);
return (c0 + 0.5 * c1 * length) * length;
}
......
corsika/detail/media/SlidingPlanarExponential.inl
+ 6
− 8
View file @ b7deecfa
......@@ -16,19 +16,17 @@ namespace corsika {
inline SlidingPlanarExponential<T>::SlidingPlanarExponential(
Point const& p0, MassDensityType rho0, LengthType lambda,
NuclearComposition const& nuclComp, LengthType referenceHeight)
: BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
: BaseExponential<SlidingPlanarExponential<T>>(p0, referenceHeight, rho0, lambda)
, nuclComp_(nuclComp)
, referenceHeight_(referenceHeight) {}
template <typename T>
inline MassDensityType SlidingPlanarExponential<T>::getMassDensity(
Point const& point) const {
auto const height =
auto const heightFull =
(point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
.getNorm() -
referenceHeight_;
return BaseExponential<SlidingPlanarExponential<T>>::getRho0() *
exp(BaseExponential<SlidingPlanarExponential<T>>::getInvLambda() * height);
.getNorm();
return BaseExponential<SlidingPlanarExponential<T>>::getMassDensity(heightFull);
}
template <typename T>
......@@ -39,11 +37,11 @@ namespace corsika {
template <typename T>
inline GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
BaseTrajectory const& traj, LengthType l) const {
BaseTrajectory const& traj) const {
auto const axis = (traj.getPosition(0) -
BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
.normalized();
return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj, l,
return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj,
axis);
}
......
corsika/detail/modules/energy_loss/BetheBlochPDG.inl
+ 1
− 2
View file @ b7deecfa
......@@ -154,8 +154,7 @@ namespace corsika {
if (p.getChargeNumber() == 0) return ProcessReturn::Ok;
GrammageType const dX =
p.getNode()->getModelProperties().getIntegratedGrammage(t, t.getLength());
GrammageType const dX = p.getNode()->getModelProperties().getIntegratedGrammage(t);
CORSIKA_LOG_TRACE("EnergyLoss pid={}, z={}, dX={} g/cm2", p.getPID(),
p.getChargeNumber(), dX / 1_g * square(1_cm));
HEPEnergyType dE = getTotalEnergyLoss(p, dX);
......
corsika/detail/modules/proposal/ContinuousProcess.inl
+ 1
− 2
View file @ b7deecfa
......@@ -97,8 +97,7 @@ namespace corsika::proposal {
if (vT.getLength() == 0_m) return ProcessReturn::Ok;
// calculate passed grammage
auto dX =
vP.getNode()->getModelProperties().getIntegratedGrammage(vT, vT.getLength());
auto dX = vP.getNode()->getModelProperties().getIntegratedGrammage(vT);
// get or build corresponding track integral calculator and solve the
// integral
......
corsika/media/BaseExponential.hpp
+ 13
− 8
View file @ b7deecfa
......@@ -23,9 +23,20 @@ namespace corsika {
*/
template <typename TDerived>
class BaseExponential {
public:
BaseExponential(Point const& point, LengthType const referenceHeight,
MassDensityType rho0, LengthType lambda);
Point const& getAnchorPoint() const { return point_; }
MassDensityType getRho0() const { return rho0_; }
InverseLengthType getInvLambda() const { return invLambda_; }
protected:
auto const& getImplementation() const;
MassDensityType getMassDensity(LengthType const height) const;
// clang-format off
/**
* For a (normalized) axis \f$ \vec{a} \f$, the grammage along a non-orthogonal line with (normalized)
......@@ -40,7 +51,7 @@ namespace corsika {
* \f]
*/
// clang-format on
GrammageType getIntegratedGrammage(BaseTrajectory const& line, LengthType vL,
GrammageType getIntegratedGrammage(BaseTrajectory const& line,
DirectionVector const& axis) const;
// clang-format off
......@@ -64,18 +75,12 @@ namespace corsika {
LengthType getArclengthFromGrammage(BaseTrajectory const& line, GrammageType grammage,
DirectionVector const& axis) const;
public:
BaseExponential(Point const& point, MassDensityType rho0, LengthType lambda);
Point const& getAnchorPoint() const { return point_; }
MassDensityType getRho0() const { return rho0_; }
InverseLengthType getInvLambda() const { return invLambda_; }
private:
MassDensityType const rho0_;
LengthType const lambda_;
InverseLengthType const invLambda_;
Point const point_;
LengthType const referenceHeight_;
}; // class BaseExponential
......
corsika/media/FlatExponential.hpp
+ 1
− 2
View file @ b7deecfa
......@@ -41,8 +41,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
GrammageType getIntegratedGrammage(BaseTrajectory const& line,
LengthType to) const override;
GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& line,
GrammageType grammage) const override;
......
corsika/media/HomogeneousMedium.hpp
+ 1
− 2
View file @ b7deecfa
......@@ -30,8 +30,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
GrammageType getIntegratedGrammage(BaseTrajectory const&,
LengthType to) const override;
GrammageType getIntegratedGrammage(BaseTrajectory const&) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const&,
GrammageType grammage) const override;
......
corsika/media/IMediumModel.hpp
+ 19
− 5
View file @ b7deecfa
......@@ -21,11 +21,25 @@ namespace corsika {
virtual MassDensityType getMassDensity(Point const&) const = 0;
// todo: think about the mixin inheritance of the trajectory vs the BaseTrajectory
// approach; for now, only lines are supported
virtual GrammageType getIntegratedGrammage(BaseTrajectory const&,
LengthType) const = 0;
/**
* Integrate the matter density along trajectory.
*
* @return GrammageType as integrated matter density along the BaseTrajectory
*
* @todo think about the mixin inheritance of the trajectory vs the BaseTrajectory
* approach; for now, only lines are supported (?).
*/
virtual GrammageType getIntegratedGrammage(BaseTrajectory const&) const = 0;
/**
* Calculates the length along the trajectory.
*
* The length along the trajectory is determined at which the integrated matter
* density is reached. If the specified matter density cannot be reached (is too
* large) the result becomes meaningless and could be "infinity" (discuss this).
*
* @return LengthType the length corresponding to grammage.
*/
virtual LengthType getArclengthFromGrammage(BaseTrajectory const&,
GrammageType) const = 0;
......
corsika/media/InhomogeneousMedium.hpp
+ 1
− 2
View file @ b7deecfa
......@@ -32,8 +32,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
GrammageType getIntegratedGrammage(BaseTrajectory const& line,
LengthType to) const override;
GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& pLine,
GrammageType grammage) const override;
......
corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+ 4
− 0
View file @ b7deecfa
......@@ -78,6 +78,10 @@ namespace corsika {
void addExponentialLayer(GrammageType b, LengthType c, LengthType upperBoundary);
void addLinearLayer(LengthType c, LengthType upperBoundary);
void addTabularLayer(GrammageType b,
std::function<MassDensityType(LengthType)> const& funcRho,
LengthType upperBoundary);
int getSize() const { return layers_.size(); }
void assemble(Environment<TMediumInterface>& env);
......
corsika/media/LinearApproximationIntegrator.hpp
+ 1
− 1
View file @ b7deecfa
......@@ -20,7 +20,7 @@ namespace corsika {
auto const& getImplementation() const;
public:
auto getIntegrateGrammage(BaseTrajectory const& line, LengthType length) const;
auto getIntegrateGrammage(BaseTrajectory const& line) const;
auto getArclengthFromGrammage(BaseTrajectory const& line,
GrammageType grammage) const;
......
corsika/media/SlidingPlanarExponential.hpp
+ 1
− 2
View file @ b7deecfa
......@@ -46,8 +46,7 @@ namespace corsika {
NuclearComposition const& getNuclearComposition() const override;
GrammageType getIntegratedGrammage(BaseTrajectory const& line,
LengthType l) const override;
GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override;
LengthType getArclengthFromGrammage(BaseTrajectory const& line,
GrammageType grammage) const override;
......
tests/media/testEnvironment.cpp
+ 40
− 29
View file @ b7deecfa
......@@ -97,33 +97,35 @@ TEST_CASE("FlatExponential") {
auto const rho0 = 1_g / static_pow<3>(1_cm);
FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
protonComposition);
auto const tEnd = 5_s;
SpeedType const speed = 20_m / second;
LengthType const length = 2_m;
TimeType const tEnd = length / speed;
SECTION("horizontal") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {20_cm / second, 0_m / second, 0_m / second}));
gCS, {speed, 0_m / second, 0_m / second}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
CHECK((medium.getIntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
CHECK((medium.getIntegratedGrammage(trajectory) / (rho0 * length)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, rho0 * length) / length) ==
Approx(1));
}
SECTION("vertical") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, 5_m / second}));
gCS, {0_m / second, 0_m / second, speed}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
LengthType const length = 2 * lambda;
GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1));
CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
SECTION("escape grammage") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 0_m / second, -5_m / second}));
gCS, {SpeedType::zero(), SpeedType::zero(), -speed}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
GrammageType const escapeGrammage = rho0 * lambda;
......@@ -134,15 +136,15 @@ TEST_CASE("FlatExponential") {
}
SECTION("inclined") {
Line const line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {0_m / second, 5_m / second, 5_m / second}));
Line const line(gOrigin,
Vector<SpeedType::dimension_type>(
gCS, {0_m / second, speed / sqrt(2.), speed / sqrt(2.)}));
setup::Trajectory const trajectory =
setup::testing::make_track<setup::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;
CHECK((medium.getIntegratedGrammage(trajectory, length) / exact) == Approx(1));
CHECK((medium.getIntegratedGrammage(trajectory) / exact) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
}
}
......@@ -179,16 +181,16 @@ TEST_CASE("SlidingPlanarExponential") {
CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude() ==
flat.getMassDensity({gCS, {0_mm, 0_m, 3_m}}).magnitude());
CHECK(medium.getIntegratedGrammage(trajectory, 2_m).magnitude() ==
flat.getIntegratedGrammage(trajectory, 2_m).magnitude());
CHECK(medium.getIntegratedGrammage(trajectory).magnitude() ==
flat.getIntegratedGrammage(trajectory).magnitude());
CHECK(medium.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude() ==
flat.getArclengthFromGrammage(trajectory, rho0 * 5_m).magnitude());
}
}
auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
MassDensityType constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
struct Exponential {
struct ExponentialTest {
auto operator()(Point const& p) const {
return exp(p.getCoordinates()[0] / 1_m) * rho0;
}
......@@ -213,41 +215,50 @@ TEST_CASE("InhomogeneousMedium") {
Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
SpeedType const speed = 20_m / second;
Line line(gOrigin, Vector<SpeedType::dimension_type>(
gCS, {20_m / second, 0_m / second, 0_m / second}));
gCS, {speed, SpeedType::zero(), SpeedType::zero()}));
auto const tEnd = 5_s;
// the tested LinearApproximationIntegrator really does a single step only. It is very
// poor for exponentials with a bit larger step-width.
TimeType const tEnd = 0.001_s;
setup::Trajectory const trajectory =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
Exponential const e;
DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
ExponentialTest const expTest;
DensityFunction<ExponentialTest, LinearApproximationIntegrator> const rho(expTest);
SECTION("DensityFunction") {
CHECK(e.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
CHECK(expTest.getDerivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
Approx(1));
CHECK(rho.evaluateAt(gOrigin) == e(gOrigin));
CHECK(rho.evaluateAt(gOrigin) == expTest(gOrigin));
}
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;
LengthType constexpr length = tEnd * speed;
NuclearComposition const composition{{Code::Proton}, {1.f}};
InhomogeneousMedium<IMediumModel, decltype(rho)> const inhMedium(composition, rho);
CORSIKA_LOG_INFO("test={} l={} {} {}", rho.getIntegrateGrammage(trajectory), length,
exactGrammage(length), 1_m * rho0 * (exp(length / 1_m) - 1));
SECTION("Integration") {
CHECK(rho.getIntegrateGrammage(trajectory, l) / exactGrammage(l) ==
CORSIKA_LOG_INFO("test={} {} {}", rho.getIntegrateGrammage(trajectory),
exactGrammage(length),
rho.getIntegrateGrammage(trajectory) / exactGrammage(length));
CHECK(rho.getIntegrateGrammage(trajectory) / exactGrammage(length) ==
Approx(1).epsilon(1e-2));
CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(l)) /
exactLength(exactGrammage(l)) ==
CHECK(rho.getArclengthFromGrammage(trajectory, exactGrammage(length)) /
exactLength(exactGrammage(length)) ==
Approx(1).epsilon(1e-2));
CHECK(rho.getMaximumLength(trajectory, 1e-2) >
l); // todo: write reasonable test when implementation is working
length); // todo: write reasonable test when implementation is working
CHECK(rho.getIntegrateGrammage(trajectory, l) ==
inhMedium.getIntegratedGrammage(trajectory, l));
CHECK(rho.getIntegrateGrammage(trajectory) ==
inhMedium.getIntegratedGrammage(trajectory));
CHECK(rho.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)) ==
inhMedium.getArclengthFromGrammage(trajectory, 20_g / (1_cm * 1_cm)));
CHECK(inhMedium.getNuclearComposition() == composition);
......
tests/media/testMagneticField.cpp
+ 0
− 5
View file @ b7deecfa
......@@ -25,7 +25,6 @@ using namespace corsika;
TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
logging::set_level(logging::level::info);
corsika_logger->set_pattern("[%n:%^%-8l%$] custom pattern: %v");
CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
......@@ -85,8 +84,4 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
// and the associated trajectory
setup::Trajectory const track =
setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.getIntegratedGrammage(track, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(track, density * 5_m) / 5_m) == Approx(1));
}
tests/media/testMedium.cpp
+ 6
− 3
View file @ b7deecfa
......@@ -92,18 +92,21 @@ TEST_CASE("MediumPropertyModel w/ Homogeneous") {
CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
medium.getNuclearComposition();
SpeedType const speed = 1_m / second;
// create a line of length 1 m
Line const line(gOrigin,
VelocityVector(gCS, {1_m / second, 0_m / second, 0_m / second}));
Line const line(gOrigin, VelocityVector(gCS, {speed, 0_m / second, 0_m / second}));
// the end time of our line
auto const tEnd = 1_s;
LengthType const length = tEnd * speed;
// and the associated trajectory
setup::Trajectory const trajectory =
corsika::setup::testing::make_track<setup::Trajectory>(line, tEnd);
// and check the integrated grammage
CHECK((medium.getIntegratedGrammage(trajectory, 3_m) / (density * 3_m)) == Approx(1));
CHECK((medium.getIntegratedGrammage(trajectory) / (density * length)) == Approx(1));
CHECK((medium.getArclengthFromGrammage(trajectory, density * 5_m) / 5_m) == Approx(1));
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum