corsika/detail/framework/core/Logging.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <spdlog/fmt/ostr.h> // will output whenerver a streaming operator is found
+#include <spdlog/sinks/stdout_color_sinks.h>
+#include <spdlog/spdlog.h>
+
+namespace corsika {
+
+  // many of these free functions are special to the logging
+  // infrastructure so we hide them in the corsika::logging namespace.
+  namespace logging {
+
+    /*
+     * The default pattern for CORSIKA8 loggers.
+     */
+    inline auto set_default_level(level::level_enum const minlevel) -> void {
+      spdlog::set_level(minlevel);
+    }
+
+    template <typename TLogger>
+    inline auto add_source_info(TLogger& logger) -> void {
+      logger->set_pattern(source_pattern);
+    }
+
+    template <typename TLogger>
+    inline auto reset_pattern(TLogger& logger) -> void {
+      logger->set_pattern(default_pattern);
+    }
+
+  } // namespace logging
+
+  inline std::shared_ptr<spdlog::logger> create_logger(std::string const& name,
+                                                       bool const defaultlog) {
+
+    // create the logger
+    // this is currently a colorized multi-threading safe logger
+    auto logger = spdlog::stdout_color_mt(name);
+
+    // set the default C8 format
+#if (!defined(_GLIBCXX_USE_CXX11_ABI) || _GLIBCXX_USE_CXX11_ABI == 1)
+    logger->set_pattern(default_pattern);
+#else
+    // special case: gcc from the software collections devtoolset
+    std::string dp(default_pattern);
+    logger->set_pattern(dp);
+#endif
+
+    // if defaultlog is True, we set this as the default spdlog logger.
+    if (defaultlog) { spdlog::set_default_logger(logger); }
+
+    return logger;
+  }
+
+  inline std::shared_ptr<spdlog::logger> get_logger(std::string const& name,
+                                                    bool const defaultlog) {
+
+    // attempt to get the logger from the registry
+    auto logger = spdlog::get(name);
+
+    // weg found the logger, so just return it
+    if (logger) {
+      return logger;
+    } else { // logger was not found so create it
+      return create_logger(name, defaultlog);
+    }
+  }
+
+} // namespace corsika

corsika/detail/framework/core/ParticleProperties.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <cmath>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/core/Logging.hpp>
+
+namespace corsika {
+
+  inline HEPEnergyType get_kinetic_energy_propagation_threshold(Code const code) {
+    if (is_nucleus(code)) return particle::detail::threshold_nuclei;
+    return particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)];
+  }
+
+  inline void set_kinetic_energy_propagation_threshold(Code const code,
+                                                       HEPEnergyType const val) {
+    if (is_nucleus(code))
+      particle::detail::threshold_nuclei = val;
+    else
+      particle::detail::propagation_thresholds[static_cast<CodeIntType>(code)] = val;
+  }
+
+  inline HEPMassType constexpr get_mass(Code const code) {
+    if (is_nucleus(code)) { return get_nucleus_mass(code); }
+    return particle::detail::masses[static_cast<CodeIntType>(code)];
+  }
+
+  inline HEPEnergyType get_energy_production_threshold(Code const p) {
+    return particle::detail::production_thresholds[static_cast<CodeIntType>(p)];
+  }
+
+  inline void set_energy_production_threshold(Code const p, HEPEnergyType const val) {
+    particle::detail::production_thresholds[static_cast<CodeIntType>(p)] = val;
+  }
+
+  inline bool constexpr is_charged(Code const c) { return get_charge_number(c) != 0; }
+
+  inline bool constexpr is_nucleus(Code const code) { return code >= Code::Nucleus; }
+
+  inline Code constexpr get_nucleus_code(size_t const A,
+                                         size_t const Z) { // 10LZZZAAAI
+    if (Z > A) { throw std::runtime_error("Z cannot be larger than A in nucleus."); }
+    return static_cast<Code>(static_cast<CodeIntType>(Code::Nucleus) + Z * 10000 +
+                             A * 10);
+  }
+
+  inline size_t constexpr get_nucleus_Z(Code const code) {
+    return (static_cast<CodeIntType>(code) % static_cast<CodeIntType>(Code::Nucleus)) /
+           10000;
+  }
+
+  inline size_t constexpr get_nucleus_A(Code const code) {
+    return (static_cast<CodeIntType>(code) % 10000) / 10;
+  }
+
+  inline PDGCode constexpr get_PDG(Code const code) {
+    if (code < Code::Nucleus) {
+      return particle::detail::pdg_codes[static_cast<CodeIntType>(code)];
+    }
+    size_t const Z = get_nucleus_Z(code);
+    size_t const A = get_nucleus_A(code);
+    return static_cast<PDGCode>(static_cast<CodeIntType>(Code::Nucleus) + Z * 10000 +
+                                A * 10); // 10LZZZAAAI
+  }
+
+  inline PDGCode constexpr get_PDG(unsigned int const A, unsigned int const Z) {
+    return PDGCode(1000000000 + Z * 10000 + A * 10);
+  }
+
+  inline int16_t constexpr get_charge_number(Code const code) {
+    if (is_nucleus(code)) return get_nucleus_Z(code);
+    return particle::detail::electric_charges[static_cast<CodeIntType>(code)];
+  }
+
+  inline ElectricChargeType constexpr get_charge(Code const code) {
+    if (code == Code::Nucleus)
+      throw std::runtime_error("charge of particle::Nucleus undefined");
+    return get_charge_number(code) * constants::e;
+  }
+
+  inline std::string_view constexpr get_name(Code const code) {
+    if (is_nucleus(code)) { return "nucleus"; }
+    return particle::detail::names[static_cast<CodeIntType>(code)];
+  }
+
+  inline std::string get_name(Code code, full_name) {
+    if (is_nucleus(code)) {
+      return fmt::format("nucleus ({},{})", get_nucleus_A(code), get_nucleus_Z(code));
+    }
+
+    return std::string{get_name(code)};
+  }
+
+  inline TimeType constexpr get_lifetime(Code const p) {
+    return particle::detail::lifetime[static_cast<CodeIntType>(p)] * second;
+  }
+
+  inline bool constexpr is_hadron(Code const code) {
+    if (is_nucleus(code)) return true;
+    return particle::detail::isHadron[static_cast<CodeIntType>(code)];
+  }
+
+  inline bool constexpr is_em(Code const c) {
+    return c == Code::Electron || c == Code::Positron || c == Code::Photon;
+  }
+
+  inline bool constexpr is_muon(Code const c) {
+    return c == Code::MuPlus || c == Code::MuMinus;
+  }
+
+  inline bool constexpr is_neutrino(Code const c) {
+    return c == Code::NuE || c == Code::NuMu || c == Code::NuTau || c == Code::NuEBar ||
+           c == Code::NuMuBar || c == Code::NuTauBar;
+  }
+
+  inline std::ostream& operator<<(std::ostream& stream, corsika::Code const code) {
+    return stream << get_name(code);
+  }
+
+  inline Code convert_from_PDG(PDGCode const p) {
+    static_assert(particle::detail::conversionArray.size() % 2 == 1);
+    // this will fail, for the strange case where the maxPDG is negative...
+    int constexpr maxPDG{(particle::detail::conversionArray.size() - 1) >> 1};
+    auto const k = static_cast<PDGCodeIntType>(p);
+    if (std::abs(k) <= maxPDG) {
+      return particle::detail::conversionArray[k + maxPDG];
+    } else {
+      if (1000000000 <= k && k <= 1009999990) { // nucleus (no L or I)
+        int const Z = (k - 1000000000) / 10000;
+        int const A = (k - 1000000000 - 10000 * Z) / 10;
+        return get_nucleus_code(A, Z);
+      }
+      return particle::detail::conversionMap.at(p);
+    }
+  }
+
+  inline HEPMassType constexpr get_nucleus_mass(Code const code) {
+    unsigned int const A = get_nucleus_A(code);
+    unsigned int const Z = get_nucleus_Z(code);
+    return get_nucleus_mass(A, Z);
+  }
+
+  inline HEPMassType constexpr get_nucleus_mass(unsigned int const A,
+                                                unsigned int const Z) {
+    return get_mass(Code::Proton) * Z + (A - Z) * get_mass(Code::Neutron);
+  }
+
+  inline std::string get_nucleus_name(Code const code) {
+    size_t const A = get_nucleus_A(code);
+    size_t const Z = get_nucleus_Z(code);
+    return fmt::format("Nucleus_A{}_Z{}", A, Z);
+  }
+
+  inline std::initializer_list<Code> constexpr get_all_particles() {
+    return particle::detail::all_particles;
+  }
+
+} // namespace corsika

corsika/detail/framework/core/SpdlogSpecializations.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <fmt/format.h>
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <boost/filesystem/path.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <Eigen/Dense>
+
+//-----------------------------------
+// STD
+//-----------------------------------
+namespace std {
+
+  auto inline format_as(std::_Put_time<char> const& arg) {
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+} // namespace std
+
+//-----------------------------------
+// CORSIKA
+//-----------------------------------
+namespace corsika {
+
+  // formatters for particle codes declared on ParticleProperties.hpp
+  auto inline format_as(Code code) { return get_name(code); }
+  auto inline format_as(PDGCode code) { return fmt::underlying(code); }
+
+  template <typename Type>
+  auto inline format_as(Type const& arg) {
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+} // namespace corsika
+
+//-----------------------------------
+// boost::filesystem
+//-----------------------------------
+namespace boost::filesystem {
+
+  auto inline format_as(path const& fname) { return fname.string().c_str(); }
+
+} // namespace boost::filesystem
+
+//-----------------------------------
+// phys::units
+//-----------------------------------
+namespace phys::units {
+
+  template <typename Dimensions>
+  auto inline format_as(phys::units::quantity<Dimensions> const& arg) {
+
+    return io::to_string(arg);
+  }
+
+} // namespace phys::units
+
+//----------------------------------
+// Eigen
+//----------------------------------
+namespace Eigen {
+
+  template <typename Scalar, int M, int N>
+  auto inline format_as(Matrix<Scalar, M, N> const& arg) {
+
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+  template <typename Matrix1, typename Matrix2>
+  auto inline format_as(Product<Matrix1, Matrix2> const& arg) {
+
+    std::ostringstream os;
+    os << arg;
+    return os.str();
+  }
+
+} // namespace Eigen

corsika/detail/framework/geometry/BaseVector.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/CoordinateSystem.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+
+namespace corsika {
+
+  template <typename TDimension>
+  inline CoordinateSystemPtr BaseVector<TDimension>::getCoordinateSystem() const {
+    return cs_;
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> const& BaseVector<TDimension>::getQuantityVector()
+      const {
+    return quantityVector_;
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension>& BaseVector<TDimension>::getQuantityVector() {
+    return quantityVector_;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Box.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+  inline Box::Box(CoordinateSystemPtr cs, LengthType const x, LengthType const y,
+                  LengthType const z)
+      : center_(Point(cs, {0_m, 0_m, 0_m}))
+      , cs_(cs)
+      , x_(x)
+      , y_(y)
+      , z_(z) {}
+
+  inline Box::Box(CoordinateSystemPtr cs, LengthType const side)
+      : center_(Point(cs, {0_m, 0_m, 0_m}))
+      , cs_(cs)
+      , x_(side / 2)
+      , y_(side / 2)
+      , z_(side / 2) {}
+
+  inline bool Box::contains(Point const& p) const {
+    if ((abs(p.getX(cs_)) < x_) && (abs(p.getY(cs_)) < y_) && (abs(p.getZ(cs_)) < z_))
+      return true;
+    else
+      return false;
+  }
+
+  inline std::string Box::asString() const {
+    std::ostringstream txt;
+    txt << "center=" << center_ << ", x-axis=" << DirectionVector{cs_, {1, 0, 0}}
+        << ", y-axis: " << DirectionVector{cs_, {0, 1, 0}}
+        << ", z-axis: " << DirectionVector{cs_, {0, 0, 1}};
+    return txt.str();
+  }
+
+  template <typename TDim>
+  inline void Box::rotate(QuantityVector<TDim> const& axis, double const angle) {
+    cs_ = make_rotation(cs_, axis, angle);
+  }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/CoordinateSystem.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <Eigen/Dense>
+
+#include <memory>
+#include <stdexcept>
+
+namespace corsika {
+
+  inline CoordinateSystemPtr CoordinateSystem::getReferenceCS() const {
+    return referenceCS_;
+  }
+
+  inline EigenTransform const& CoordinateSystem::getTransform() const { return transf_; }
+
+  inline bool CoordinateSystem::operator==(CoordinateSystem const& cs) const {
+    return referenceCS_ == cs.referenceCS_ && transf_.matrix() == cs.transf_.matrix();
+  }
+
+  inline bool CoordinateSystem::operator!=(CoordinateSystem const& cs) const {
+    return !(cs == *this);
+  }
+
+  /// find transformation between two CS, using most optimal common base
+  inline EigenTransform get_transformation(CoordinateSystem const& pFrom,
+                                           CoordinateSystem const& pTo) {
+    CoordinateSystem const* a{&pFrom};
+    CoordinateSystem const* b{&pTo};
+
+    while (a != b && b) {
+
+      // traverse pFrom
+      a = &pFrom;
+      while (a != b && a) { a = a->getReferenceCS().get(); }
+
+      if (a == b) break;
+
+      b = b->getReferenceCS().get();
+    }
+
+    if (a != b || a == nullptr) {
+      throw std::runtime_error("no connection between coordinate systems found!");
+    }
+
+    CoordinateSystem const* commonBase = a;
+    CoordinateSystem const* p = &pFrom;
+    EigenTransform t = EigenTransform::Identity();
+    while ((*p) != (*commonBase)) {
+      t = p->getTransform() * t;
+      p = p->getReferenceCS().get();
+    }
+
+    p = &pTo;
+
+    while (*p != *commonBase) {
+      t = t * p->getTransform().inverse(Eigen::TransformTraits::Isometry);
+      p = p->getReferenceCS().get();
+    }
+
+    return t;
+  }
+
+  inline CoordinateSystemPtr make_translation(CoordinateSystemPtr const& cs,
+                                              QuantityVector<length_d> const& vector) {
+    EigenTransform const translation{EigenTranslation(vector.getEigenVector())};
+    return CoordinateSystemPtr{new CoordinateSystem(cs, translation)};
+  }
+
+  template <typename TDim>
+  inline CoordinateSystemPtr make_rotationToZ(CoordinateSystemPtr const& cs,
+                                              Vector<TDim> const& vVec) {
+    auto const vVecComp = vVec.getComponents(cs);
+    if (vVecComp.getX().magnitude() == 0 && vVecComp.getY().magnitude() == 0 &&
+        vVecComp.getZ().magnitude() == 0) {
+      return cs;
+    }
+
+    auto const a = vVecComp.normalized().getEigenVector();
+    auto const a1 = a(0), a2 = a(1), a3 = a(2);
+
+    Eigen::Matrix3d A, B;
+
+    if (a3 > 0) {
+      auto const c = 1 / (1 + a3);
+      A << 1, 0, a1,                      // comment to prevent clang-format
+          0, 1, a2,                       // .
+          -a1, -a2, 1;                    // .
+      B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+          -a1 * a2 * c, -a2 * a2 * c, 0,  // .
+          0, 0, -(a1 * a1 + a2 * a2) * c; // .
+
+    } else {
+      auto const c = 1 / (1 - a3);
+      A << 1, 0, a1,                      // .
+          0, -1, a2,                      // .
+          a1, -a2, -1;                    // .
+      B << -a1 * a1 * c, +a1 * a2 * c, 0, // .
+          -a1 * a2 * c, +a2 * a2 * c, 0,  // .
+          0, 0, (a1 * a1 + a2 * a2) * c;  // .
+    }
+
+    return CoordinateSystemPtr{new CoordinateSystem{cs, EigenTransform{A + B}}};
+  }
+
+  template <typename TDim>
+  inline CoordinateSystemPtr make_rotation(CoordinateSystemPtr const& cs,
+                                           QuantityVector<TDim> const& axis,
+                                           double const angle) {
+    if (axis.getEigenVector().isZero()) {
+      throw std::runtime_error("null-vector given as axis parameter");
+    }
+
+    EigenTransform const rotation{
+        Eigen::AngleAxisd(angle, axis.getEigenVector().normalized())};
+
+    return CoordinateSystemPtr{new CoordinateSystem{cs, rotation}};
+  }
+
+  template <typename TDim>
+  inline CoordinateSystemPtr make_translationAndRotation(
+      CoordinateSystemPtr const& cs, QuantityVector<length_d> const& translation,
+      QuantityVector<TDim> const& axis, double const angle) {
+    if (axis.getEigenVector().isZero()) {
+      throw std::runtime_error("null-vector given as axis parameter");
+    }
+
+    EigenTransform const transf{
+        Eigen::AngleAxisd(angle, axis.getEigenVector().normalized()) *
+        EigenTranslation(translation.getEigenVector())};
+
+    return CoordinateSystemPtr{new CoordinateSystem{cs, transf}};
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/FourVector.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <type_traits>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline TTimeType FourVector<TTimeType, TSpaceVecType>::getTimeLikeComponent() const {
+    return timeLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline TSpaceVecType& FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents() {
+    return spaceLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline TSpaceVecType const&
+  FourVector<TTimeType, TSpaceVecType>::getSpaceLikeComponents() const {
+    return spaceLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_square_type
+  FourVector<TTimeType, TSpaceVecType>::getNormSqr() const {
+    return getTimeSquared() - spaceLike_.getSquaredNorm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
+  FourVector<TTimeType, TSpaceVecType>::getNorm() const {
+
+    return sqrt(abs(getNormSqr()));
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline bool FourVector<TTimeType, TSpaceVecType>::isTimelike() const {
+    return getTimeSquared() < spaceLike_.getSquaredNorm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline bool FourVector<TTimeType, TSpaceVecType>::isSpacelike() const {
+    return getTimeSquared() > spaceLike_.getSquaredNorm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator+=(FourVector const& b) {
+    timeLike_ += b.timeLike_;
+    spaceLike_ += b.spaceLike_;
+
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator-=(FourVector const& b) {
+    timeLike_ -= b.timeLike_;
+    spaceLike_ -= b.spaceLike_;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator*=(double const b) {
+    timeLike_ *= b;
+    spaceLike_ *= b;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator/=(double const b) {
+    timeLike_ /= b;
+    spaceLike_.getComponents() /= b;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline FourVector<TTimeType, TSpaceVecType>& FourVector<TTimeType, TSpaceVecType>::
+  operator/(double const b) {
+    *this /= b;
+    return *this;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_type
+      FourVector<TTimeType, TSpaceVecType>::operator*(FourVector const& b) {
+    if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
+                                                   second)>::value)
+      return timeLike_ * b.timeLike_ * constants::cSquared - spaceLike_.norm();
+    else
+      return timeLike_ * timeLike_ - spaceLike_.norm();
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline typename FourVector<TTimeType, TSpaceVecType>::norm_square_type
+  FourVector<TTimeType, TSpaceVecType>::getTimeSquared() const {
+    if constexpr (std::is_same<time_type, decltype(std::declval<space_type>() / meter *
+                                                   second)>::value)
+      return timeLike_ * timeLike_ * constants::cSquared;
+    else
+      return timeLike_ * timeLike_;
+  }
+
+  template <typename TTimeType, typename TSpaceVecType>
+  inline std::ostream& operator<<(
+      std::ostream& os, corsika::FourVector<TTimeType, TSpaceVecType> const qv) {
+
+    os << '(' << qv.timeLike_ << ", " << qv.spaceLike_ << ") ";
+    return os;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Helix.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+#include <cmath>
+
+namespace corsika {
+
+  inline LengthType Helix::getRadius() const { return radius_; }
+
+  inline Point Helix::getPosition(TimeType const t) const {
+    return r0_ + vPar_ * t +
+           (vPerp_ * (std::cos(omegaC_ * t) - 1) + uPerp_ * std::sin(omegaC_ * t)) /
+               omegaC_;
+  }
+
+  inline Point Helix::getPositionFromArclength(LengthType const l) const {
+    return getPosition(getTimeFromArclength(l));
+  }
+
+  inline LengthType Helix::getArcLength(TimeType const t1, TimeType const t2) const {
+    return (vPar_ + vPerp_).getNorm() * (t2 - t1);
+  }
+
+  inline TimeType Helix::getTimeFromArclength(LengthType const l) const {
+    return l / (vPar_ + vPerp_).getNorm();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/LeapFrogTrajectory.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/QuarticSolver.hpp>
+
+namespace corsika {
+
+  inline Line LeapFrogTrajectory::getLine() const {
+    auto D = getPosition(1) - getPosition(0);
+    auto d = D.getNorm();
+    auto v = initialVelocity_;
+    if (d > 1_um) { // if trajectory is ultra-short, we do not
+                    // re-calculate velocity, just use initial
+                    // value. Otherwise, this is numerically unstable
+      v = D / d * getVelocity(0).getNorm();
+    }
+    return Line(getPosition(0), v);
+  }
+
+  inline Point LeapFrogTrajectory::getPosition(double const u) const {
+    if (u == 0) return initialPosition_;
+    Point const position = initialPosition_ + initialVelocity_ * timeStep_ * u / 2;
+    VelocityVector const velocity =
+        initialVelocity_ + initialVelocity_.cross(magneticfield_) * timeStep_ * u * k_;
+    return position + velocity * timeStep_ * u / 2;
+  }
+
+  inline VelocityVector LeapFrogTrajectory::getVelocity(double const u) const {
+    return getDirection(u) * initialVelocity_.getNorm();
+  }
+
+  inline DirectionVector LeapFrogTrajectory::getDirection(double const u) const {
+    if (u == 0) return initialDirection_;
+
+    return (initialDirection_ +
+            initialDirection_.cross(magneticfield_) * timeStep_ * u * k_)
+        .normalized();
+  }
+
+  inline TimeType LeapFrogTrajectory::getDuration(double const u) const {
+    TimeType const step = timeStep_ * u;
+    double const correction = 1;
+    // the eventual (delta-L to L) correction factor is:
+    //    (initialDirection_ + initialDirection_.cross(magneticfield_) * step *
+    //    k_).getNorm();
+    return step / 2 * (correction + 1);
+  }
+
+  template <typename Particle>
+  inline TimeType LeapFrogTrajectory::getTime(Particle const& particle,
+                                              double const u) const {
+    return particle.getTime() + getDuration(u);
+  }
+
+  inline LengthType LeapFrogTrajectory::getLength(double const u) const {
+    return getDuration(u) * initialVelocity_.getNorm();
+  }
+
+  inline void LeapFrogTrajectory::setLength(LengthType const limit) {
+    if (initialVelocity_.getNorm() == SpeedType::zero()) setDuration(0_s);
+    setDuration(limit / initialVelocity_.getNorm());
+  }
+
+  inline void LeapFrogTrajectory::setDuration(TimeType const limit) {
+    /*
+    initial attempt to calculate delta-L from assumed full-leap-frog-length L:
+
+    Note: often return 0. Not good enough yet.
+
+    LengthType const L = initialVelocity_.getNorm() * limit; // distance
+    double const a = (initialVelocity_.cross(magneticfield_) * k_).getSquaredNorm() / 4 /
+                     square(1_m) * static_pow<4>(1_s);
+    double const d = L * initialVelocity_.getNorm() / square(1_m) * 1_s;
+    double const e = -square(L) / square(1_m);
+    std::vector<double> solutions = solve_quartic_real(a, 0, 0, d, e);
+    CORSIKA_LOG_DEBUG("setDuration limit={} L={} solution={}", limit, L,
+                      fmt::join(solutions, ", "));
+    */
+    timeStep_ = limit;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Line.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  inline Point Line::getPosition(TimeType const t) const {
+    return start_point_ + velocity_ * t;
+  }
+
+  inline VelocityVector const& Line::getVelocity(TimeType const) const {
+    return velocity_;
+  }
+
+  inline Point Line::getPositionFromArclength(LengthType const l) const {
+    return start_point_ + velocity_.normalized() * l;
+  }
+
+  inline LengthType Line::getArcLength(TimeType const t1, TimeType const t2) const {
+    return velocity_.getNorm() * (t2 - t1);
+  }
+
+  inline TimeType Line::getTimeFromArclength(LengthType const t) const {
+    return t / velocity_.getNorm();
+  }
+
+  inline Point const& Line::getStartPoint() const { return start_point_; }
+
+  inline DirectionVector Line::getDirection() const { return velocity_.normalized(); }
+
+  inline VelocityVector const& Line::getVelocity() const { return velocity_; }
+
+} // namespace corsika

corsika/detail/framework/geometry/Path.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <deque>
+
+namespace corsika {
+
+  inline Path::Path(Point const& point) { points_.push_front(point); }
+
+  inline Path::Path(std::deque<Point> const& points)
+      : points_(points) {
+    int dequesize_ = points.size();
+    if (dequesize_ == 0 || dequesize_ == 1) {
+      length_ = LengthType::zero();
+    } else if (dequesize_ == 2) {
+      length_ = (points.back() - points.front()).getNorm();
+    } else {
+      for (auto point = points.begin(); point != points.end() - 1; ++point) {
+        auto point_next = *(point + 1);
+        auto point_now = *(point);
+        length_ += (point_next - point_now).getNorm();
+      }
+    }
+  }
+
+  inline void Path::addToEnd(Point const& point) {
+    length_ += (point - points_.back()).getNorm();
+    points_.push_back(point);
+  }
+
+  inline void Path::removeFromEnd() {
+    auto lastpoint_ = points_.back();
+    points_.pop_back();
+    int dequesize_ = points_.size();
+    if (dequesize_ == 0 || dequesize_ == 1) {
+      length_ = LengthType::zero();
+    } else if (dequesize_ == 2) {
+      length_ = (points_.back() - points_.front()).getNorm();
+    } else {
+      length_ -= (lastpoint_ - points_.back()).getNorm();
+    }
+  }
+
+  inline LengthType Path::getLength() const { return length_; }
+
+  inline Point const& Path::getStart() const { return points_.front(); }
+
+  inline Point const& Path::getEnd() const { return points_.back(); }
+
+  inline Point const& Path::getPoint(std::size_t const index) const {
+    return points_.at(index);
+  }
+
+  inline Path::const_iterator Path::begin() const { return points_.cbegin(); }
+
+  inline Path::const_iterator Path::end() const { return points_.cend(); }
+
+  inline Path::iterator Path::begin() { return points_.begin(); }
+
+  inline Path::iterator Path::end() { return points_.end(); }
+
+  inline int Path::getNSegments() const { return points_.size() - 1; }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/Plane.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+#include <sstream>
+
+namespace corsika {
+
+  inline bool Plane::isAbove(Point const& vP) const {
+    return normal_.dot(vP - center_) > LengthType::zero();
+  }
+
+  inline LengthType Plane::getDistanceTo(Point const& vP) const {
+    return (normal_ * (vP - center_).dot(normal_)).getNorm();
+  }
+
+  inline Point const& Plane::getCenter() const { return center_; }
+
+  inline DirectionVector const& Plane::getNormal() const { return normal_; }
+
+  inline std::string Plane::asString() const {
+    std::ostringstream txt;
+    txt << "center=" << center_ << ", normal=" << normal_;
+    return txt.str();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Point.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/BaseVector.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Vector.hpp>
+
+namespace corsika {
+
+  inline QuantityVector<length_d> const& Point::getCoordinates() const {
+    return BaseVector<length_d>::getQuantityVector();
+  }
+
+  inline QuantityVector<length_d>& Point::getCoordinates() {
+    return BaseVector<length_d>::getQuantityVector();
+  }
+
+  inline LengthType Point::getX(CoordinateSystemPtr const& pCS) const {
+    return getCoordinates(pCS).getX();
+  }
+
+  inline LengthType Point::getY(CoordinateSystemPtr const& pCS) const {
+    return getCoordinates(pCS).getY();
+  }
+
+  inline LengthType Point::getZ(CoordinateSystemPtr const& pCS) const {
+    return getCoordinates(pCS).getZ();
+  }
+
+  /// this always returns a QuantityVector as triple
+  inline QuantityVector<length_d> Point::getCoordinates(
+      CoordinateSystemPtr const& pCS) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    if (pCS == cs) {
+      return BaseVector<length_d>::getQuantityVector();
+    } else {
+      return QuantityVector<length_d>(
+          get_transformation(*cs.get(), *pCS.get()) *
+          BaseVector<length_d>::getQuantityVector().eigenVector_);
+    }
+  }
+
+  /// this always returns a QuantityVector as triple
+  inline QuantityVector<length_d>& Point::getCoordinates(CoordinateSystemPtr const& pCS) {
+    if (pCS != BaseVector<length_d>::getCoordinateSystem()) { rebase(pCS); }
+    return BaseVector<length_d>::getQuantityVector();
+  }
+
+  inline void Point::rebase(CoordinateSystemPtr const& pCS) {
+    BaseVector<length_d>::setQuantityVector(QuantityVector<length_d>(
+        get_transformation(*BaseVector<length_d>::getCoordinateSystem().get(),
+                           *pCS.get()) *
+        BaseVector<length_d>::getQuantityVector().eigenVector_));
+    BaseVector<length_d>::setCoordinateSystem(pCS);
+  }
+
+  inline Point Point::operator+(Vector<length_d> const& pVec) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    return Point(cs, getCoordinates() + pVec.getComponents(cs));
+  }
+
+  inline Point Point::operator-(Vector<length_d> const& pVec) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    return Point(cs, getCoordinates() - pVec.getComponents(cs));
+  }
+
+  inline Vector<length_d> Point::operator-(Point const& pB) const {
+    CoordinateSystemPtr const& cs = BaseVector<length_d>::getCoordinateSystem();
+    return Vector<length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
+  }
+  inline std::ostream& operator<<(std::ostream& os, corsika::Point const& p) {
+    auto const& qv = p.getCoordinates();
+    os << qv << " (ref:" << fmt::ptr(p.getCoordinateSystem()) << ")";
+    return os;
+  }
+
+  inline LengthType distance(Point const& p1, Point const& p2) {
+    return (p1 - p2).getNorm();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/QuantityVector.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <Eigen/Dense>
+
+#include <iostream>
+#include <utility>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type QuantityVector<TDimension>::
+  operator[](size_t const index) const {
+    return quantity_type(phys::units::detail::magnitude_tag, eigenVector_[index]);
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getX() const {
+    return (*this)[0];
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getY() const {
+    return (*this)[1];
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getZ() const {
+    return (*this)[2];
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_type
+  QuantityVector<TDimension>::getNorm() const {
+    return quantity_type(phys::units::detail::magnitude_tag, eigenVector_.norm());
+  }
+
+  template <typename TDimension>
+  inline typename QuantityVector<TDimension>::quantity_square_type
+  QuantityVector<TDimension>::getSquaredNorm() const {
+    using QuantitySquared =
+        decltype(std::declval<quantity_type>() * std::declval<quantity_type>());
+    return QuantitySquared(phys::units::detail::magnitude_tag,
+                           eigenVector_.squaredNorm());
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> QuantityVector<TDimension>::operator+(
+      QuantityVector<TDimension> const& pQVec) const {
+    return QuantityVector<TDimension>(eigenVector_ + pQVec.eigenVector_);
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> QuantityVector<TDimension>::operator-(
+      QuantityVector<TDimension> const& pQVec) const {
+    return QuantityVector<TDimension>(eigenVector_ - pQVec.eigenVector_);
+  }
+
+  template <typename TDimension>
+  template <typename TScalarDim>
+  inline auto QuantityVector<TDimension>::operator*(
+      phys::units::quantity<TScalarDim, double> const p) const {
+    using ResQuantity =
+        phys::units::detail::Product<TScalarDim, TDimension, double, double>;
+
+    if constexpr (std::is_same<ResQuantity, double>::value) // result dimensionless, not
+                                                            // a "quantity_type" anymore
+    {
+      return QuantityVector<phys::units::dimensionless_d>(eigenVector_ * p.magnitude());
+    } else {
+      return QuantityVector<typename ResQuantity::dimension_type>(eigenVector_ *
+                                                                  p.magnitude());
+    }
+  }
+
+  template <typename TDimension>
+  template <typename TScalarDim>
+  inline auto QuantityVector<TDimension>::operator/(
+      phys::units::quantity<TScalarDim, double> const p) const {
+    return (*this) * (1 / p);
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator*(double const p) const {
+    return QuantityVector<TDimension>(eigenVector_ * p);
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator/(double const p) const {
+    return QuantityVector<TDimension>(eigenVector_ / p);
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator/=(double const p) {
+    eigenVector_ /= p;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator*=(double const p) {
+    eigenVector_ *= p;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator+=(
+      QuantityVector<TDimension> const& pQVec) {
+    eigenVector_ += pQVec.eigenVector_;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& QuantityVector<TDimension>::operator-=(
+      QuantityVector<TDimension> const& pQVec) {
+    eigenVector_ -= pQVec.eigenVector_;
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator-() const {
+    return QuantityVector<TDimension>(-eigenVector_);
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::normalized() const {
+    return QuantityVector<TDimension>(eigenVector_.normalized());
+  }
+
+  template <typename TDimension>
+  inline auto QuantityVector<TDimension>::operator==(
+      QuantityVector<TDimension> const& p) const {
+    return eigenVector_ == p.eigenVector_;
+  }
+
+  template <typename TDimension>
+  inline std::ostream& operator<<(std::ostream& os,
+                                  corsika::QuantityVector<TDimension> const& qv) {
+    using quantity_type = phys::units::quantity<TDimension, double>;
+
+    os << '(' << qv.eigenVector_(0) << ' ' << qv.eigenVector_(1) << ' '
+       << qv.eigenVector_(2) << ") "
+       << phys::units::to_unit_symbol<TDimension, double>(
+              quantity_type(phys::units::detail::magnitude_tag, 1));
+    return os;
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/SeparationPlane.inl

+/*
+ * (c) Copyright 2023 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+namespace corsika {
+
+  inline SeparationPlane::SeparationPlane(Plane const& plane)
+      : plane_(plane) {}
+
+  inline bool SeparationPlane::contains(Point const& p) const {
+    return !plane_.isAbove(p);
+  }
+
+  inline std::string SeparationPlane::asString() const { return plane_.asString(); }
+
+} // namespace corsika
\ No newline at end of file

corsika/detail/framework/geometry/Sphere.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  inline bool Sphere::contains(Point const& p) const {
+    return radius_ * radius_ > (center_ - p).getSquaredNorm();
+  }
+
+  inline Point const& Sphere::getCenter() const { return center_; }
+
+  inline void Sphere::setCenter(Point const& p) { center_ = p; }
+
+  inline LengthType Sphere::getRadius() const { return radius_; }
+
+  inline void Sphere::setRadius(LengthType const r) { radius_ = r; }
+
+  inline CoordinateSystemPtr const Sphere::getCoordinateSystem() const {
+    return center_.getCoordinateSystem();
+  }
+
+  inline std::string Sphere::asString() const {
+    std::ostringstream txt;
+    txt << "center=" << center_ << ", radius=" << radius_;
+    return txt.str();
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/StraightTrajectory.inl

+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/Line.hpp>
+#include <corsika/framework/geometry/Point.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+namespace corsika {
+
+  inline VelocityVector StraightTrajectory::getVelocity(double const u) const {
+    return initialVelocity_ * (1 - u) + finalVelocity_ * u;
+  }
+
+  inline TimeType StraightTrajectory::getDuration(double const u) const {
+    return u * timeStep_;
+  }
+
+  template <typename Particle>
+  inline TimeType StraightTrajectory::getTime(Particle const& particle,
+                                              double const u) const {
+    return particle.getTime() + getDuration(u); // timeStep_ * u;
+  }
+
+  inline LengthType StraightTrajectory::getLength(double const u) const {
+    if (timeLength_ == 0_s) return 0_m;
+    if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s)
+      return std::numeric_limits<LengthType::value_type>::infinity() * 1_m;
+    return getDistance(u) * timeStep_ / timeLength_;
+  }
+
+  inline void StraightTrajectory::setLength(LengthType const limit) {
+    setDuration(line_.getTimeFromArclength(limit));
+  }
+
+  inline void StraightTrajectory::setDuration(TimeType const limit) {
+    if (timeStep_ == 0_s) {
+      timeLength_ = 0_s;
+      setFinalVelocity(getVelocity(0));
+      timeStep_ = limit;
+    } else {
+      // for infinite steps there can't be a difference between
+      // curved and straight trajectory, this is fundamentally
+      // undefined: assume they are the same (which, i.e. is always correct for a
+      // straight line trajectory).
+      //
+      // Final note: only straight-line trajectories should have
+      // infinite steps! Everything else is ill-defined.
+      if (timeStep_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s ||
+          timeLength_ == std::numeric_limits<TimeType::value_type>::infinity() * 1_s) {
+        timeLength_ = limit;
+        timeStep_ = limit;
+        // ...and don't touch velocity
+      } else {
+        const double scale = limit / timeStep_;
+        timeLength_ *= scale;
+        setFinalVelocity(getVelocity(scale));
+        timeStep_ = limit;
+      }
+    }
+  }
+
+  inline LengthType StraightTrajectory::getDistance(double const u) const {
+    assert(u <= 1);
+    assert(u >= 0);
+    return line_.getArcLength(0 * second, u * timeLength_);
+  }
+
+} // namespace corsika

corsika/detail/framework/geometry/Vector.inl

+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/geometry/BaseVector.hpp>
+#include <corsika/framework/geometry/QuantityVector.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+namespace corsika {
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> const& Vector<TDimension>::getComponents() const {
+    return BaseVector<TDimension>::getQuantityVector();
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension>& Vector<TDimension>::getComponents() {
+    return BaseVector<TDimension>::getQuantityVector();
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension> Vector<TDimension>::getComponents(
+      CoordinateSystemPtr const& pCS) const {
+    if (pCS == BaseVector<TDimension>::getCoordinateSystem()) {
+      return BaseVector<TDimension>::getQuantityVector();
+    } else {
+      return QuantityVector<TDimension>(
+          get_transformation(*BaseVector<TDimension>::getCoordinateSystem().get(),
+                             *pCS.get())
+              .linear() *
+          BaseVector<TDimension>::getQuantityVector().eigenVector_);
+    }
+  }
+
+  template <typename TDimension>
+  inline QuantityVector<TDimension>& Vector<TDimension>::getComponents(
+      CoordinateSystemPtr const& pCS) {
+    if (*pCS != *BaseVector<TDimension>::getCoordinateSystem()) { rebase(pCS); }
+    return BaseVector<TDimension>::getQuantityVector();
+  }
+
+  template <typename TDimension>
+  inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getX(
+      CoordinateSystemPtr const& pCS) const {
+    if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+      return BaseVector<TDimension>::getQuantityVector()[0];
+    } else {
+      return QuantityVector<TDimension>(
+          get_transformation(*BaseVector<TDimension>::getCoordinateSystem().get(),
+                             *pCS.get())
+              .linear() *
+          BaseVector<TDimension>::getQuantityVector().eigenVector_)[0];
+    }
+  }
+
+  template <typename TDimension>
+  inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getY(
+      CoordinateSystemPtr const& pCS) const {
+    if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+      return BaseVector<TDimension>::getQuantityVector()[1];
+    } else {
+      return QuantityVector<TDimension>(
+          get_transformation(*BaseVector<TDimension>::getCoordinateSystem().get(),
+                             *pCS.get())
+              .linear() *
+          BaseVector<TDimension>::getQuantityVector().eigenVector_)[1];
+    }
+  }
+
+  template <typename TDimension>
+  inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getZ(
+      CoordinateSystemPtr const& pCS) const {
+    if (*pCS == *BaseVector<TDimension>::getCoordinateSystem()) {
+      return BaseVector<TDimension>::getQuantityVector()[2];
+    } else {
+      return QuantityVector<TDimension>(
+          get_transformation(*BaseVector<TDimension>::getCoordinateSystem().get(),
+                             *pCS.get())
+              .linear() *
+          BaseVector<TDimension>::getQuantityVector().eigenVector_)[2];
+    }
+  }
+
+  template <typename TDimension>
+  inline void Vector<TDimension>::rebase(CoordinateSystemPtr const& pCS) {
+    BaseVector<TDimension>::setQuantityVector(QuantityVector<TDimension>(
+        get_transformation(*BaseVector<TDimension>::getCoordinateSystem().get(),
+                           *pCS.get())
+            .linear() *
+        BaseVector<TDimension>::getQuantityVector().eigenVector_));
+    BaseVector<TDimension>::setCoordinateSystem(pCS);
+  }
+
+  template <typename TDimension>
+  inline typename Vector<TDimension>::quantity_type Vector<TDimension>::getNorm() const {
+    return BaseVector<TDimension>::getQuantityVector().getNorm();
+  }
+
+  template <typename TDimension>
+  inline typename Vector<TDimension>::quantity_square_type
+  Vector<TDimension>::getSquaredNorm() const {
+    return BaseVector<TDimension>::getQuantityVector().getSquaredNorm();
+  }
+
+  template <typename TDimension>
+  template <typename TDimension2>
+  inline auto Vector<TDimension>::getParallelProjectionOnto(
+      Vector<TDimension2> const& pVec, CoordinateSystemPtr const& pCS) const {
+    auto const ourCompVec = getComponents(pCS);
+    auto const otherCompVec = pVec.getComponents(pCS);
+    auto const& a = ourCompVec.eigenVector_;
+    auto const& b = otherCompVec.eigenVector_;
+
+    return Vector<TDimension>(
+        pCS, QuantityVector<TDimension>(b * ((a.dot(b)) / b.squaredNorm())));
+  }
+
+  template <typename TDimension>
+  template <typename TDimension2>
+  inline auto Vector<TDimension>::getParallelProjectionOnto(
+      Vector<TDimension2> const& pVec) const {
+    return getParallelProjectionOnto<TDimension2>(
+        pVec, BaseVector<TDimension>::getCoordinateSystem());
+  }
+
+  template <typename TDimension>
+  inline Vector<TDimension> Vector<TDimension>::operator+(
+      Vector<TDimension> const& pVec) const {
+    CoordinateSystemPtr const& cs = BaseVector<TDimension>::getCoordinateSystem();
+    auto const components = getComponents(cs) + pVec.getComponents(cs);
+    return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(), components);
+  }
+
+  template <typename TDimension>
+  inline Vector<TDimension> Vector<TDimension>::operator-(
+      Vector<TDimension> const& pVec) const {
+    CoordinateSystemPtr const& cs = BaseVector<TDimension>::getCoordinateSystem();
+    return Vector<TDimension>(cs, getComponents() - pVec.getComponents(cs));
+  }
+
+  template <typename TDimension>
+  inline auto& Vector<TDimension>::operator*=(double const p) {
+    BaseVector<TDimension>::getQuantityVector() *= p;
+    return *this;
+  }
+
+  template <typename TDimension>
+  template <typename TScalarDim>
+  inline auto Vector<TDimension>::operator*(
+      phys::units::quantity<TScalarDim, double> const p) const {
+    using ProdDim = phys::units::detail::product_d<TDimension, TScalarDim>;
+
+    return Vector<ProdDim>(BaseVector<TDimension>::getCoordinateSystem(),
+                           BaseVector<TDimension>::getQuantityVector() * p);
+  }
+
+  template <typename TDimension>
+  template <typename TScalarDim>
+  inline auto Vector<TDimension>::operator/(
+      phys::units::quantity<TScalarDim, double> const p) const {
+    return (*this) * (1 / p);
+  }
+
+  template <typename TDimension>
+  inline auto Vector<TDimension>::operator*(double const p) const {
+    return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
+                              BaseVector<TDimension>::getQuantityVector() * p);
+  }
+
+  template <typename TDimension>
+  inline auto Vector<TDimension>::operator/(double const p) const {
+    return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
+                              BaseVector<TDimension>::getQuantityVector() / p);
+  }
+
+  template <typename TDimension>
+  inline auto& Vector<TDimension>::operator+=(Vector<TDimension> const& pVec) {
+    BaseVector<TDimension>::getQuantityVector() +=
+        pVec.getComponents(BaseVector<TDimension>::getCoordinateSystem());
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto& Vector<TDimension>::operator-=(Vector<TDimension> const& pVec) {
+    BaseVector<TDimension>::getQuantityVector() -=
+        pVec.getComponents(BaseVector<TDimension>::getCoordinateSystem());
+    return *this;
+  }
+
+  template <typename TDimension>
+  inline auto Vector<TDimension>::operator-() const {
+    return Vector<TDimension>(BaseVector<TDimension>::getCoordinateSystem(),
+                              -BaseVector<TDimension>::getQuantityVector());
+  }
+
+  template <typename TDimension>
+  inline auto Vector<TDimension>::normalized() const {
+    return (*this) * (1 / getNorm());
+  }
+
+  template <typename TDimension>
+  template <typename TDimension2>
+  inline auto Vector<TDimension>::cross(Vector<TDimension2> const& pV) const {
+    auto const c1 = getComponents().eigenVector_;
+    auto const c2 =
+        pV.getComponents(BaseVector<TDimension>::getCoordinateSystem()).eigenVector_;
+    auto const bareResult = c1.cross(c2);
+
+    using ProdDim = phys::units::detail::product_d<TDimension, TDimension2>;
+    return Vector<ProdDim>(BaseVector<TDimension>::getCoordinateSystem(), bareResult);
+  }
+
+  template <typename TDimension>
+  template <typename TDimension2>
+  inline auto Vector<TDimension>::dot(Vector<TDimension2> const& pV) const {
+    auto const c1 = getComponents().eigenVector_;
+    auto const c2 =
+        pV.getComponents(BaseVector<TDimension>::getCoordinateSystem()).eigenVector_;
+    auto const bareResult = c1.dot(c2);
+
+    using ProdDim = phys::units::detail::product_d<TDimension, TDimension2>;
+
+    return phys::units::quantity<ProdDim, double>(phys::units::detail::magnitude_tag,
+                                                  bareResult);
+  }
+
+  template <typename TDimension>
+  inline std::ostream& operator<<(std::ostream& os,
+                                  corsika::Vector<TDimension> const& v) {
+    auto const& qv = v.getComponents();
+    os << qv << " (ref:" << fmt::ptr(v.getCoordinateSystem()) << ")";
+    return os;
+  }
+
+  /*
+   * scalar * vector multiplication
+   */
+  template <typename TDimension, typename UDimension>
+  inline Vector<phys::units::detail::product_d<TDimension, UDimension>> operator*(
+      quantity<UDimension> const n, Vector<TDimension> const& vec) {
+    return vec * n;
+  }
+
+  template <typename TDimension>
+  inline Vector<TDimension> operator*(double const n, Vector<TDimension> const& vec) {
+    return vec * n;
+  }
+} // namespace corsika

corsika/detail/framework/process/BoundaryCrossingProcess.hpp

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+namespace corsika {
+
+  /**
+     traits test for BoundaryCrossingProcess::doBoundaryCrossing method
+  */
+  template <class TProcess, typename TReturn, typename TParticle>
+  struct has_method_doBoundaryCrossing
+      : public detail::has_method_signature<TReturn, TParticle&,
+                                            typename TParticle::node_type const&,
+                                            typename TParticle::node_type const&> {
+
+    ///! method signature
+    using detail::has_method_signature<
+        TReturn, TParticle&, typename TParticle::node_type const&,
+        typename TParticle::node_type const&>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template doBoundaryCrossing<TParticle>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::doBoundaryCrossing)) test(std::nullptr_t);
+
+  public:
+    /**
+        @name traits results
+        @{
+    */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  //! @file BoundaryCrossingProcess.hpp
+  //! value traits type
+
+  template <class TProcess, typename TReturn, typename TParticle>
+  bool constexpr has_method_doBoundaryCrossing_v =
+      has_method_doBoundaryCrossing<TProcess, TReturn, TParticle>::value;
+
+} // namespace corsika

corsika/detail/framework/process/CascadeEquationsProcess.hpp

+/*
+ * (c) Copyright 2021 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the 3-clause BSD license.
+ * See file LICENSE for a full version of the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/process/ProcessTraits.hpp>
+#include <corsika/framework/utility/HasMethodSignature.hpp>
+
+/**
+ * @file CascadeEquationsProcess.hpp
+ */
+
+namespace corsika {
+
+  /**
+   * traits test for CascadeEquationsProcess::doCascadeEquations method.
+   */
+  template <class TProcess, typename TReturn, typename... TArg>
+  struct has_method_doCascadeEquations
+      : public detail::has_method_signature<TReturn, TArg...> {
+
+    //! method signature
+    using detail::has_method_signature<TReturn, TArg...>::testSignature;
+
+    //! the default value
+    template <class T>
+    static std::false_type test(...);
+
+    //! templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::template doCascadeEquations<TArg...>)) test(
+        std::nullptr_t);
+
+    //! non templated parameter option
+    template <class T>
+    static decltype(testSignature(&T::doCascadeEquations)) test(std::nullptr_t);
+
+  public:
+    /**
+     *  @name traits results
+     *  @{
+     */
+    using type = decltype(test<std::decay_t<TProcess>>(nullptr));
+    static const bool value = type::value;
+    //! @}
+  };
+
+  /**
+   * value traits type.
+   */
+  template <class TProcess, typename TReturn, typename... TArg>
+  bool constexpr has_method_doCascadeEquations_v =
+      has_method_doCascadeEquations<TProcess, TReturn, TArg...>::value;
+
+} // namespace corsika