IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 1dbb9245 authored by Remy Prechelt's avatar Remy Prechelt Committed by ralfulrich
Browse files

Port over exponential atmosphere profiles.

parent 676c72e6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 211 additions and 150 deletions
corsika/detail/media/BaseExponential.inl
+ 35
− 22
View file @ 1dbb9245
...@@ -10,49 +10,62 @@ ...@@ -10,49 +10,62 @@
#pragma once #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 { namespace corsika {
template <class TDerived> template <typename TDerived>
auto const& BaseExponential<TDerived>::GetImplementation() const { auto const& BaseExponential<TDerived>::getImplementation() const {
return *static_cast<TDerived const*>(this); return *static_cast<TDerived const*>(this);
} }
template <class TDerived> template <typename TDerived>
GrammageType BaseExponential<TDerived>::IntegratedGrammage( units::si::GrammageType BaseExponential<TDerived>::integratedGrammage(
Trajectory<Line> const& vLine, LengthType vL, Trajectory<Line> const& line, units::si::LengthType vL,
Vector<dimensionless_d> const& vAxis) const { Vector<units::si::dimensionless_d> const& axis) const {
if (vL == LengthType::zero()) { return GrammageType::zero(); } if (vL == units::si::LengthType::zero()) { return units::si::GrammageType::zero(); }
auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude(); auto const uDotA = line.NormalizedDirection().dot(axis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0()); auto const rhoStart = getImplementation().getMassDensity(line.GetR0());
if (uDotA == 0) { if (uDotA == 0) {
return vL * rhoStart; return vL * rhoStart;
} else { } else {
return rhoStart * (fLambda / uDotA) * (exp(uDotA * vL * fInvLambda) - 1); return rhoStart * (lambda_ / uDotA) * (exp(uDotA * vL * invLambda_) - 1);
} }
} }
template <class TDerived> template <typename TDerived>
LengthType BaseExponential<TDerived>::ArclengthFromGrammage( units::si::LengthType BaseExponential<TDerived>::arclengthFromGrammage(
Trajectory<Line> const& vLine, GrammageType vGrammage, Trajectory<Line> const& line, units::si::GrammageType grammage,
Vector<dimensionless_d> const& vAxis) const { Vector<units::si::dimensionless_d> const& axis) const {
auto const uDotA = vLine.NormalizedDirection().dot(vAxis).magnitude(); auto const uDotA = line.NormalizedDirection().dot(axis).magnitude();
auto const rhoStart = GetImplementation().GetMassDensity(vLine.GetR0()); auto const rhoStart = getImplementation().getMassDensity(line.GetR0());
if (uDotA == 0) { if (uDotA == 0) {
return vGrammage / rhoStart; return grammage / rhoStart;
} else { } else {
auto const logArg = vGrammage * fInvLambda * uDotA / rhoStart + 1; auto const logArg = grammage * invLambda_ * uDotA / rhoStart + 1;
if (logArg > 0) { if (logArg > 0) {
return fLambda / uDotA * log(logArg); return lambda_ / uDotA * log(logArg);
} else { } else {
return std::numeric_limits<typename decltype(vGrammage)::value_type>::infinity() * return std::numeric_limits<typename decltype(grammage)::value_type>::infinity() *
meter; 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 } // namespace corsika
corsika/detail/media/FlatExponential.inl
+ 35
− 16
View file @ 1dbb9245
/* /*
* (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. * See file AUTHORS for a list of contributors.
* *
...@@ -10,30 +10,49 @@ ...@@ -10,30 +10,49 @@
#pragma once #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 { namespace corsika {
template <class T> template <typename T>
MassDensityType FlatExponential<T>::GetMassDensity(Point const& vP) const { FlatExponential<T>::FlatExponential(Point const& point,
return Base::fRho0 * exp(Base::fInvLambda * (vP - Base::fP0).dot(fAxis)); 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> template <typename T>
NuclearComposition const& FlatExponential<T>::GetNuclearComposition() const { NuclearComposition const& FlatExponential<T>::getNuclearComposition() const {
return fNuclComp; return nuclComp_;
} }
template <class T> template <typename T>
GrammageType FlatExponential<T>::IntegratedGrammage(Trajectory<Line> const& vLine, units::si::GrammageType FlatExponential<T>::integratedGrammage(
LengthType vTo) const { Trajectory<Line> const& line, units::si::LengthType to) const {
return Base::IntegratedGrammage(vLine, vTo, fAxis); return BaseExponential<FlatExponential<T>>::integratedGrammage(line, to, axis_);
} }
template <class T> template <typename T>
LengthType FlatExponential<T>::ArclengthFromGrammage(Trajectory<Line> const& vLine, units::si::LengthType FlatExponential<T>::arclengthFromGrammage(
GrammageType vGrammage) const { Trajectory<Line> const& line, units::si::GrammageType grammage) const {
return Base::ArclengthFromGrammage(vLine, vGrammage, fAxis); return BaseExponential<FlatExponential<T>>::arclengthFromGrammage(line, grammage,
axis_);
} }
} // namespace corsika } // namespace corsika
......
corsika/detail/media/LinearApproximationIntegrator.inl
+ 36
− 40
View file @ 1dbb9245
...@@ -14,43 +14,39 @@ ...@@ -14,43 +14,39 @@
namespace corsika { namespace corsika {
template <typename TDerived>
template <class TDerived> auto const& LinearApproximationIntegrator<TDerived>::getImplementation() const {
auto const& LinearApproximationIntegrator<TDerived>::GetImplementation() const { return *static_cast<TDerived const*>(this);
return *static_cast<TDerived const*>(this); }
}
template <typename TDerived>
auto LinearApproximationIntegrator<TDerived>::integrateGrammage(
template <class TDerived> Trajectory<Line> const& line, units::si::LengthType length) const {
auto LinearApproximationIntegrator<TDerived>::IntegrateGrammage( auto const c0 = getImplementation().evaluateAt(line.GetPosition(0));
corsika::Trajectory<corsika::Line> const& line, auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0),
LengthType length) const { line.NormalizedDirection());
auto const c0 = GetImplementation().EvaluateAt(line.GetPosition(0)); return (c0 + 0.5 * c1 * length) * length;
auto const c1 = GetImplementation().fRho.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 {
template <class TDerived> auto const c0 = getImplementation().rho_(line.GetPosition(0));
auto LinearApproximationIntegrator<TDerived>::ArclengthFromGrammage( auto const c1 = getImplementation().rho_.FirstDerivative(line.GetPosition(0),
corsika::Trajectory<corsika::Line> const& line, line.NormalizedDirection());
GrammageType grammage) const {
auto const c0 = GetImplementation().fRho(line.GetPosition(0)); return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0;
auto const c1 = GetImplementation().fRho.FirstDerivative( }
line.GetPosition(0), line.NormalizedDirection());
template <typename TDerived>
return (1 - 0.5 * grammage * c1 / (c0 * c0)) * grammage / c0; auto LinearApproximationIntegrator<TDerived>::maximumLength(
} Trajectory<Line> const& line, [[maybe_unused]] double relError) const {
using namespace units::si;
template <class TDerived> [[maybe_unused]] auto const c1 = getImplementation().rho_.SecondDerivative(
auto LinearApproximationIntegrator<TDerived>::MaximumLength(corsika::Trajectory<corsika::Line> const& line, line.GetPosition(0), line.NormalizedDirection());
[[maybe_unused]] double relError) const {
[[maybe_unused]] auto const c1 = GetImplementation().fRho.SecondDerivative( // todo: provide a real, working implementation
line.GetPosition(0), line.NormalizedDirection()); return 1_m * std::numeric_limits<double>::infinity();
}
// todo: provide a real, working implementation
return 1_m * std::numeric_limits<double>::infinity(); } // namespace corsika
}
}
\ No newline at end of file
corsika/detail/media/SlidingPlanarExponential.inl
+ 36
− 15
View file @ 1dbb9245
...@@ -14,24 +14,45 @@ ...@@ -14,24 +14,45 @@
namespace corsika { namespace corsika {
template <class T> template <typename T>
MassDensityType SlidingPlanarExponential<T>::GetMassDensity(Point const& p) const { SlidingPlanarExponential<T>::SlidingPlanarExponential(
auto const height = (p - Base::fP0).norm() - referenceHeight_; Point const& p0, units::si::MassDensityType rho0, units::si::LengthType lambda,
return Base::fRho0 * exp(Base::fInvLambda * height); 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> template <typename T>
GrammageType SlidingPlanarExponential<T>::IntegratedGrammage( units::si::GrammageType SlidingPlanarExponential<T>::integratedGrammage(
Trajectory<Line> const& line, LengthType l) const { Trajectory<Line> const& line, units::si::LengthType l) const {
auto const axis = (line.GetR0() - Base::fP0).normalized(); auto const axis =
return Base::IntegratedGrammage(line, l, axis); (line.GetR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
return BaseExponential<SlidingPlanarExponential<T>>::integratedGrammage(line, l,
axis);
} }
template <class T> template <typename T>
LengthType SlidingPlanarExponential<T>::ArclengthFromGrammage( units::si::LengthType SlidingPlanarExponential<T>::arclengthFromGrammage(
Trajectory<Line> const& line, GrammageType grammage) const { Trajectory<Line> const& line, units::si::GrammageType const grammage) const {
auto const axis = (line.GetR0() - Base::fP0).normalized(); auto const axis =
return Base::ArclengthFromGrammage(line, grammage, axis); (line.GetR0() - BaseExponential<SlidingPlanarExponential<T>>::point_).normalized();
return BaseExponential<SlidingPlanarExponential<T>>::arclengthFromGrammage(
line, grammage, axis);
} }
} // namespace corsika } // namespace corsika
\ No newline at end of file
corsika/media/BaseExponential.hpp
+ 22
− 24
View file @ 1dbb9245
...@@ -22,23 +22,15 @@ namespace corsika { ...@@ -22,23 +22,15 @@ namespace corsika {
* This class provides the grammage/length conversion functionality for * This class provides the grammage/length conversion functionality for
* (locally) flat exponential atmospheres. * (locally) flat exponential atmospheres.
*/ */
template <class TDerived> template <typename TDerived>
class BaseExponential { class BaseExponential {
public:
BaseExponential(Point const& vP0, MassDensityType vRho, LengthType vLambda)
: fRho0(vRho)
, fLambda(vLambda)
, fInvLambda(1 / vLambda)
, fP0(vP0) {}
protected: protected:
MassDensityType const fRho0; units::si::MassDensityType const rho0_;
LengthType const fLambda; units::si::LengthType const lambda_;
InverseLengthType const fInvLambda; units::si::InverseLengthType const invLambda_;
Point const fP0; Point const point_;
auto const& GetImplementation() const; auto const& getImplementation() const;
// clang-format off // clang-format off
/** /**
...@@ -47,15 +39,16 @@ namespace corsika { ...@@ -47,15 +39,16 @@ namespace corsika {
* \f[ * \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) * 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. * \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: * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[ * \f[
* X = \varrho_0 l; * X = \varrho_0 l;
* \f] * \f]
*/ */
// clang-format on // clang-format on
GrammageType IntegratedGrammage(Trajectory<Line> const& vLine, LengthType vL, units::si::GrammageType integratedGrammage(
Vector<dimensionless_d> const& vAxis) const; Trajectory<Line> const& line, units::si::LengthType vL,
Vector<units::si::dimensionless_d> const& axis) const;
// clang-format off // clang-format off
/** /**
...@@ -64,22 +57,27 @@ namespace corsika { ...@@ -64,22 +57,27 @@ namespace corsika {
* \f[ * \f[
* l = \begin{cases} * l = \begin{cases}
* \frac{\lambda}{\vec{u} \cdot \vec{a}} \log\left(Y \right), & \text{if} Y := 0 > 1 + * \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,} * \infty & \text{else,}
* \end{cases} * \end{cases}
* \f] where \f$ \varrho_0 \f$ is the density at the starting point. * \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: * If \f$ \vec{u} \cdot \vec{a} = 0 \f$, the calculation is just like with a homogeneous density:
* \f[ * \f[
* l = \frac{X}{\varrho_0} * l = \frac{X}{\varrho_0}
* \f] * \f]
*/ */
// clang-format on // clang-format on
LengthType ArclengthFromGrammage(Trajectory<Line> const& vLine, units::si::LengthType arclengthFromGrammage(
GrammageType vGrammage, Trajectory<Line> const& line, units::si::GrammageType grammage,
Vector<dimensionless_d> const& vAxis) const; 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> #include <corsika/detail/media/BaseExponential.inl>
corsika/media/FlatExponential.hpp
+ 23
− 16
View file @ 1dbb9245
...@@ -8,7 +8,6 @@ n/* ...@@ -8,7 +8,6 @@ n/*
#pragma once #pragma once
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp> #include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/Line.hpp> #include <corsika/framework/geometry/Line.hpp>
#include <corsika/framework/geometry/Point.hpp> #include <corsika/framework/geometry/Point.hpp>
...@@ -18,30 +17,38 @@ n/* ...@@ -18,30 +17,38 @@ n/*
namespace corsika { 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 { class FlatExponential : public BaseExponential<FlatExponential<T>>, public T {
Vector<dimensionless_d> const fAxis; Vector<units::si::dimensionless_d> const axis_;
NuclearComposition const fNuclComp; NuclearComposition const nuclComp_;
using Base = BaseExponential<FlatExponential<T>>; using Base = BaseExponential<FlatExponential<T>>;
public: public:
FlatExponential(Point const& vP0, Vector<dimensionless_d> const& vAxis, FlatExponential(Point const& point, Vector<units::si::dimensionless_d> const& axis,
MassDensityType vRho, LengthType vLambda, units::si::MassDensityType rho, units::si::LengthType lambda,
NuclearComposition vNuclComp) NuclearComposition nuclComp);
: Base(vP0, vRho, vLambda)
, fAxis(vAxis)
, fNuclComp(vNuclComp) {}
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, units::si::GrammageType integratedGrammage(Trajectory<Line> const& line,
LengthType vTo) const override; units::si::LengthType to) const;
LengthType ArclengthFromGrammage(Trajectory<Line> const& vLine, units::si::LengthType arclengthFromGrammage(Trajectory<Line> const& line,
GrammageType vGrammage) const override; units::si::GrammageType grammage) const;
}; };
} // namespace corsika } // namespace corsika
......
corsika/media/SlidingPlanarExponential.hpp
+ 24
− 17
View file @ 1dbb9245
...@@ -19,34 +19,41 @@ ...@@ -19,34 +19,41 @@
namespace corsika { 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>>, class SlidingPlanarExponential : public BaseExponential<SlidingPlanarExponential<T>>,
public T { public T {
NuclearComposition const nuclComp_; NuclearComposition const nuclComp_;
LengthType const referenceHeight_; units::si::LengthType const referenceHeight_;
using Base = BaseExponential<SlidingPlanarExponential<T>>; using Base = BaseExponential<SlidingPlanarExponential<T>>;
public: public:
SlidingPlanarExponential(Point const& p0, MassDensityType rho0, LengthType lambda, SlidingPlanarExponential(
NuclearComposition nuclComp, Point const& p0, units::si::MassDensityType rho0, units::si::LengthType lambda,
LengthType referenceHeight = LengthType::zero()) NuclearComposition nuclComp,
: Base(p0, rho0, lambda) units::si::LengthType referenceHeight = units::si::LengthType::zero());
, nuclComp_(nuclComp)
, referenceHeight_(referenceHeight) {}
inline MassDensityType GetMassDensity(Point const& p) const override; units::si::MassDensityType getMassDensity(Point const& point) const override;
inline NuclearComposition const& GetNuclearComposition() const override { NuclearComposition const& getNuclearComposition() const override;
return nuclComp_;
}
inline GrammageType IntegratedGrammage(Trajectory<Line> const& line, units::si::GrammageType integratedGrammage(Trajectory<Line> const& line,
LengthType l) const override; units::si::LengthType l) const override;
inline LengthType ArclengthFromGrammage(Trajectory<Line> const& line, units::si::LengthType arclengthFromGrammage(
GrammageType grammage) const override; Trajectory<Line> const& line, units::si::GrammageType grammage) const override;
}; };
} // namespace corsika } // namespace corsika
......
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