From 1dbb9245457cadffc91e88c1fea4d73824220e9d Mon Sep 17 00:00:00 2001
From: Remy Prechelt <prechelt@hawaii.edu>
Date: Mon, 16 Nov 2020 23:18:11 -1000
Subject: [PATCH] Port over exponential atmosphere profiles.
---
corsika/detail/media/BaseExponential.inl | 57 ++++++++------
corsika/detail/media/FlatExponential.inl | 51 +++++++++----
.../media/LinearApproximationIntegrator.inl | 76 +++++++++----------
.../detail/media/SlidingPlanarExponential.inl | 51 +++++++++----
corsika/media/BaseExponential.hpp | 46 ++++++-----
corsika/media/FlatExponential.hpp | 39 ++++++----
corsika/media/SlidingPlanarExponential.hpp | 41 +++++-----
7 files changed, 211 insertions(+), 150 deletions(-)
diff --git a/corsika/detail/media/BaseExponential.inl b/corsika/detail/media/BaseExponential.inl
index 4e954264e..7cfac790e 100644
--- a/corsika/detail/media/BaseExponential.inl
+++ b/corsika/detail/media/BaseExponential.inl
@@ -10,49 +10,62 @@
#pragma once
-#include <corsika/media/BaseExponential.hpp>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Trajectory.hpp>
namespace corsika {
- template <class TDerived>
- auto const& BaseExponential<TDerived>::GetImplementation() const {
+ template <typename TDerived>
+ auto const& BaseExponential<TDerived>::getImplementation() const {
return *static_cast<TDerived const*>(this);
}
- template <class TDerived>
- GrammageType BaseExponential<TDerived>::IntegratedGrammage(
- Trajectory<Line> const& vLine, LengthType vL,
- Vector<dimensionless_d> const& vAxis) const {
- if (vL == LengthType::zero()) { return GrammageType::zero(); }
+ template <typename TDerived>
+ units::si::GrammageType BaseExponential<TDerived>::integratedGrammage(
+ Trajectory<Line> const& line, units::si::LengthType vL,
+ Vector<units::si::dimensionless_d> const& axis) const {
+ if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
- auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
- auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+ auto const uDotA = line.NormalizedDirection().dot(axis).magnitude();
+ auto const rhoStart = getImplementation().getMassDensity(line.GetR0());
if (uDotA == 0) {
return vL * rhoStart;
} else {
- return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1);
+ return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1);
}
}
- template <class TDerived>
- LengthType BaseExponential<TDerived>::ArclengthFromGrammage(
- Trajectory<Line> const& vLine, GrammageType vGrammage,
- Vector<dimensionless_d> const& vAxis) const {
- auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude();
- auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0());
+ template <typename TDerived>
+ units::si::LengthType BaseExponential<TDerived>::arclengthFromGrammage(
+ Trajectory<Line> const& line, units::si::GrammageType grammage,
+ Vector<units::si::dimensionless_d> const& axis) const {
+ auto const uDotA = line.NormalizedDirection().dot(axis).magnitude();
+ auto const rhoStart = getImplementation().getMassDensity(line.GetR0());
if (uDotA == 0) {
- return vGrammage / rhoStart;
+ return grammage / rhoStart;
} else {
- auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1;
+ auto const logArg = grammage * invLambda_ * uDotA / rhoStart + 1;
if (logArg > 0) {
- return fLambda / uDotA * log(logArg);
+ return lambda_ / uDotA * log(logArg);
} else {
- return std::numeric_limits<typename decltype(vGrammage)::value_type>::infinity() *
- meter;
+ return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
+ units::si::meter;
}
}
}
+ template <typename TDerived>
+ BaseExponential<TDerived>::BaseExponential(Point const& point,
+ units::si::MassDensityType rho0,
+ units::si::LengthType lambda)
+ : rho0_(rho0)
+ , lambda_(lambda)
+ , invLambda_(1 / lambda)
+ , point_(point) {}
+
} // namespace corsika
diff --git a/corsika/detail/media/FlatExponential.inl b/corsika/detail/media/FlatExponential.inl
index 31453d7b9..34b94e0a2 100644
--- a/corsika/detail/media/FlatExponential.inl
+++ b/corsika/detail/media/FlatExponential.inl
@@ -1,5 +1,5 @@
/*
- * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
*
* See file AUTHORS for a list of contributors.
*
@@ -10,30 +10,49 @@
#pragma once
-#include <corsika/media/FlatExponential.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Trajectory.hpp>
+#include <corsika/media/BaseExponential.hpp>
+#include <corsika/media/NuclearComposition.hpp>
namespace corsika {
- template <class T>
- MassDensityType FlatExponential<T>::GetMassDensity(Point const& vP) const {
- return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis));
+ template <typename T>
+ FlatExponential<T>::FlatExponential(Point const& point,
+ Vector<units::si::dimensionless_d> const& axis,
+ units::si::MassDensityType rho,
+ units::si::LengthType lambda,
+ NuclearComposition nuclComp)
+ : BaseExponential<FlatExponential<T>>(point, rho, lambda)
+ , axis_(axis)
+ , nuclComp_(nuclComp) {}
+
+ template <typename T>
+ units::si::MassDensityType FlatExponential<T>::getMassDensity(
+ Point const& point) const {
+ return BaseExponential<FlatExponential<T>>::rho0_ *
+ exp(BaseExponential<FlatExponential<T>>::invLambda_ *
+ (point - BaseExponential<FlatExponential<T>>::point_).dot(axis_));
}
- template <class T>
- NuclearComposition const& FlatExponential<T>::GetNuclearComposition() const {
- return fNuclComp;
+ template <typename T>
+ NuclearComposition const& FlatExponential<T>::getNuclearComposition() const {
+ return nuclComp_;
}
- template <class T>
- GrammageType FlatExponential<T>::IntegratedGrammage(Trajectory<Line> const& vLine,
- LengthType vTo) const {
- return Base::IntegratedGrammage(vLine, vTo, fAxis);
+ template <typename T>
+ units::si::GrammageType FlatExponential<T>::integratedGrammage(
+ Trajectory<Line> const& line, units::si::LengthType to) const {
+ return BaseExponential<FlatExponential<T>>::integratedGrammage(line, to, axis_);
}
- template <class T>
- LengthType FlatExponential<T>::ArclengthFromGrammage(Trajectory<Line> const& vLine,
- GrammageType vGrammage) const {
- return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis);
+ template <typename T>
+ units::si::LengthType FlatExponential<T>::arclengthFromGrammage(
+ Trajectory<Line> const& line, units::si::GrammageType grammage) const {
+ return BaseExponential<FlatExponential<T>>::arclengthFromGrammage(line, grammage,
+ axis_);
}
} // namespace corsika
diff --git a/corsika/detail/media/LinearApproximationIntegrator.inl b/corsika/detail/media/LinearApproximationIntegrator.inl
index 324bdbb96..98d52db1c 100644
--- a/corsika/detail/media/LinearApproximationIntegrator.inl
+++ b/corsika/detail/media/LinearApproximationIntegrator.inl
@@ -14,43 +14,39 @@
namespace corsika {
-
- template <class TDerived>
- auto const& LinearApproximationIntegrator<TDerived>::GetImplementation() const {
- return *static_cast<TDerived const*>(this);
- }
-
-
- template <class TDerived>
- auto LinearApproximationIntegrator<TDerived>::IntegrateGrammage(
- corsika::Trajectory<corsika::Line> const& line,
- LengthType length) const {
- 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;
- }
-
- template <class TDerived>
- auto LinearApproximationIntegrator<TDerived>::ArclengthFromGrammage(
- corsika::Trajectory<corsika::Line> const& line,
- GrammageType grammage) const {
- auto const c0 = GetImplementation().fRho(line.GetPosition(0));
- auto const c1 = GetImplementation().fRho.FirstDerivative(
- line.GetPosition(0), line.NormalizedDirection());
-
- return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
- }
-
- template <class TDerived>
- auto LinearApproximationIntegrator<TDerived>::MaximumLength(corsika::Trajectory<corsika::Line> const& line,
- [[maybe_unused]] double relError) const {
- [[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative(
- line.GetPosition(0), line.NormalizedDirection());
-
- // todo: provide a real, working implementation
- return 1_m * std::numeric_limits<double>::infinity();
- }
-
-
-}
\ No newline at end of file
+ template <typename TDerived>
+ auto const& LinearApproximationIntegrator<TDerived>::getImplementation() const {
+ return *static_cast<TDerived const*>(this);
+ }
+
+ template <typename TDerived>
+ auto LinearApproximationIntegrator<TDerived>::integrateGrammage(
+ Trajectory<Line> const& line, units::si::LengthType length) const {
+ auto const c0 = getImplementation().evaluateAt(line.GetPosition(0));
+ auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0),
+ line.NormalizedDirection());
+ return (c0 + 0.5 * c1 * length) * length;
+ }
+
+ template <typename TDerived>
+ auto LinearApproximationIntegrator<TDerived>::arclengthFromGrammage(
+ Trajectory<Line> const& line, units::si::GrammageType grammage) const {
+ auto const c0 = getImplementation().rho_(line.GetPosition(0));
+ auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0),
+ line.NormalizedDirection());
+
+ return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
+ }
+
+ template <typename TDerived>
+ auto LinearApproximationIntegrator<TDerived>::maximumLength(
+ Trajectory<Line> const& line, [[maybe_unused]] double relError) const {
+ using namespace units::si;
+ [[maybe_unused]] auto const c1 = getImplementation().rho_.SecondDerivative(
+ line.GetPosition(0), line.NormalizedDirection());
+
+ // todo: provide a real, working implementation
+ return 1_m * std::numeric_limits<double>::infinity();
+ }
+
+} // namespace corsika
diff --git a/corsika/detail/media/SlidingPlanarExponential.inl b/corsika/detail/media/SlidingPlanarExponential.inl
index 61043f072..991820e1c 100644
--- a/corsika/detail/media/SlidingPlanarExponential.inl
+++ b/corsika/detail/media/SlidingPlanarExponential.inl
@@ -14,24 +14,45 @@
namespace corsika {
- template <class T>
- MassDensityType SlidingPlanarExponential<T>::GetMassDensity(Point const& p) const {
- auto const height = (p - Base::fP0).norm() - referenceHeight_;
- return Base::fRho0 * exp(Base::fInvLambda * height);
+ template <typename T>
+ SlidingPlanarExponential<T>::SlidingPlanarExponential(
+ Point const& p0, units::si::MassDensityType rho0, units::si::LengthType lambda,
+ NuclearComposition nuclComp, units::si::LengthType referenceHeight)
+ : BaseExponential<SlidingPlanarExponential<T>>(p0, rho0, lambda)
+ , nuclComp_(nuclComp)
+ , referenceHeight_(referenceHeight) {}
+
+ template <typename T>
+ units::si::MassDensityType SlidingPlanarExponential<T>::getMassDensity(
+ Point const& point) const {
+ auto const height =
+ (point - BaseExponential<SlidingPlanarExponential<T>>::point_).norm() -
+ referenceHeight_;
+ return BaseExponential<SlidingPlanarExponential<T>>::rho0_ *
+ exp(BaseExponential<SlidingPlanarExponential<T>>::invLambda_ * height);
+ }
+
+ template <typename T>
+ NuclearComposition const& SlidingPlanarExponential<T>::getNuclearComposition() const {
+ return nuclComp_;
}
- template <class T>
- GrammageType SlidingPlanarExponential<T>::IntegratedGrammage(
- Trajectory<Line> const& line, LengthType l) const {
- auto const axis = (line.GetR0() - Base::fP0).normalized();
- return Base::IntegratedGrammage(line, l, axis);
+ template <typename T>
+ units::si::GrammageType SlidingPlanarExponential<T>::integratedGrammage(
+ Trajectory<Line> const& line, units::si::LengthType l) const {
+ auto const axis =
+ (line.GetR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
+ return BaseExponential<SlidingPlanarExponential<T>>::integratedGrammage(line, l,
+ axis);
}
- template <class T>
- LengthType SlidingPlanarExponential<T>::ArclengthFromGrammage(
- Trajectory<Line> const& line, GrammageType grammage) const {
- auto const axis = (line.GetR0() - Base::fP0).normalized();
- return Base::ArclengthFromGrammage(line, grammage, axis);
+ template <typename T>
+ units::si::LengthType SlidingPlanarExponential<T>::arclengthFromGrammage(
+ Trajectory<Line> const& line, units::si::GrammageType const grammage) const {
+ auto const axis =
+ (line.GetR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
+ return BaseExponential<SlidingPlanarExponential<T>>::arclengthFromGrammage(
+ line, grammage, axis);
}
-} // namespace corsika
\ No newline at end of file
+} // namespace corsika
diff --git a/corsika/media/BaseExponential.hpp b/corsika/media/BaseExponential.hpp
index 577629db2..2a85a0b4f 100644
--- a/corsika/media/BaseExponential.hpp
+++ b/corsika/media/BaseExponential.hpp
@@ -22,23 +22,15 @@ namespace corsika {
* This class provides the grammage/length conversion functionality for
* (locally) flat exponential atmospheres.
*/
- template <class TDerived>
+ template <typename TDerived>
class BaseExponential {
-
- public:
- BaseExponential(Point const& vP0, MassDensityType vRho, LengthType vLambda)
- : fRho0(vRho)
- , fLambda(vLambda)
- , fInvLambda(1 / vLambda)
- , fP0(vP0) {}
-
protected:
- MassDensityType const fRho0;
- LengthType const fLambda;
- InverseLengthType const fInvLambda;
- Point const fP0;
+ units::si::MassDensityType const rho0_;
+ units::si::LengthType const lambda_;
+ units::si::InverseLengthType const invLambda_;
+ Point const point_;
- auto const& GetImplementation() const;
+ auto const& getImplementation() const;
// clang-format off
/**
@@ -47,15 +39,16 @@ namespace corsika {
* \f[
* X = \frac{\varrho_0 \lambda}{\vec{u} \cdot \vec{a}} \left( \exp\left( \vec{u} \cdot \vec{a} \frac{l}{\lambda} \right) - 1 \right)
* \f], where \f$ \varrho_0 \f$ is the density at the starting point.
- *
+ *
* If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[
* X = \varrho_0 l;
* \f]
*/
// clang-format on
- GrammageType IntegratedGrammage(Trajectory<Line> const& vLine, LengthType vL,
- Vector<dimensionless_d> const& vAxis) const;
+ units::si::GrammageType integratedGrammage(
+ Trajectory<Line> const& line, units::si::LengthType vL,
+ Vector<units::si::dimensionless_d> const& axis) const;
// clang-format off
/**
@@ -64,22 +57,27 @@ namespace corsika {
* \f[
* l = \begin{cases}
* \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 +
- * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda}
+ * \vec{u} \cdot \vec{a} \frac{X}{\rho_0 \lambda}
* \infty & \text{else,}
* \end{cases}
* \f] where \f$ \varrho_0 \f$ is the density at the starting point.
- *
+ *
* If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[
* l = \frac{X}{\varrho_0}
* \f]
*/
// clang-format on
- LengthType ArclengthFromGrammage(Trajectory<Line> const& vLine,
- GrammageType vGrammage,
- Vector<dimensionless_d> const& vAxis) const;
- };
+ units::si::LengthType arclengthFromGrammage(
+ Trajectory<Line> const& line, units::si::GrammageType grammage,
+ Vector<units::si::dimensionless_d> const& axis) const;
+
+ public:
+ BaseExponential(Point const& point, units::si::MassDensityType rho0,
+ units::si::LengthType lambda);
+
+ }; // class BaseExponential
-} // namespace corsika
+ } // namespace corsika
#include <corsika/detail/media/BaseExponential.inl>
diff --git a/corsika/media/FlatExponential.hpp b/corsika/media/FlatExponential.hpp
index c0757dfc1..acb456d12 100644
--- a/corsika/media/FlatExponential.hpp
+++ b/corsika/media/FlatExponential.hpp
@@ -8,7 +8,6 @@ n/*
#pragma once
-#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp>
@@ -18,30 +17,38 @@ n/*
namespace corsika {
- template <class T>
+ // clang-format off
+ /**
+ * flat exponential density distribution with
+ * \f[
+ * \varrho(r) = \varrho_0 \exp\left( \frac{1}{\lambda} (r - p) \cdot
+ * \vec{a} \right).
+ * \f]
+ * \f$ \vec{a} \f$ denotes the axis and should be normalized to avoid degeneracy
+ * with the scale parameter \f$ \lambda \f$.
+ */
+ // clang-format on
+ template <typename T>
class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
- Vector<dimensionless_d> const fAxis;
- NuclearComposition const fNuclComp;
+ Vector<units::si::dimensionless_d> const axis_;
+ NuclearComposition const nuclComp_;
using Base = BaseExponential<FlatExponential<T>>;
public:
- FlatExponential(Point const& vP0, Vector<dimensionless_d> const& vAxis,
- MassDensityType vRho, LengthType vLambda,
- NuclearComposition vNuclComp)
- : Base(vP0, vRho, vLambda)
- , fAxis(vAxis)
- , fNuclComp(vNuclComp) {}
+ FlatExponential(Point const& point, Vector<units::si::dimensionless_d> const& axis,
+ units::si::MassDensityType rho, units::si::LengthType lambda,
+ NuclearComposition nuclComp);
- MassDensityType GetMassDensity(Point const& vP) const override;
+ units::si::MassDensityType getMassDensity(Point const& point) const override;
- NuclearComposition const& GetNuclearComposition() const override;
+ NuclearComposition const& getNuclearComposition() const override;
- GrammageType IntegratedGrammage(Trajectory<Line> const& vLine,
- LengthType vTo) const override;
+ units::si::GrammageType integratedGrammage(Trajectory<Line> const& line,
+ units::si::LengthType to) const;
- LengthType ArclengthFromGrammage(Trajectory<Line> const& vLine,
- GrammageType vGrammage) const override;
+ units::si::LengthType arclengthFromGrammage(Trajectory<Line> const& line,
+ units::si::GrammageType grammage) const;
};
} // namespace corsika
diff --git a/corsika/media/SlidingPlanarExponential.hpp b/corsika/media/SlidingPlanarExponential.hpp
index c34949722..013135491 100644
--- a/corsika/media/SlidingPlanarExponential.hpp
+++ b/corsika/media/SlidingPlanarExponential.hpp
@@ -19,34 +19,41 @@
namespace corsika {
- template <class T>
+ // clang-format off
+ /**
+ * The SlidingPlanarExponential models mass density as
+ * \f[
+ * \varrho(r) = \varrho_0 \exp\left( \frac{|p_0 - r|}{\lambda} \right).
+ * \f]
+ * For grammage/length conversion, the density distribution is approximated as
+ * locally flat at the starting point \f$ r_0 \f$ of the trajectory with the axis pointing
+ * from \f$ p_0 \f$ to \f$ r_0 \f$.
+ */
+ // clang-format on
+
+ template <typename T>
class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
public T {
-
NuclearComposition const nuclComp_;
- LengthType const referenceHeight_;
+ units::si::LengthType const referenceHeight_;
using Base = BaseExponential<SlidingPlanarExponential<T>>;
public:
- SlidingPlanarExponential(Point const& p0, MassDensityType rho0, LengthType lambda,
- NuclearComposition nuclComp,
- LengthType referenceHeight = LengthType::zero())
- : Base(p0, rho0, lambda)
- , nuclComp_(nuclComp)
- , referenceHeight_(referenceHeight) {}
+ SlidingPlanarExponential(
+ Point const& p0, units::si::MassDensityType rho0, units::si::LengthType lambda,
+ NuclearComposition nuclComp,
+ units::si::LengthType referenceHeight = units::si::LengthType::zero());
- inline MassDensityType GetMassDensity(Point const& p) const override;
+ units::si::MassDensityType getMassDensity(Point const& point) const override;
- inline NuclearComposition const& GetNuclearComposition() const override {
- return nuclComp_;
- }
+ NuclearComposition const& getNuclearComposition() const override;
- inline GrammageType IntegratedGrammage(Trajectory<Line> const& line,
- LengthType l) const override;
+ units::si::GrammageType integratedGrammage(Trajectory<Line> const& line,
+ units::si::LengthType l) const override;
- inline LengthType ArclengthFromGrammage(Trajectory<Line> const& line,
- GrammageType grammage) const override;
+ units::si::LengthType arclengthFromGrammage(
+ Trajectory<Line> const& line, units::si::GrammageType grammage) const override;
};
} // namespace corsika
--
GitLab