diff --git a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
index 791d57972d8496fb130bd87f1673a22686ff8557..17371ff3b775986099e6436894c1213a945e4733 100644
--- a/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
+++ b/corsika/detail/media/LayeredSphericalAtmosphereBuilder.inl
@@ -13,7 +13,7 @@
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
-//#include <corsika/media/SlidingPlanarTabular.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
namespace corsika {
@@ -42,7 +42,7 @@ namespace corsika {
TModelArgs...>::addExponentialLayer(GrammageType b, LengthType c,
LengthType upperBoundary) {
- auto const radius = earthRadius_ + upperBoundary;
+ auto const radius = planetRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
@@ -56,7 +56,7 @@ namespace corsika {
auto lastBound = [&](auto... argPack) {
return std::make_shared<
TMediumModelExtra<SlidingPlanarExponential<TMediumInterface>>>(
- argPack..., center_, rho0, -c, *composition_, earthRadius_);
+ argPack..., center_, rho0, -c, *composition_, planetRadius_);
};
// now unpack the additional arguments
@@ -64,7 +64,7 @@ namespace corsika {
node->setModelProperties(std::move(model));
} else {
node->template setModelProperties<SlidingPlanarExponential<TMediumInterface>>(
- center_, rho0, -c, *composition_, earthRadius_);
+ center_, rho0, -c, *composition_, planetRadius_);
}
layers_.push(std::move(node));
@@ -75,14 +75,14 @@ namespace corsika {
inline void LayeredSphericalAtmosphereBuilder<
TMediumInterface, TMediumModelExtra,
TModelArgs...>::addLinearLayer(LengthType c, LengthType upperBoundary) {
- auto const radius = earthRadius_ + upperBoundary;
+ auto const radius = planetRadius_ + upperBoundary;
checkRadius(radius);
previousRadius_ = radius;
auto node = std::make_unique<VolumeTreeNode<TMediumInterface>>(
std::make_unique<Sphere>(center_, radius));
- units::si::GrammageType constexpr b = 1 * 1_g / (1_cm * 1_cm);
+ units::si::GrammageType constexpr b = 1_g / (1_cm * 1_cm);
auto const rho0 = b / c;
if constexpr (detail::has_extra_models<TMediumModelExtra>::value) {
@@ -94,7 +94,6 @@ namespace corsika {
// now unpack the additional arguments
auto model = std::apply(lastBound, additionalModelArgs_);
-
node->setModelProperties(std::move(model));
} else {
node->template setModelProperties<HomogeneousMedium<TMediumInterface>>(
@@ -108,33 +107,33 @@ namespace corsika {
typename... TModelArgs>
inline void
LayeredSphericalAtmosphereBuilder<TMediumInterface, TMediumModelExtra, TModelArgs...>::
- addTabularLayer(GrammageType b,
- std::function<MassDensityType(LengthType)> const& funcRho,
- LengthType upperBoundary) {
+ addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+ unsigned int const nBins, LengthType const deltaHeight,
+ LengthType const upperBoundary) {
- auto const radius = earthRadius_ + upperBoundary;
+ auto const radius = planetRadius_ + 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_);
- }
- */
+
+ 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_, funcRho, nBins, deltaHeight, *composition_,
+ planetRadius_);
+ };
+
+ // now unpack the additional arguments
+ auto model = std::apply(lastBound, additionalModelArgs_);
+ node->setModelProperties(std::move(model));
+ } else {
+ node->template setModelProperties<SlidingPlanarTabular<TMediumInterface>>(
+ center_, funcRho, nBins, deltaHeight, *composition_, planetRadius_);
+ }
layers_.push(std::move(node));
}
@@ -167,10 +166,10 @@ namespace corsika {
template <typename TMediumInterface, template <typename> typename MExtraEnvirnoment>
struct make_layered_spherical_atmosphere_builder {
template <typename... TArgs>
- static auto create(Point const& center, LengthType earthRadius, TArgs... args) {
+ static auto create(Point const& center, LengthType planetRadius, TArgs... args) {
return LayeredSphericalAtmosphereBuilder<TMediumInterface, MExtraEnvirnoment,
TArgs...>{std::forward<TArgs>(args)...,
- center, earthRadius};
+ center, planetRadius};
}
};
diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl
index 80fc9bf6b77624e771b4df0637903660e9a707d8..9b921872811e157c9d0d232286b3714444983727 100644
--- a/corsika/detail/media/SlidingPlanarExponential.inl
+++ b/corsika/detail/media/SlidingPlanarExponential.inl
@@ -8,50 +8,51 @@
#pragma once
-#include <corsika/media/SlidingPlanarExponential.hpp>
-
namespace corsika {
- template <typename T>
- inline SlidingPlanarExponential<T>::SlidingPlanarExponential(
+ template <typename TDerived>
+ inline SlidingPlanarExponential<TDerived>::SlidingPlanarExponential(
Point const& p0, MassDensityType rho0, LengthType lambda,
NuclearComposition const& nuclComp, LengthType referenceHeight)
- : BaseExponential<SlidingPlanarExponential<T>>(p0, referenceHeight, rho0, lambda)
- , nuclComp_(nuclComp)
- , referenceHeight_(referenceHeight) {}
+ : BaseExponential<SlidingPlanarExponential<TDerived>>(p0, referenceHeight, rho0,
+ lambda)
+ , nuclComp_(nuclComp) {}
- template <typename T>
- inline MassDensityType SlidingPlanarExponential<T>::getMassDensity(
+ template <typename TDerived>
+ inline MassDensityType SlidingPlanarExponential<TDerived>::getMassDensity(
Point const& point) const {
auto const heightFull =
- (point - BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
+ (point - BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
.getNorm();
- return BaseExponential<SlidingPlanarExponential<T>>::getMassDensity(heightFull);
+ return BaseExponential<SlidingPlanarExponential<TDerived>>::getMassDensity(
+ heightFull);
}
- template <typename T>
- inline NuclearComposition const& SlidingPlanarExponential<T>::getNuclearComposition()
- const {
+ template <typename TDerived>
+ inline NuclearComposition const&
+ SlidingPlanarExponential<TDerived>::getNuclearComposition() const {
return nuclComp_;
}
- template <typename T>
- inline GrammageType SlidingPlanarExponential<T>::getIntegratedGrammage(
+ template <typename TDerived>
+ inline GrammageType SlidingPlanarExponential<TDerived>::getIntegratedGrammage(
BaseTrajectory const& traj) const {
- auto const axis = (traj.getPosition(0) -
- BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
- .normalized();
- return BaseExponential<SlidingPlanarExponential<T>>::getIntegratedGrammage(traj,
- axis);
+ auto const axis =
+ (traj.getPosition(0) -
+ BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+ .normalized();
+ return BaseExponential<SlidingPlanarExponential<TDerived>>::getIntegratedGrammage(
+ traj, axis);
}
- template <typename T>
- inline LengthType SlidingPlanarExponential<T>::getArclengthFromGrammage(
+ template <typename TDerived>
+ inline LengthType SlidingPlanarExponential<TDerived>::getArclengthFromGrammage(
BaseTrajectory const& traj, GrammageType const grammage) const {
- auto const axis = (traj.getPosition(0) -
- BaseExponential<SlidingPlanarExponential<T>>::getAnchorPoint())
- .normalized();
- return BaseExponential<SlidingPlanarExponential<T>>::getArclengthFromGrammage(
+ auto const axis =
+ (traj.getPosition(0) -
+ BaseExponential<SlidingPlanarExponential<TDerived>>::getAnchorPoint())
+ .normalized();
+ return BaseExponential<SlidingPlanarExponential<TDerived>>::getArclengthFromGrammage(
traj, grammage, axis);
}
diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp
index 27a267df6061947843d55bac73362adea9ca0f6f..507ed640fa130eda007a5e71a3dec3437a530d9b 100644
--- a/corsika/media/BaseExponential.hpp
+++ b/corsika/media/BaseExponential.hpp
@@ -29,8 +29,6 @@ namespace corsika {
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;
diff --git a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
index 41ad158011dfb7215abca026a84b082736ccc0e0..6c178f4589537fff2a79c7e166662791c30a1306 100644
--- a/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
+++ b/corsika/media/LayeredSphericalAtmosphereBuilder.hpp
@@ -22,6 +22,7 @@
#include <stack>
#include <tuple>
#include <type_traits>
+#include <functional>
namespace corsika {
@@ -68,9 +69,9 @@ namespace corsika {
protected:
LayeredSphericalAtmosphereBuilder(TModelArgs... args, Point const& center,
- LengthType earthRadius)
+ LengthType planetRadius)
: center_(center)
- , earthRadius_(earthRadius)
+ , planetRadius_(planetRadius)
, additionalModelArgs_{args...} {}
public:
@@ -78,9 +79,9 @@ 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);
+ void addTabularLayer(std::function<MassDensityType(LengthType)> const& funcRho,
+ unsigned int const nBins, LengthType const deltaHeight,
+ LengthType const upperBoundary);
int getSize() const { return layers_.size(); }
@@ -88,9 +89,9 @@ namespace corsika {
Environment<TMediumInterface> assemble();
/**
- * Get the current Earth radius.
+ * Get the current planet radius.
*/
- LengthType getEarthRadius() const { return earthRadius_; }
+ LengthType getPlanetRadius() const { return planetRadius_; }
private:
void checkRadius(LengthType r) const;
@@ -98,7 +99,7 @@ namespace corsika {
std::unique_ptr<NuclearComposition> composition_;
Point center_;
LengthType previousRadius_{LengthType::zero()};
- LengthType earthRadius_;
+ LengthType planetRadius_;
std::tuple<TModelArgs...> const additionalModelArgs_;
std::stack<typename VolumeTreeNode<TMediumInterface>::VTNUPtr>
diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp
index 63f8238ba2139899957fb4ff12d6e1e1b6b0f7e0..4ef91bb4e1ca338241c61d805f1fa3dd68c47ba4 100644
--- a/corsika/media/SlidingPlanarExponential.hpp
+++ b/corsika/media/SlidingPlanarExponential.hpp
@@ -13,7 +13,6 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/random/RNGManager.hpp>
-#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/framework/geometry/BaseTrajectory.hpp>
@@ -53,7 +52,6 @@ namespace corsika {
private:
NuclearComposition const nuclComp_;
- LengthType const referenceHeight_;
};
} // namespace corsika
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index e3d056e22202e7af786c3f1eec00fcc36e077a61..70816c0d30968b67caa65f98b204c4ab0b09c338 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -60,4 +60,6 @@ CORSIKA_REGISTER_EXAMPLE (hybrid_MC RUN_OPTIONS 4 2 10000.)
add_executable (corsika corsika.cpp)
target_link_libraries (corsika CORSIKA8::CORSIKA8)
+add_executable (mars mars.cpp)
+target_link_libraries (mars CORSIKA8::CORSIKA8)
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index fd048470df3c3d4c598ab35c0be67b3db0ebf0e1..ee424c1de386a288b934c1f6990e1cb44ff1e214 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -103,8 +103,6 @@ using MyExtraEnv = MediumPropertyModel<UniformMagneticField<T>>;
int main(int argc, char** argv) {
- logging::set_level(logging::level::info);
-
// the main command line description
CLI::App app{"Simulate standard (downgoing) showers with CORSIKA 8."};
@@ -154,10 +152,31 @@ int main(int argc, char** argv) {
->group("Misc.");
app.add_flag("--force-interaction", "Force the location of the first interaction.")
->group("Misc.");
+ app.add_option("-v,--verbosity", "Verbosity level: warn, info, debug, trace.")
+ ->default_val("info")
+ ->check(CLI::IsMember({"warn", "info", "debug", "trace"}))
+ ->group("Misc.");
// parse the command line options into the variables
CLI11_PARSE(app, argc, argv);
+ if (app.count("--verbosity")) {
+ string const loglevel = app["verbosity"]->as<string>();
+ if (loglevel == "warn") {
+ logging::set_level(logging::level::warn);
+ } else if (loglevel == "info") {
+ logging::set_level(logging::level::info);
+ } else if (loglevel == "debug") {
+ logging::set_level(logging::level::debug);
+ } else if (loglevel == "trace") {
+#ifndef DEBUG
+ CORSIKA_LOG_ERROR("trace log level requires a Debug build.");
+ return 1;
+#endif
+ logging::set_level(logging::level::trace);
+ }
+ }
+
// check that we got either PDG or A/Z
// this can be done with option_groups but the ordering
// gets all messed up
@@ -185,17 +204,29 @@ int main(int argc, char** argv) {
50_uT, 0_T});
builder.setNuclearComposition(
{{Code::Nitrogen, Code::Oxygen},
- {0.7847f, 1.f - 0.7847f}}); // values taken from AIRES manual, Ar removed for now
+ {0.7847, 1. - 0.7847}}); // values taken from AIRES manual, Ar removed for now
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 2_km);
builder.addExponentialLayer(1222.6562_g / (1_cm * 1_cm), 994186.38_cm, 4_km);
builder.addExponentialLayer(1144.9069_g / (1_cm * 1_cm), 878153.55_cm, 10_km);
builder.addExponentialLayer(1305.5948_g / (1_cm * 1_cm), 636143.04_cm, 40_km);
builder.addExponentialLayer(540.1778_g / (1_cm * 1_cm), 772170.16_cm, 100_km);
- builder.addLinearLayer(1e9_cm, 112.8_km + constants::EarthRadius::Mean);
+ builder.addLinearLayer(1e9_cm, 112.8_km);
builder.assemble(env);
/* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
+ ofstream atmout("earth.dat");
+ for (LengthType h = 0_m; h < 110_km; h += 100_m) {
+ Point const ptest{rootCS, 0_m, 0_m, builder.getPlanetRadius() + h};
+ auto rho =
+ env.getUniverse()->getContainingNode(ptest)->getModelProperties().getMassDensity(
+ ptest);
+ atmout << h / 1_m << " " << rho / 1_kg * cube(1_m) << "\n";
+ }
+ atmout.close();
+
+ /* === START: CONSTRUCT PRIMARY PARTICLE === */
+
/* === START: CONSTRUCT PRIMARY PARTICLE === */
// parse the primary ID as a PDG or A/Z code
@@ -231,8 +262,8 @@ int main(int argc, char** argv) {
/* === END: CONSTRUCT PRIMARY PARTICLE === */
/* === START: CONSTRUCT GEOMETRY === */
- auto const observationHeight = 0_km + builder.getEarthRadius();
- auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+ auto const observationHeight = 0_km + builder.getPlanetRadius();
+ auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
@@ -281,9 +312,9 @@ int main(int argc, char** argv) {
Code::KStar0_1430_MinusBar,
}};
- decaySibyll.printDecayConfig();
+ // decaySibyll.printDecayConfig();
- ParticleCut cut{50_GeV, 50_GeV, 50_GeV, 50_GeV, false};
+ ParticleCut cut{1_GeV, 1_GeV, 1_GeV, 1_GeV, false};
corsika::proposal::Interaction emCascade(env);
corsika::proposal::ContinuousProcess emContinuous(env);
InteractionCounter emCascadeCounted(emCascade);
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index 0ff4d28b7425d1fd988af0718ccecbe84b4b8ba1..d8b8385f3696ff1fee4f48518253bd24e97fbdb0 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -128,8 +128,8 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
- auto const observationHeight = 1.4_km + builder.getEarthRadius();
- auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+ auto const observationHeight = 1.4_km + builder.getPlanetRadius();
+ auto const injectionHeight = 112.75_km + builder.getPlanetRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/examples/hybrid_MC.cpp b/examples/hybrid_MC.cpp
index e065c4ee8ddc787e2bd56d2f8502f37cfb7d0bd1..4674512cc4d0f204538999c6bf6a3eb384704ea2 100644
--- a/examples/hybrid_MC.cpp
+++ b/examples/hybrid_MC.cpp
@@ -178,8 +178,8 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
- auto const observationHeight = 0_km + builder.getEarthRadius();
- auto const injectionHeight = 112.75_km + builder.getEarthRadius();
+ auto const observationHeight = 0_km + builder.getPlanetRadius();
+ auto const injectionHeight = 112.75_km + builder.getPlanetRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 6ae766bb197e2e6610aa503442cea05704652b3b..6800a84be8858d31ed01d36d213ef56754a265b1 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -177,8 +177,8 @@ int main(int argc, char** argv) {
cout << "input momentum: " << plab.getComponents() / 1_GeV
<< ", norm = " << plab.getNorm() << endl;
- auto const observationHeight = 0_km + builder.getEarthRadius();
- auto const injectionHeight = 111.75_km + builder.getEarthRadius();
+ auto const observationHeight = 0_km + builder.getPlanetRadius();
+ auto const injectionHeight = 111.75_km + builder.getPlanetRadius();
auto const t = -observationHeight * cos(thetaRad) +
sqrt(-static_pow<2>(sin(thetaRad) * observationHeight) +
static_pow<2>(injectionHeight));
diff --git a/tests/media/testEnvironment.cpp b/tests/media/testEnvironment.cpp
index 88e143b7be1d66e0e349284882cedb6002636fff..fffdff7fe04389a0feb7d9471a83ed3b7740d145 100644
--- a/tests/media/testEnvironment.cpp
+++ b/tests/media/testEnvironment.cpp
@@ -11,6 +11,7 @@
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+
#include <corsika/media/DensityFunction.hpp>
#include <corsika/media/FlatExponential.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
@@ -26,6 +27,7 @@
#include <corsika/media/LinearApproximationIntegrator.hpp>
#include <corsika/media/NuclearComposition.hpp>
#include <corsika/media/SlidingPlanarExponential.hpp>
+#include <corsika/media/SlidingPlanarTabular.hpp>
#include <corsika/media/VolumeTreeNode.hpp>
#include <SetupTestTrajectory.hpp>
@@ -94,7 +96,7 @@ TEST_CASE("FlatExponential") {
Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
LengthType const lambda = 3_m;
- auto const rho0 = 1_g / static_pow<3>(1_cm);
+ auto const rho0 = 1_g / cube(1_cm);
FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
protonComposition);
SpeedType const speed = 20_m / second;
@@ -188,6 +190,256 @@ TEST_CASE("SlidingPlanarExponential") {
}
}
+struct RhoFuncConst {
+ MassDensityType operator()(LengthType) const { return 1_g / cube(1_cm); }
+ static GrammageType integrate(LengthType dL) { return dL * 1_g / cube(1_cm); }
+};
+
+struct RhoFuncExp {
+ MassDensityType operator()(LengthType height) const {
+ return 1_g / cube(1_cm) * exp(-height / 1000_m);
+ }
+ static GrammageType integrate(BaseTrajectory const& traj, Point const& origin,
+ LengthType const& refH) {
+ LengthType height1 = (traj.getPosition(0) - origin).getNorm() - refH;
+ LengthType height2 = (traj.getPosition(1) - origin).getNorm() - refH;
+ if (height1 > height2) { std::swap(height1, height2); }
+
+ DirectionVector const axis(
+ (traj.getPosition(0) - origin).normalized()); // to gravity center
+ double const cosTheta = axis.dot(traj.getDirection(0));
+
+ CORSIKA_LOG_INFO("h1={} h2={} cT={} rho1={}, rho2={}", height1, height2, cosTheta,
+ 1_g / cube(1_cm) * exp(-height1 / 1000_m),
+ 1_g / cube(1_cm) * exp(-height2 / 1000_m));
+ return (1_km * 1_g / cube(1_cm) * exp(-height1 / 1000_m) -
+ 1_km * 1_g / cube(1_cm) * exp(-height2 / 1000_m)) /
+ cosTheta;
+ }
+ static GrammageType integrate(BaseTrajectory const& traj, LengthType const& length,
+ Point const& origin, LengthType const& refH) {
+ LengthType height1 = (traj.getPosition(0) - origin).getNorm() - refH;
+ LengthType height2 =
+ (traj.getPosition(0) + traj.getDirection(0) * length - origin).getNorm() - refH;
+ if (height1 > height2) { std::swap(height1, height2); }
+
+ DirectionVector const axis(
+ (traj.getPosition(0) - origin).normalized()); // to gravity center
+ double const cosTheta = axis.dot(traj.getDirection(0));
+
+ CORSIKA_LOG_INFO("h1={} h2={} cT={}", height1, height2, cosTheta);
+ return (1_km * 1_g / cube(1_cm) * exp(-height1 / 1000_m) -
+ 1_km * 1_g / cube(1_cm) * exp(-height2 / 1000_m)) /
+ cosTheta;
+ }
+};
+
+TEST_CASE("SlidingPlanarTabular") {
+
+ logging::set_level(logging::level::info);
+
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ RhoFuncConst rhoFunc;
+ SlidingPlanarTabular<IMediumModel> const medium(gOrigin, rhoFunc, 1000, 10_m,
+ protonComposition);
+
+ SECTION("density") {
+ CHECK(medium.getMassDensity({gCS, {0_m, 0_m, 3_m}}) /
+ medium.getMassDensity({gCS, {0_m, 3_m, 0_m}}) ==
+ Approx(1));
+ CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 3_m}}) == 1_g / cube(1_cm));
+ CHECK(medium.getMassDensity({gCS, {0_mm, 0_m, 300_m}}) == 1_g / cube(1_cm));
+ }
+
+ SECTION("vertical") {
+ SpeedType const speed = 5_m / second;
+ TimeType const tEnd1 = 1_s;
+ LengthType const length1 = speed * tEnd1;
+ TimeType const tEnd2 = 300_s;
+ LengthType const length2 = speed * tEnd2;
+ Line const line(
+ {gCS, {0_m, 0_m, 1_m}},
+ Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, speed}));
+ setup::Trajectory const trajectory1 =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd1);
+ Line const line1Reverse(
+ trajectory1.getPosition(1),
+ Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, -speed}));
+ setup::Trajectory const trajectory1Reverse =
+ setup::testing::make_track<setup::Trajectory>(line1Reverse, tEnd1);
+
+ setup::Trajectory const trajectory2 =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd2);
+ Line const line2Reverse(
+ trajectory2.getPosition(0),
+ Vector<SpeedType::dimension_type>(gCS, {0_m / second, 0_m / second, -speed}));
+ setup::Trajectory const trajectory2Reverse =
+ setup::testing::make_track<setup::Trajectory>(line2Reverse, tEnd2);
+
+ // failures
+ CHECK_THROWS(medium.getArclengthFromGrammage(trajectory1, -1_kg / square(1_cm)));
+
+ MassDensityType const rho0 = 1_g / cube(1_cm);
+
+ // short track
+ CHECK(medium.getIntegratedGrammage(trajectory1) == length1 * rho0);
+ LengthType const testD1 = length1 / 200; // within bin
+ CHECK(medium.getArclengthFromGrammage(trajectory1, rho0 * testD1) / testD1 ==
+ Approx(1));
+ // short track, reverse
+ CHECK(medium.getIntegratedGrammage(trajectory1Reverse) == length1 * rho0);
+ CHECK(medium.getArclengthFromGrammage(trajectory1Reverse, rho0 * testD1) / testD1 ==
+ Approx(1));
+
+ // long track
+ CHECK(medium.getIntegratedGrammage(trajectory2) == length2 * 1_g / cube(1_cm));
+ LengthType const testD2 = length2 / 25; // multi bin
+ CHECK(medium.getArclengthFromGrammage(trajectory2, rho0 * testD2) == testD2);
+ }
+
+ SECTION("inclined") {
+ SpeedType const speed = 5_m / second;
+ TimeType const tEnd1 = 1_s;
+ LengthType const length1 = speed * tEnd1;
+ TimeType const tEnd2 = 300_s;
+ LengthType const length2 = speed * tEnd2;
+ Line const line({gCS, {0_m, 0_m, 1_m}},
+ Vector<SpeedType::dimension_type>(
+ gCS, {speed / sqrt(2.), 0_m / second, speed / sqrt(2.)}));
+ setup::Trajectory const trajectory1 =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd1);
+ Line const line1Reverse(
+ trajectory1.getPosition(1),
+ Vector<SpeedType::dimension_type>(
+ gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+ setup::Trajectory const trajectory1Reverse =
+ setup::testing::make_track<setup::Trajectory>(line1Reverse, tEnd1);
+
+ setup::Trajectory const trajectory2 =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd2);
+ Line const line2Reverse(
+ trajectory2.getPosition(1),
+ Vector<SpeedType::dimension_type>(
+ gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+ setup::Trajectory const trajectory2Reverse =
+ setup::testing::make_track<setup::Trajectory>(line2Reverse, tEnd2);
+
+ MassDensityType const rho0 = 1_g / cube(1_cm);
+
+ // short track
+ CHECK(medium.getIntegratedGrammage(trajectory1) / (length1 * rho0) == Approx(1));
+ LengthType const testD1 = length1 / 200; // within bin
+ CHECK(medium.getArclengthFromGrammage(trajectory1, RhoFuncConst::integrate(testD1)) /
+ testD1 ==
+ Approx(1));
+ // short track, reverse
+ CHECK(medium.getIntegratedGrammage(trajectory1Reverse) / (length1 * rho0) ==
+ Approx(1));
+ CHECK(medium.getArclengthFromGrammage(trajectory1Reverse, rho0 * testD1) / testD1 ==
+ Approx(1));
+
+ // long track
+ CHECK(medium.getIntegratedGrammage(trajectory2) / (length2 * rho0) ==
+ Approx(1).epsilon(0.01));
+ LengthType const testD2 = length2 / 25; // multi bin
+ CHECK(medium.getArclengthFromGrammage(trajectory2, rho0 * testD2) / testD2 ==
+ Approx(1).epsilon(0.01));
+ // long track reverse
+ CORSIKA_LOG_INFO("length2={}", length2);
+ CHECK(medium.getIntegratedGrammage(trajectory2Reverse) / (length2 * rho0) ==
+ Approx(1).epsilon(0.01));
+ CHECK(medium.getArclengthFromGrammage(trajectory2Reverse, rho0 * testD2) / testD2 ==
+ Approx(1).epsilon(0.01));
+ }
+
+ /*The exponential test is taken over phase-space where the exponential is not so steep
+ * and is samples in sufficient substeps. An reference-height offset of 1000_km is used.
+ * Thus, density is given from 1000 to 1010 km. And curvature effects are small.
+ */
+
+ RhoFuncExp rhoFuncExp;
+ SlidingPlanarTabular<IMediumModel> const mediumExp(gOrigin, rhoFuncExp, 1000, 10_m,
+ protonComposition, 1000_km);
+
+ SECTION("exponential") {
+
+ SpeedType const speed = 5_m / second;
+ TimeType const tEnd1 = 1_s;
+ LengthType const length1 = speed * tEnd1;
+ TimeType const tEnd2 = 300_s;
+ LengthType const length2 = speed * tEnd2;
+ Line const line({gCS, {0_m, 0_m, 1000.005_km}},
+ Vector<SpeedType::dimension_type>(
+ gCS, {speed / sqrt(2.), 0_m / second, speed / sqrt(2.)}));
+ setup::Trajectory const trajectory1 =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd1);
+ Line const line1Reverse(
+ trajectory1.getPosition(1),
+ Vector<SpeedType::dimension_type>(
+ gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+ setup::Trajectory const trajectory1Reverse =
+ setup::testing::make_track<setup::Trajectory>(line1Reverse, tEnd1);
+
+ setup::Trajectory const trajectory2 =
+ setup::testing::make_track<setup::Trajectory>(line, tEnd2);
+
+ CORSIKA_LOG_INFO("{} {}", RhoFuncExp::integrate(trajectory1, gOrigin, 1000_km),
+ length1);
+
+ // short track
+ GrammageType const testShortX = RhoFuncExp::integrate(trajectory1, gOrigin, 1000_km);
+ CHECK(mediumExp.getIntegratedGrammage(trajectory1) / testShortX ==
+ Approx(1).epsilon(0.01));
+ LengthType const testD1 = length1 / 200; // within bin
+ GrammageType const testD1X =
+ RhoFuncExp::integrate(trajectory1, testD1, gOrigin, 1000_km);
+ CHECK(mediumExp.getArclengthFromGrammage(trajectory1, testD1X) / testD1 ==
+ Approx(1).epsilon(0.01));
+ // short track, reverse
+ CHECK(mediumExp.getIntegratedGrammage(trajectory1Reverse) / testShortX ==
+ Approx(1).epsilon(0.01));
+ CHECK(mediumExp.getArclengthFromGrammage(trajectory1Reverse, testD1X) / testD1 ==
+ Approx(1).epsilon(0.01));
+
+ // long track
+ GrammageType const testLongX = RhoFuncExp::integrate(trajectory2, gOrigin, 1000_km);
+ CORSIKA_LOG_INFO("testLongX={}", testLongX);
+ CHECK(mediumExp.getIntegratedGrammage(trajectory2) / testLongX ==
+ Approx(1).epsilon(0.01));
+ LengthType const testD2 = length2 / 25; // multi bin
+ GrammageType const testD2X =
+ RhoFuncExp::integrate(trajectory2, testD2, gOrigin, 1000_km);
+ CHECK(mediumExp.getArclengthFromGrammage(trajectory2, testD2X) / testD2 ==
+ Approx(1).epsilon(0.01));
+ // long track, reverse
+
+ // first full trajectory2 reverse
+ Line line2Reverse(trajectory2.getPosition(1),
+ Vector<SpeedType::dimension_type>(
+ gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+ setup::Trajectory trajectory2Reverse =
+ setup::testing::make_track<setup::Trajectory>(line2Reverse, tEnd2);
+
+ CHECK(mediumExp.getIntegratedGrammage(trajectory2Reverse) / testLongX ==
+ Approx(1).epsilon(0.01));
+
+ // but now shorter trajectory2 reversed to correspond 100% to testD2
+
+ line2Reverse = Line(trajectory2.getPosition(0) + trajectory2.getDirection(0) * testD2,
+ Vector<SpeedType::dimension_type>(
+ gCS, {-speed / sqrt(2.), 0_m / second, -speed / sqrt(2.)}));
+ auto const trajectory2ReverseShort =
+ setup::testing::make_track<setup::Trajectory>(line2Reverse, testD2 / speed);
+
+ CORSIKA_LOG_INFO("here {} {} {}", trajectory2ReverseShort.getLength(), testD2,
+ testD2X / 1_g * square(1_cm));
+ CHECK(mediumExp.getArclengthFromGrammage(trajectory2ReverseShort, testD2X) / testD2 ==
+ Approx(1).epsilon(0.01));
+ }
+}
+
MassDensityType constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
struct ExponentialTest {
@@ -289,7 +541,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder") {
CHECK(builder.getSize() == 0);
- auto const R = builder.getEarthRadius();
+ auto const R = builder.getPlanetRadius();
CHECK(univ->getChildNodes().size() == 1);
@@ -334,7 +586,7 @@ TEST_CASE("LayeredSphericalAtmosphereBuilder w/ magnetic field") {
CHECK(builder.getSize() == 0);
CHECK(univ->getChildNodes().size() == 1);
- auto const R = builder.getEarthRadius();
+ auto const R = builder.getPlanetRadius();
// check magnetic field at several locations
const Point pTest(gCS, -10_m, 4_m, R + 35_m);