Framework/Geometry/Volume.h

+#ifndef _include_VOLUME_H_
+#define _include_VOLUME_H_
+
+#include <corsika/geometry/Point.h>
+
+namespace corsika::geometry {
+
+  class Volume {
+
+  public:
+    //! returns true if the Point p is within the volume
+    virtual bool Contains(Point const& p) const = 0;
+
+    virtual ~Volume() = default;
+  };
+
+} // namespace corsika::geometry
+
+#endif

Framework/Geometry/testGeometry.cc

@@ -132,7 +132,15 @@ TEST_CASE("Sphere") {
   Point center(rootCS, {0_m, 3_m, 4_m});
   Sphere sphere(center, 5_m);
 
-  SECTION("isInside") {
+  SECTION("GetCenter") {
+    CHECK((sphere.GetCenter().GetCoordinates(rootCS) -
+           QuantityVector<length_d>(0_m, 3_m, 4_m))
+              .norm()
+              .magnitude() == Approx(0).margin(absMargin));
+    CHECK(sphere.GetRadius() / 5_m == Approx(1));
+  }
+
+  SECTION("Contains") {
     REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
     REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
   }
@@ -152,11 +160,11 @@ TEST_CASE("Trajectories") {
             .norm()
             .magnitude() == Approx(0).margin(absMargin));
 
-    Trajectory<Line> base(line, 1_s);
-    CHECK(line.GetPosition(2_s).GetCoordinates() ==
-          base.GetPosition(2_s).GetCoordinates());
+    auto const t = 1_s;
+    Trajectory<Line> base(line, t);
+    CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
-    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(1));
+    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(1));
   }
 
   SECTION("Helix") {
@@ -180,10 +188,10 @@ TEST_CASE("Trajectories") {
               .norm()
               .magnitude() == Approx(0).margin(absMargin));
 
-    Trajectory<Helix> const base(helix, 1_s);
-    CHECK(helix.GetPosition(1234_s).GetCoordinates() ==
-          base.GetPosition(1234_s).GetCoordinates());
+    auto const t = 1234_s;
+    Trajectory<Helix> const base(helix, t);
+    CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
 
-    CHECK(base.GetDistanceBetween(1_s, 2_s) / 1_m == Approx(5));
+    CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(5));
   }
 }

Framework/Units/PhysicalUnits.h

@@ -51,6 +51,8 @@ namespace corsika::units::si {
       phys::units::quantity<phys::units::electric_charge_d, double>;
   using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
   using MassType = phys::units::quantity<phys::units::mass_d, double>;
+  using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
+  using GrammageType = phys::units::quantity<phys::units::dimensions<-2, 1, 0>, double>;
   using MomentumType = phys::units::quantity<momentum_d, double>;
   using CrossSectionType = phys::units::quantity<sigma_d, double>;
 

Processes/StackInspector/StackInspector.cc

@@ -12,6 +12,7 @@
 #include <corsika/geometry/RootCoordinateSystem.h>
 #include <corsika/process/stack_inspector/StackInspector.h>
 #include <corsika/units/PhysicalUnits.h>
+#include <corsika/particles/ParticleProperties.h>
 
 #include <corsika/logging/Logger.h>
 

Processes/TrackingLine/CMakeLists.txt

@@ -22,17 +22,16 @@ target_include_directories (
 
 install (FILES ${MODEL_HEADERS} DESTINATION include/${MODEL_NAMESPACE})
 
-
-
-# # --------------------
-# # code unit testing
-# add_executable (testStackInspector testStackInspector.cc)
-
-# target_link_libraries (
-#   testStackInspector
-#   CORSIKAunits
-#   CORSIKAthirdparty # for catch2
-#   )
-
-# add_test (NAME testStackInspector COMMAND testStackInspector)
-
+# #-- -- -- -- -- -- -- -- -- --
+# #code unit testing
+add_executable (testTrackingLine testTrackingLine.cc)
+
+target_link_libraries (
+   testTrackingLine
+   CORSIKAunits
+   CORSIKAenvironment
+   CORSIKAgeometry
+   CORSIKAthirdparty # for catch2
+)
+
+add_test (NAME testTrackingLine COMMAND testTrackingLine)

Processes/TrackingLine/TrackingLine.h

-#ifndef _include_corsika_processes_TrackinLine_h_
-#define _include_corsika_processes_TrackinLine_h_
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _include_corsika_processes_TrackingLine_h_
+#define _include_corsika_processes_TrackingLine_h_
 
 #include <corsika/geometry/Point.h>
+#include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Vector.h>
 
 #include <corsika/units/PhysicalUnits.h>
 
+#include <corsika/environment/Environment.h>
 #include <corsika/setup/SetupStack.h>
 #include <corsika/setup/SetupTrajectory.h>
 
+#include <algorithm>
+#include <iostream>
+#include <optional>
+#include <stdexcept>
+#include <utility>
+
 using namespace corsika;
 
 namespace corsika::process {
@@ -16,15 +34,96 @@ namespace corsika::process {
   namespace tracking_line {
 
     template <typename Stack>
-    class TrackingLine { // Newton-step, naja.. not yet
+    class TrackingLine { //
       typedef typename Stack::ParticleType Particle;
 
+      corsika::environment::Environment const& fEnvironment;
+
     public:
+      std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+      TimeOfIntersection(corsika::geometry::Line const& traj,
+                         geometry::Sphere const& sphere) {
+        using namespace corsika::units::si;
+        auto const& cs = fEnvironment.GetCoordinateSystem();
+        geometry::Point const origin(cs, 0_m, 0_m, 0_m);
+
+        auto const r0 = (traj.GetR0() - origin);
+        auto const v0 = traj.GetV0();
+        auto const c0 = (sphere.GetCenter() - origin);
+
+        auto const alpha = r0.dot(v0) - 2 * v0.dot(c0);
+        auto const beta = c0.squaredNorm() + r0.squaredNorm() + 2 * c0.dot(r0) -
+                          sphere.GetRadius() * sphere.GetRadius();
+
+        auto const discriminant = alpha * alpha - 4 * beta * v0.squaredNorm();
+
+        //~ std::cout << "discriminant: " << discriminant << std::endl;
+        //~ std::cout << "alpha: " << alpha << std::endl;
+        //~ std::cout << "beta: " << beta << std::endl;
+
+        if (discriminant.magnitude() > 0) {
+          (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
+          return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
+                                (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
+        } else {
+          return {};
+        }
+      }
+
+      TrackingLine(corsika::environment::Environment const& pEnv)
+          : fEnvironment(pEnv) {}
       void Init() {}
-      setup::Trajectory GetTrack(Particle& p) {
-        geometry::Vector<SpeedType::dimension_type> v = p.GetDirection();
-        geometry::Line traj(p.GetPosition(), v);
-        return geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
+
+      auto GetTrack(Particle& p) {
+        using namespace corsika::units::si;
+        geometry::Vector<SpeedType::dimension_type> const velocity =
+            p.GetMomentum() / p.GetEnergy() * corsika::units::si::constants::cSquared;
+
+        auto const currentPosition = p.GetPosition();
+        geometry::Line line(currentPosition, velocity);
+
+        auto const* currentVolumeNode =
+            fEnvironment.GetUniverse()->GetContainingNode(currentPosition);
+        auto const& children = currentVolumeNode->GetChildNodes();
+        auto const& excluded = currentVolumeNode->GetExcludedNodes();
+
+        std::vector<TimeType> intersectionTimes;
+
+        auto addIfIntersects = [&](auto& vtn) {
+          auto const& volume = vtn.GetVolume();
+          auto const& sphere = dynamic_cast<geometry::Sphere const&>(
+              volume); // for the moment we are a bit bold here and assume
+                       // everything is a sphere, crashes with exception if not
+
+          if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
+            auto const [t1, t2] = *opt;
+            if (t1.magnitude() >= 0)
+              intersectionTimes.push_back(t1);
+            else if (t2.magnitude() >= 0)
+              intersectionTimes.push_back(t2);
+          }
+        };
+
+        for (auto const& child : children) { addIfIntersects(*child); }
+
+        for (auto const* child : excluded) { addIfIntersects(*child); }
+
+        addIfIntersects(*currentVolumeNode);
+
+        auto const minIter =
+            std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
+
+        TimeType min;
+
+        if (minIter == intersectionTimes.cend()) {
+          min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+                     // to handle the numerics properly
+          //~ throw std::runtime_error("no intersection with anything!");
+        } else {
+          min = *minIter;
+        }
+
+        return geometry::Trajectory<corsika::geometry::Line>(line, min);
       }
     };
 

Processes/TrackingLine/testTrackingLine.cc

+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#include <corsika/environment/Environment.h>
+
+#include <corsika/process/tracking_line/TrackingLine.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/geometry/Sphere.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+using corsika::setup::Trajectory;
+
+#define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
+                          // cpp file
+#include <catch2/catch.hpp>
+
+using namespace corsika;
+using namespace corsika::process;
+using namespace corsika::units;
+using namespace corsika::geometry;
+
+#include <iostream>
+using namespace std;
+using namespace corsika::units::si;
+
+struct DummyParticle {
+    EnergyType fEnergy;
+    Vector<momentum_d> fMomentum;
+    Point fPosition;
+    
+    DummyParticle(EnergyType pEnergy, Vector<momentum_d> pMomentum, Point pPosition) : fEnergy(pEnergy), fMomentum(pMomentum), fPosition(pPosition) {}
+    
+    auto GetEnergy() const { return fEnergy; }
+    auto GetMomentum() const { return fMomentum; }
+    auto GetPosition() const { return fPosition; }
+};
+
+struct DummyStack {
+    using ParticleType = DummyParticle;
+};
+
+TEST_CASE("TrackingLine") {
+  corsika::environment::Environment env; // dummy environment
+  auto const& cs = env.GetCoordinateSystem();
+  
+  tracking_line::TrackingLine<DummyStack> tracking(env);
+
+  SECTION("intersection with sphere") {
+    Point const origin(cs, {0_m, 0_m, 0_m});
+    Point const center(cs, {0_m, 0_m, 10_m});
+    Sphere const sphere(center, 1_m);
+    Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+    Line line(origin, v);
+    
+    geometry::Trajectory<Line> traj(line, 12345_s);
+    
+    auto const opt = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
+    REQUIRE(opt.has_value());
+    
+    auto [t1, t2] = opt.value();
+    REQUIRE(t1 / 9_s == Approx(1));
+    REQUIRE(t2 / 11_s == Approx(1));
+    
+    auto const optNoIntersection = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {5_m, 0_m, 10_m}), 1_m));
+    REQUIRE_FALSE(optNoIntersection.has_value());
+  }
+  
+  SECTION("maximally possible propagation") {
+      auto& universe = *(env.GetUniverse());
+      
+      //~ std::cout << env.GetUniverse().get() << std::endl;
+      
+      DummyParticle p(1_J, Vector<momentum_d>(cs, 0*kilogram*meter/second, 0*kilogram*meter/second, 1*kilogram*meter/second), Point(cs, 0_m, 0_m,0_m));
+      
+      auto const radius = 20_m;
+      
+      auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);      
+      universe.AddChild(std::move(theMedium));
+      
+      Point const origin(cs, {0_m, 0_m, 0_m});
+      Vector<corsika::units::si::SpeedType::dimension_type> v(cs, 0_m/second, 0_m/second, 1_m/second);
+      Line line(origin, v);
+      
+      auto const traj = tracking.GetTrack(p);
+      
+      REQUIRE((traj.GetPosition(1.) - Point(cs, 0_m, 0_m, radius)).GetComponents(cs).norm().magnitude() == Approx(0).margin(1e-4));
+  }
+}

Setup/SetupEnvironment.h

@@ -12,6 +12,10 @@
 #ifndef _include_corsika_setup_environment_h_
 #define _include_corsika_setup_environment_h_
 
-namespace corsika {}
+#include <corsika/environment/IMediumModel.h>
+
+namespace corsika::setup {
+    using IEnvironmentModel = corsika::environment::IMediumModel;
+}
 
 #endif