diff --git a/corsika/detail/media/BaseTabular.inl b/corsika/detail/media/BaseTabular.inl new file mode 100644 index 0000000000000000000000000000000000000000..916fb03202c17aee1b2e4f105399bcbb63786bb8 --- /dev/null +++ b/corsika/detail/media/BaseTabular.inl @@ -0,0 +1,201 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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. + */ + +#pragma once + +#include <corsika/framework/core/ParticleProperties.hpp> +#include <corsika/framework/core/PhysicalUnits.hpp> +#include <corsika/framework/geometry/Point.hpp> + +#include <exception> + +namespace corsika { + + template <typename TDerived> + inline BaseTabular<TDerived>::BaseTabular( + Point const& point, LengthType const referenceHeight, + std::function<MassDensityType(LengthType)> const& rho, unsigned int const nBins, + LengthType const deltaHeight) + : nBins_(nBins) + , deltaHeight_(deltaHeight) + , point_(point) + , referenceHeight_(referenceHeight) { + density_.resize(nBins_); + for (unsigned int bin = 0; bin < nBins; ++bin) { + density_[bin] = rho(deltaHeight_ * bin); + CORSIKA_LOG_DEBUG("new tabulated atm bin={} h={} rho={}", bin, deltaHeight_ * bin, + density_[bin]); + } + } + + template <typename TDerived> + inline auto const& BaseTabular<TDerived>::getImplementation() const { + return *static_cast<TDerived const*>(this); + } + + template <typename TDerived> + inline MassDensityType BaseTabular<TDerived>::getMassDensity( + LengthType const height) const { + double const fbin = (height - referenceHeight_) / deltaHeight_; + int const bin = int(fbin); + if (bin < 0) return MassDensityType::zero(); + if (bin >= int(nBins_ - 1)) { + CORSIKA_LOG_ERROR( + "invalid height {} (corrected {}) in BaseTabular atmosphere. Min 0, max {}. If " + "max is too low: increase!", + height, height - referenceHeight_, nBins_ * deltaHeight_); + throw std::runtime_error("invalid height"); + } + return density_[bin] + (fbin - bin) * (density_[bin + 1] - density_[bin]); + } + + template <typename TDerived> + inline GrammageType BaseTabular<TDerived>::getIntegratedGrammage( + BaseTrajectory const& traj) const { + + Point pCurr = traj.getPosition(0); + DirectionVector dCurr = traj.getDirection(0); + LengthType height1 = (traj.getPosition(0) - point_).getNorm() - referenceHeight_; + LengthType height2 = (traj.getPosition(1) - point_).getNorm() - referenceHeight_; + + LengthType const fullLength = traj.getLength(1); + + int sign = +1; // normal direction + if (height1 > height2) { + std::swap(height1, height2); + pCurr = traj.getPosition(1); + dCurr = traj.getDirection(1); + sign = -1; // inverted direction + } + + double const fbin1 = height1 / deltaHeight_; + unsigned int const bin1 = int(fbin1); + + double const fbin2 = height2 / deltaHeight_; + unsigned int const bin2 = int(fbin2); + + if (fbin1 == fbin2) { return GrammageType::zero(); } + + if (bin1 >= nBins_ - 1 || bin2 >= nBins_ - 1) { + CORSIKA_LOG_ERROR("invalid height {} {} in BaseTabular atmosphere integration", + height1, height2); + throw std::runtime_error("invalid height"); + } + + // interpolated start/end densities + MassDensityType const rho1 = getMassDensity(height1 + referenceHeight_); + MassDensityType const rho2 = getMassDensity(height2 + referenceHeight_); + + // within first bin + if (bin1 == bin2) { return fullLength * (rho2 + rho1) / 2; } + + // inclination of trajectory (local) + DirectionVector axis((pCurr - point_).normalized()); // to gravity center + double cosTheta = axis.dot(dCurr); + + // distance to next height bin + unsigned int bin = bin1; + LengthType dD = (deltaHeight_ * (bin + 1) - height1) / cosTheta * sign; + LengthType distance = dD; + + GrammageType X = dD * (rho1 + density_[bin + 1]) / 2; + double frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength); + pCurr = traj.getPosition(frac); + dCurr = traj.getDirection(frac); + + for (++bin; bin < bin2; ++bin) { + // inclination of trajectory + axis = (pCurr - point_).normalized(); + cosTheta = axis.dot(dCurr); + // distance to next height bin + dD = deltaHeight_ / cosTheta * sign; + distance += dD; + GrammageType const dX = dD * (density_[bin] + density_[bin + 1]) / 2; + X += dX; + frac = (sign > 0 ? distance / fullLength : 1 - distance / fullLength); + pCurr = traj.getPosition(frac); + dCurr = traj.getDirection(frac); + } + + // inclination of trajectory + axis = ((pCurr - point_).normalized()); + cosTheta = axis.dot(dCurr); + // distance to next height bin + dD = (height2 - deltaHeight_ * bin2) / cosTheta * sign; + X += dD * (rho2 + density_[bin2]) / 2; + distance += dD; + return X; + } + + template <typename TDerived> + inline LengthType BaseTabular<TDerived>::getArclengthFromGrammage( + BaseTrajectory const& traj, GrammageType const grammage) const { + + if (grammage < GrammageType::zero()) { + CORSIKA_LOG_ERROR("cannot integrate negative grammage"); + throw std::runtime_error("negative grammage error"); + } + + LengthType const height = (traj.getPosition(0) - point_).getNorm() - referenceHeight_; + + double const fbin = height / deltaHeight_; + int bin = int(fbin); + + if (bin >= int(nBins_ - 1)) { + CORSIKA_LOG_ERROR("invalid height {} in BaseTabular atmosphere integration", + height); + throw std::runtime_error("invalid height"); + } + + // interpolated start/end densities + MassDensityType const rho = getMassDensity(height + referenceHeight_); + + // inclination of trajectory + Point pCurr = traj.getPosition(0); + DirectionVector dCurr = traj.getDirection(0); + DirectionVector axis((pCurr - point_).normalized()); + double cosTheta = axis.dot(dCurr); + int sign = +1; // height increasing along traj + if (cosTheta < 0) { + cosTheta = -cosTheta; // absolute value only + sign = -1; // height decreasing along traj + } + + // height -> distance + LengthType const deltaDistance = deltaHeight_ / cosTheta; + + // start with 0 g/cm2 + GrammageType X(GrammageType::zero()); + LengthType distance(LengthType::zero()); + + // within first bin + distance = + (sign > 0 ? deltaDistance * (bin + 1 - fbin) : deltaDistance * (fbin - bin)); + GrammageType binGrammage = (sign > 0 ? distance * (rho + density_[bin + 1]) / 2 + : distance * (rho + density_[bin]) / 2); + if (X + binGrammage > grammage) { + double const binFraction = (grammage - X) / binGrammage; + return distance * binFraction; + } + X += binGrammage; + + // the following bins (along trajectory) + for (bin += sign; bin < int(nBins_) && bin >= 0; bin += sign) { + + binGrammage = deltaDistance * (density_[bin] + density_[bin + 1]) / 2; + if (X + binGrammage > grammage) { + double const binFraction = (grammage - X) / binGrammage; + return distance + deltaDistance * binFraction; + } + X += binGrammage; + distance += deltaDistance; + } + return std::numeric_limits<double>::infinity() * meter; + } + +} // namespace corsika diff --git a/corsika/detail/media/SlidingPlanarTabular.inl b/corsika/detail/media/SlidingPlanarTabular.inl new file mode 100644 index 0000000000000000000000000000000000000000..c8e15792cafea62f7ba7bb05b8388fdca2fa0608 --- /dev/null +++ b/corsika/detail/media/SlidingPlanarTabular.inl @@ -0,0 +1,47 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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. + */ + +#pragma once + +namespace corsika { + + template <typename T> + inline SlidingPlanarTabular<T>::SlidingPlanarTabular( + Point const& p0, std::function<MassDensityType(LengthType)> const& rho, + unsigned int const nBins, LengthType const deltaHeight, + NuclearComposition const& nuclComp, LengthType referenceHeight) + : BaseTabular<SlidingPlanarTabular<T>>(p0, referenceHeight, rho, nBins, deltaHeight) + , nuclComp_(nuclComp) {} + + template <typename T> + inline MassDensityType SlidingPlanarTabular<T>::getMassDensity( + Point const& point) const { + auto const heightFull = + (point - BaseTabular<SlidingPlanarTabular<T>>::getAnchorPoint()).getNorm(); + return BaseTabular<SlidingPlanarTabular<T>>::getMassDensity(heightFull); + } + + template <typename T> + inline NuclearComposition const& SlidingPlanarTabular<T>::getNuclearComposition() + const { + return nuclComp_; + } + + template <typename T> + inline GrammageType SlidingPlanarTabular<T>::getIntegratedGrammage( + BaseTrajectory const& traj) const { + return BaseTabular<SlidingPlanarTabular<T>>::getIntegratedGrammage(traj); + } + + template <typename T> + inline LengthType SlidingPlanarTabular<T>::getArclengthFromGrammage( + BaseTrajectory const& traj, GrammageType const grammage) const { + return BaseTabular<SlidingPlanarTabular<T>>::getArclengthFromGrammage(traj, grammage); + } + +} // namespace corsika diff --git a/corsika/media/BaseTabular.hpp b/corsika/media/BaseTabular.hpp new file mode 100644 index 0000000000000000000000000000000000000000..a552b66483d305f4893f5e457c3e68bcd12d53a1 --- /dev/null +++ b/corsika/media/BaseTabular.hpp @@ -0,0 +1,89 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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. + */ + +#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> +#include <corsika/framework/geometry/BaseTrajectory.hpp> + +#include <limits> +#include <vector> + +namespace corsika { + + /** + * This class provides the grammage/length conversion functionality for + * (locally) flat exponential atmospheres. + */ + template <typename TDerived> + class BaseTabular { + + public: + BaseTabular(Point const& point, LengthType const referenceHeight, + std::function<MassDensityType(LengthType)> const& rho, + unsigned int const nBins, LengthType const deltaHeight); + + Point const& getAnchorPoint() const { return point_; } + + 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) + * direction \f$ \vec{u} \f$ is given by + * \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 getIntegratedGrammage(BaseTrajectory const& line) const; + + // clang-format off + /** + * For a (normalized) axis \f$ \vec{a} \f$, the length of a non-orthogonal line with (normalized) + * direction \f$ \vec{u} \f$ corresponding to grammage \f$ X \f$ is given by + * \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} + * \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 getArclengthFromGrammage(BaseTrajectory const& line, + GrammageType const grammage) const; + + private: + unsigned int const nBins_; + LengthType deltaHeight_; + std::vector<MassDensityType> density_; + Point const point_; + LengthType const referenceHeight_; + + }; // class BaseTabular + +} // namespace corsika + +#include <corsika/detail/media/BaseTabular.inl> diff --git a/corsika/media/SlidingPlanarTabular.hpp b/corsika/media/SlidingPlanarTabular.hpp new file mode 100644 index 0000000000000000000000000000000000000000..4550792aa9fcb8e7b46f2b901c979d689c11b682 --- /dev/null +++ b/corsika/media/SlidingPlanarTabular.hpp @@ -0,0 +1,62 @@ +/* + * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu + * + * 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. + */ + +#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> +#include <corsika/framework/random/RNGManager.hpp> +#include <corsika/framework/geometry/BaseTrajectory.hpp> +#include <corsika/media/NuclearComposition.hpp> +#include <corsika/media/BaseTabular.hpp> + +namespace corsika { + + // clang-format off + /** + * The SlidingPlanarTabular 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 TDerived> + class SlidingPlanarTabular : public BaseTabular<SlidingPlanarTabular<TDerived>>, + public TDerived { + + using Base = BaseTabular<SlidingPlanarTabular<TDerived>>; + + public: + SlidingPlanarTabular(Point const& p0, + std::function<MassDensityType(LengthType)> const& rho, + unsigned int const nBins, LengthType const deltaHeight, + NuclearComposition const& nuclComp, + LengthType referenceHeight = LengthType::zero()); + + MassDensityType getMassDensity(Point const& point) const override; + + NuclearComposition const& getNuclearComposition() const override; + + GrammageType getIntegratedGrammage(BaseTrajectory const& line) const override; + + LengthType getArclengthFromGrammage(BaseTrajectory const& line, + GrammageType grammage) const override; + + private: + NuclearComposition const nuclComp_; + }; + +} // namespace corsika + +#include <corsika/detail/media/SlidingPlanarTabular.inl>