propagate only in forward direction, note that testCascade is failing currently

Documentation/Examples/cascade_example.cc

@@ -245,11 +245,10 @@ double s_rndm_(int&) {
 }
 
 int main() {
-  Environment env;
-  //~ auto& universe = env.GetUniverse();
+  corsika::environment::Environment env; // dummy environment  
   auto& universe = *(env.GetUniverse());
 
-  auto const theMedium = Environment::CreateNode<Sphere>(
+  auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
       Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
       
   universe.AddChild(std::move(theMedium));

Environment/Environment.h

@@ -17,13 +17,19 @@
 #include <corsika/setup/SetupEnvironment.h>
 
 namespace corsika::environment {
+  struct Universe : public corsika::geometry::Volume {
+      bool Contains(corsika::geometry::Point const&) const override {return true;}
+  };
 
   class Environment {
   public:
+    Environment() : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+          std::make_unique<Universe>())) {}
+  
     using IEnvironmentModel = corsika::setup::IEnvironmentModel;
 
-    auto& GetUniverse() { return universe; }
-    auto const& GetUniverse() const { return universe; }
+    auto& GetUniverse() { return fUniverse; }
+    auto const& GetUniverse() const { return fUniverse; }
 
     auto const& GetCoordinateSystem() const {
       return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
@@ -41,7 +47,7 @@ namespace corsika::environment {
     }
 
   private:
-    VolumeTreeNode<IEnvironmentModel>::VTNUPtr universe;
+    VolumeTreeNode<IEnvironmentModel>::VTNUPtr fUniverse;
   };
 
 } // namespace corsika::environment

Framework/Cascade/testCascade.cc

@@ -78,8 +78,13 @@ private:
 
 TEST_CASE("Cascade", "[Cascade]") {
   corsika::environment::Environment env; // dummy environment
+  auto& universe = *(env.GetUniverse());
+  auto const radius = 20_km;
+  auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
+    Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);      
+  universe.AddChild(std::move(theMedium));
 
-  tracking_line::TrackingLine<setup::Stack> tracking(env);
+  tracking_line::TrackingLine<setup::Stack> tracking(env);  
 
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;

Framework/Geometry/Vector.h

@@ -190,7 +190,7 @@ namespace corsika::geometry {
 
       using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
 
-      return ProdQuantity(bareResult);
+      return ProdQuantity(phys::units::detail::magnitude_tag, bareResult);
     }
   };
 

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/TrackingLine.h

@@ -14,6 +14,7 @@
 #include <algorithm>
 #include <iostream>
 #include <optional>
+#include <utility>
 
 using namespace corsika;
 
@@ -28,7 +29,7 @@ namespace corsika::process {
       corsika::environment::Environment const& fEnvironment;
 
     public:
-      std::optional<corsika::units::si::TimeType> TimeOfIntersection(
+      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();
@@ -49,8 +50,8 @@ namespace corsika::process {
         //~ std::cout << "beta: " << beta << std::endl;
 
         if (discriminant.magnitude() > 0) {
-          auto const t = (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
-          return t;
+          (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
+          return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()), (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
         } else {
           return {};
         }
@@ -82,7 +83,13 @@ namespace corsika::process {
                        // everything is a sphere, crashes with exception if not
 
           if (auto opt = TimeOfIntersection(line, sphere); opt.has_value())
-            intersectionTimes.push_back(*opt);
+          {
+            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); }

Processes/TrackingLine/testTrackingLine.cc

@@ -32,12 +32,27 @@ using namespace corsika::geometry;
 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<setup::Stack> tracking(env);
+  tracking_line::TrackingLine<DummyStack> tracking(env);
 
   SECTION("intersection with sphere") {
     Point const origin(cs, {0_m, 0_m, 0_m});
@@ -50,7 +65,10 @@ TEST_CASE("TrackingLine") {
     
     auto const opt = tracking.TimeOfIntersection(traj, Sphere(Point(cs, {0_m, 0_m, 10_m}), 1_m));
     REQUIRE(opt.has_value());
-    REQUIRE(opt.value() / 9_s == Approx(1));
+    
+    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());
@@ -59,9 +77,22 @@ TEST_CASE("TrackingLine") {
   SECTION("maximally possible propagation") {
       auto& universe = *(env.GetUniverse());
       
-      auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
-      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
+      //~ 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));
   }
 }