IAP GITLAB

Skip to content
Snippets Groups Projects
Intersect.hpp 5.82 KiB
/*
 * (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/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <corsika/logging/Logging.h>
#include <corsika/geometry/Intersections.hpp>

#include <limits>

namespace corsika::process::tracking {

  /**
   * \class Intersect
   *
   * This is a CRTP class to provide a generic volume-tree
   * intersection for the purpose of tracking.
   *
   * It return the closest distance in time to the next geometric
   * intersection, as well as a pointer to the corresponding new
   * volume.
   *
   * User may provide an optional global step-length limit as
   * parameter. This may be needd for (simpler) algorithms in magnetic
   * fields, where tracking errors grow linearly with step-length.
   * Obviously, in the case of the step-length limit, the returend
   * "next" volume is just the current one.
   *
   **/

  template <typename TDerived>
  class Intersect {

  protected:
    template <typename TParticle>
    auto nextIntersect(
        const TParticle& particle,
        corsika::units::si::TimeType step_limit =
            std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() *
            corsika::units::si::second) const {
      using namespace corsika::units::si;
      using namespace corsika::geometry;

      typedef
          typename std::remove_reference<decltype(*particle.GetNode())>::type node_type;
      node_type& volumeNode =
          *particle.GetNode(); // current "logical" node, from previous tracking step
      C8LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode),
                  volumeNode.GetVolume().Contains(particle.GetPosition()));

      // start values:
      TimeType minTime = step_limit;
      node_type* minNode = &volumeNode;

      // determine the first geometric collision with any other Volume boundary

      // first check, where we leave the current volume
      // this assumes our convention, that all Volume primitives must be convex
      // thus, the last entry is always the exit point
      const Intersections time_intersections_curr =
          TDerived::Intersect(particle, volumeNode);
      C8LOG_TRACE("curr node {}, parent node {}, hasIntersections={} ",
                  fmt::ptr(&volumeNode), fmt::ptr(volumeNode.GetParent()),
                  time_intersections_curr.hasIntersections());
      if (time_intersections_curr.hasIntersections()) {
        C8LOG_DEBUG("intersection times with currentLogicalVolumeNode: {} s and {} s",
                    time_intersections_curr.getEntry() / 1_s,
                    time_intersections_curr.getExit() / 1_s);
        if (time_intersections_curr.getExit() <= minTime) {
          minTime =
              time_intersections_curr.getExit(); // we exit currentLogicalVolumeNode here
          minNode = volumeNode.GetParent();
        }
      }

      // where do we collide with any of the next-tree-level volumes
      // entirely contained by currentLogicalVolumeNode
      for (const auto& node : volumeNode.GetChildNodes()) {

        const Intersections time_intersections = TDerived::Intersect(particle, *node);
        if (!time_intersections.hasIntersections()) { continue; }
        C8LOG_DEBUG("intersection times with child volume {} : enter {} s, exit {} s",
                    fmt::ptr(node), time_intersections.getEntry() / 1_s,
                    time_intersections.getExit() / 1_s);

        const auto t_entry = time_intersections.getEntry();
        const auto t_exit = time_intersections.getExit();
        C8LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
                    t_entry <= minTime);
        // note, theoretically t can even be smaller than 0 since we
        // KNOW we can't yet be in this volume yet, so we HAVE TO
        // enter it IF exit point is not also in the "past", AND if
        // extry point is [[much much]] closer than exit point
        // (because we might have already numerically "exited" it)!
        if (t_exit > 0_s && t_entry <= minTime &&
            -t_entry < t_exit) { // protection agains numerical problem, when we already
                                 // _exited_ before
                                 // enter chile volume here
          minTime = t_entry;
          minNode = node.get();
        }
      }

      // these are volumes from the previous tree-level that are cut-out partly from the
      // current volume
      for (node_type* node : volumeNode.GetExcludedNodes()) {

        const Intersections time_intersections = TDerived::Intersect(particle, *node);
        if (!time_intersections.hasIntersections()) { continue; }
        C8LOG_DEBUG("intersection times with exclusion volume {} : enter {} s, exit {} s",
                    fmt::ptr(node), time_intersections.getEntry() / 1_s,
                    time_intersections.getExit() / 1_s);
        const auto t_entry = time_intersections.getEntry();
        const auto t_exit = time_intersections.getExit();
        C8LOG_TRACE("children t-entry: {}, t-exit: {}, smaller? {} ", t_entry, t_exit,
                    t_entry <= minTime);
        // note, theoretically t can even be smaller than 0 since we
        // KNOW we can't yet be in this volume yet, so we HAVE TO
        // enter it IF exit point is not also in the "past"!
        if (t_exit > 0_s && t_entry <= minTime) { // enter volumen child here
          minTime = t_entry;
          minNode = node;
        }
      }
      C8LOG_TRACE("t-intersect: {}, node {} ", minTime, fmt::ptr(minNode));
      return std::make_tuple(minTime, minNode);
    }
  }; // namespace corsika::process::tracking
} // namespace corsika::process::tracking