diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt
index ef1d24035859fc6668d902ae89ee7aebc994a95b..47b973cb138cbc7a52ba11e36b8347ec8497e883 100644
--- a/Environment/CMakeLists.txt
+++ b/Environment/CMakeLists.txt
@@ -9,6 +9,7 @@ set (
LinearApproximationIntegrator.h
DensityFunction.h
Environment.h
+ FlatExponential.h
)
set (
diff --git a/Environment/FlatExponential.h b/Environment/FlatExponential.h
new file mode 100644
index 0000000000000000000000000000000000000000..7079a22bfa39f4b061e1bfe111553c34e595d8c7
--- /dev/null
+++ b/Environment/FlatExponential.h
@@ -0,0 +1,85 @@
+
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#ifndef _include_Environment_FlatExponential_h_
+#define _include_Environment_FlatExponential_h_
+
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Line.h>
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Trajectory.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <cassert>
+
+/**
+ *
+ */
+
+namespace corsika::environment {
+
+ template <class T>
+ class FlatExponential : public T {
+ corsika::units::si::MassDensityType const fRho0;
+ units::si::LengthType const fLambda;
+ units::si::InverseLengthType const fInvLambda;
+ NuclearComposition const fNuclComp;
+ geometry::Vector<units::si::dimensionless_d> const fAxis;
+ geometry::Point const fP0;
+
+ public:
+ FlatExponential(geometry::Point const& p0,
+ geometry::Vector<units::si::dimensionless_d> const& axis,
+ units::si::MassDensityType rho, units::si::LengthType lambda,
+ NuclearComposition pNuclComp)
+ : fRho0(rho)
+ , fLambda(lambda)
+ , fInvLambda(1 / lambda)
+ , fNuclComp(pNuclComp)
+ , fAxis(axis)
+ , fP0(p0) {}
+
+ corsika::units::si::MassDensityType GetMassDensity(
+ corsika::geometry::Point const& p) const override {
+ return fRho0 * exp(fInvLambda * (p - fP0).dot(fAxis));
+ }
+ NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+
+ corsika::units::si::GrammageType IntegratedGrammage(
+ corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::units::si::LengthType pTo) const override {
+ auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
+
+ if (vDotA == 0) {
+ return pTo * GetMassDensity(line.GetR0());
+ } else {
+ return GetMassDensity(line.GetR0()) * (fLambda / vDotA) *
+ (exp(vDotA * pTo / fLambda) - 1);
+ }
+ }
+
+ corsika::units::si::LengthType ArclengthFromGrammage(
+ corsika::geometry::Trajectory<corsika::geometry::Line> const& line,
+ corsika::units::si::GrammageType pGrammage) const override {
+ auto const vDotA = line.NormalizedDirection().dot(fAxis).magnitude();
+
+ if (vDotA == 0) {
+ return pGrammage / GetMassDensity(line.GetR0());
+ } else {
+ return fLambda / vDotA * log(pGrammage * vDotA / (fRho0 * fLambda) + 1);
+ }
+ }
+ };
+
+} // namespace corsika::environment
+#endif
diff --git a/Environment/testEnvironment.cc b/Environment/testEnvironment.cc
index 5bb18bfe09e48d4865f5aafcd8ad4aa01f4afeb8..5beab6587a20199a24bb8946b96ccedfe4b8d293 100644
--- a/Environment/testEnvironment.cc
+++ b/Environment/testEnvironment.cc
@@ -13,6 +13,7 @@
// cpp file
#include <corsika/environment/DensityFunction.h>
+#include <corsika/environment/FlatExponential.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/IMediumModel.h>
#include <corsika/environment/InhomogeneousMedium.h>
@@ -30,6 +31,12 @@ using namespace corsika::geometry;
using namespace corsika::environment;
using namespace corsika::particles;
using namespace corsika::units::si;
+using namespace corsika;
+
+CoordinateSystem const& gCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+Point const gOrigin(gCS, {0_m, 0_m, 0_m});
TEST_CASE("HomogeneousMedium") {
NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
@@ -37,6 +44,51 @@ TEST_CASE("HomogeneousMedium") {
HomogeneousMedium<IMediumModel> const medium(19.2_g / cube(1_cm), protonComposition);
}
+TEST_CASE("FlatExponential") {
+ NuclearComposition const protonComposition(std::vector<Code>{Code::Proton},
+ std::vector<float>{1.f});
+
+ Vector const axis(gCS, QuantityVector<dimensionless_d>(0, 0, 1));
+ LengthType const lambda = 3_m;
+ auto const rho0 = 1_g / units::si::detail::static_pow<3>(1_cm);
+ FlatExponential<IMediumModel> const medium(gOrigin, axis, rho0, lambda,
+ protonComposition);
+ auto const tEnd = 5_s;
+
+ SECTION("horizontal") {
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {20_cm / second, 0_m / second, 0_m / second}));
+ Trajectory<Line> const trajectory(line, tEnd);
+
+ REQUIRE((medium.IntegratedGrammage(trajectory, 2_m) / (rho0 * 2_m)) == Approx(1));
+ REQUIRE((medium.ArclengthFromGrammage(trajectory, rho0 * 5_m) / 5_m) == Approx(1));
+ }
+
+ SECTION("vertical") {
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {0_m / second, 0_m / second, 5_m / second}));
+ Trajectory<Line> const trajectory(line, tEnd);
+ LengthType const length = 2 * lambda;
+ GrammageType const exact = rho0 * lambda * (exp(length / lambda) - 1);
+
+ REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
+ REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
+ }
+
+ SECTION("vertical") {
+ Line const line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {0_m / second, 5_m / second, 5_m / second}));
+ Trajectory<Line> const 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;
+
+ REQUIRE((medium.IntegratedGrammage(trajectory, length) / exact) == Approx(1));
+ REQUIRE((medium.ArclengthFromGrammage(trajectory, exact) / length) == Approx(1));
+ }
+}
+
auto constexpr rho0 = 1_kg / 1_m / 1_m / 1_m;
struct Exponential {
@@ -60,15 +112,10 @@ struct Exponential {
};
TEST_CASE("InhomogeneousMedium") {
- CoordinateSystem const& cs =
- RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
- Point const origin(cs, {0_m, 0_m, 0_m});
-
- Vector direction(cs, QuantityVector<dimensionless_d>(1, 0, 0));
+ Vector direction(gCS, QuantityVector<dimensionless_d>(1, 0, 0));
- Line line(origin, Vector<SpeedType::dimension_type>(
- cs, {20_m / second, 0_m / second, 0_m / second}));
+ Line line(gOrigin, Vector<SpeedType::dimension_type>(
+ gCS, {20_m / second, 0_m / second, 0_m / second}));
auto const tEnd = 5_s;
Trajectory<Line> const trajectory(line, tEnd);
@@ -77,9 +124,9 @@ TEST_CASE("InhomogeneousMedium") {
DensityFunction<decltype(e), LinearApproximationIntegrator> const rho(e);
SECTION("DensityFunction") {
- REQUIRE(e.Derivative<1>(origin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
+ REQUIRE(e.Derivative<1>(gOrigin, direction) / (1_kg / 1_m / 1_m / 1_m / 1_m) ==
Approx(1));
- REQUIRE(rho.EvaluateAt(origin) == e(origin));
+ REQUIRE(rho.EvaluateAt(gOrigin) == e(gOrigin));
}
auto const exactGrammage = [](auto l) { return 1_m * rho0 * (exp(l / 1_m) - 1); };
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index 737b7830821df05cdf0fa5fd93918741f0cf731e..ca97abb50feec457527751d79cd9e76b760cb4f2 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -125,17 +125,9 @@ namespace corsika::geometry {
template <typename ScalarDim>
auto operator*(phys::units::quantity<ScalarDim, double> const p) const {
- using ProdQuantity = phys::units::detail::Product<dim, ScalarDim, double, double>;
+ using ProdDim = phys::units::detail::product_d<dim, ScalarDim>;
- if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
- // not a "Quantity" anymore
- {
- return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs,
- BaseVector<dim>::qVector * p);
- } else {
- return Vector<typename ProdQuantity::dimension_type>(
- *BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
- }
+ return Vector<ProdDim>(*BaseVector<dim>::cs, BaseVector<dim>::qVector * p);
}
template <typename ScalarDim>
@@ -173,16 +165,8 @@ namespace corsika::geometry {
auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
auto const bareResult = c1.cross(c2);
- using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
-
- if constexpr (std::is_same<ProdQuantity, double>::value) // result dimensionless,
- // not a "Quantity" anymore
- {
- return Vector<phys::units::dimensionless_d>(*BaseVector<dim>::cs, bareResult);
- } else {
- return Vector<typename ProdQuantity::dimension_type>(*BaseVector<dim>::cs,
- bareResult);
- }
+ using ProdDim = phys::units::detail::product_d<dim, dim2>;
+ return Vector<ProdDim>(*BaseVector<dim>::cs, bareResult);
}
template <typename dim2>
@@ -191,9 +175,10 @@ namespace corsika::geometry {
auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
auto const bareResult = c1.dot(c2);
- using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+ using ProdDim = phys::units::detail::product_d<dim, dim2>;
- return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
+ return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
+ bareResult);
}
};
diff --git a/ThirdParty/phys/units/quantity.hpp b/ThirdParty/phys/units/quantity.hpp
index bec55318685ffb0e69bfeeb19260ef20e6d1eed1..06fbb4b2c760dfc8641d58392aa30062ff7c4e84 100644
--- a/ThirdParty/phys/units/quantity.hpp
+++ b/ThirdParty/phys/units/quantity.hpp
@@ -144,9 +144,18 @@ namespace phys {
d8 = DX::dim8 + DY::dim8
};
- typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+ typedef Collapse<dim, T> type;
};
+ /**
+ * convenience type for product dimension
+ */
+
+ template <typename DX, typename DY>
+ using product_d = typename product<DX, DY, Rep>::dim;
+
template <typename DX, typename DY, typename X, typename Y>
using Product = typename product<DX, DY, PromoteMul<X, Y>>::type;
@@ -166,9 +175,18 @@ namespace phys {
d8 = DX::dim8 - DY::dim8
};
- typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+ typedef Collapse<dim, T> type;
};
+ /**
+ * convenience type for quotient dimension
+ */
+
+ template <typename DX, typename DY>
+ using quotient_d = typename quotient<DX, DY, Rep>::dim;
+
template <typename DX, typename DY, typename X, typename Y>
using Quotient = typename quotient<DX, DY, PromoteMul<X, Y>>::type;
@@ -188,9 +206,18 @@ namespace phys {
d8 = -D::dim8
};
- typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+ typedef Collapse<dim, T> type;
};
+ /**
+ * convenience type for reciprocal dimension
+ */
+
+ template <typename DX, typename DY>
+ using reciprocal_d = typename reciprocal<DX, Rep>::dim;
+
template <typename D, typename X, typename Y>
using Reciprocal = typename reciprocal<D, PromoteMul<X, Y>>::type;
@@ -210,11 +237,20 @@ namespace phys {
d8 = N * D::dim8
};
- typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+ typedef Collapse<dim, T> type;
};
+ /**
+ * convenience type for power dimension
+ */
+
+ template <typename DX, int N>
+ using power_d = typename power<DX, N, Rep>::dim;
+
template <typename D, int N, typename T>
- using Power = typename detail::power<D, N, T>::type;
+ using Power = typename power<D, N, T>::type;
/**
* root type generator.
@@ -238,9 +274,17 @@ namespace phys {
d8 = D::dim8 / N
};
- typedef Collapse<dimensions<d1, d2, d3, d4, d5, d6, d7, d8>, T> type;
+ typedef dimensions<d1, d2, d3, d4, d5, d6, d7, d8> dim;
+
+ typedef Collapse<dim, T> type;
};
+ /**
+ * convenience type for root dimension
+ */
+ template <typename D, int N>
+ using root_d = typename root<D, N, Rep>::dim;
+
template <typename D, int N, typename T>
using Root = typename detail::root<D, N, T>::type;