added Plane class

Framework/Geometry/CMakeLists.txt

@@ -10,6 +10,7 @@ set (
   Point.h
   Line.h
   Sphere.h
+  Plane.h
   Volume.h
   CoordinateSystem.h
   RootCoordinateSystem.h

Framework/Geometry/Plane.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_Framework_Geometry_Plane_h_
+#define _include_Framework_Geometry_Plane_h_
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/units/PhysicalUnits.h>
+
+namespace corsika::geometry {
+  class Plane {
+
+    using DimLessVec = Vector<corsika::units::si::dimensionless_d>;
+
+    Point const fCenter;
+    DimLessVec const fNormal;
+
+  public:
+    Plane(Point const& vCenter, DimLessVec const& vNormal)
+        : fCenter(vCenter)
+        , fNormal(vNormal) {}
+    bool IsAbove(Point const& vP) const {
+      return fNormal.dot(vP - fCenter) > corsika::units::si::LengthType::zero();
+    }
+
+    Point const& GetCenter() const { return fCenter; }
+    DimLessVec const& GetNormal() const { return fNormal; }
+  };
+
+} // namespace corsika::geometry
+
+#endif

Processes/TrackingLine/TrackingLine.cc

@@ -16,31 +16,27 @@
 #include <corsika/geometry/Vector.h>
 #include <corsika/process/tracking_line/TrackingLine.h>
 
-#include <algorithm>
-#include <iostream>
+#include <limits>
 #include <stdexcept>
 #include <utility>
 
-using namespace corsika;
+using namespace corsika::geometry;
+using namespace corsika::units::si;
 
 namespace corsika::process::tracking_line {
 
-  std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
-  TimeOfIntersection(corsika::geometry::Line const& line,
-                     geometry::Sphere const& sphere) {
+  std::optional<std::pair<TimeType, TimeType>> TimeOfIntersection(Line const& line,
+                                                                  Sphere const& sphere) {
     auto const delta = line.GetR0() - sphere.GetCenter();
     auto const v = line.GetV0();
-    auto const vSqNorm = v.squaredNorm();
+    auto const vSqNorm =
+        v.squaredNorm(); // todo: get rid of this by having V0 normalized always
     auto const R = sphere.GetRadius();
 
     auto const vDotDelta = v.dot(delta);
     auto const discriminant =
         vDotDelta * vDotDelta - vSqNorm * (delta.squaredNorm() - R * R);
 
-    //~ std::cout << "discriminant: " << discriminant << std::endl;
-    //~ std::cout << "alpha: " << alpha << std::endl;
-    //~ std::cout << "beta: " << beta << std::endl;
-
     if (discriminant.magnitude() > 0) {
       auto const sqDisc = sqrt(discriminant);
       auto const invDenom = 1 / vSqNorm;
@@ -50,4 +46,17 @@ namespace corsika::process::tracking_line {
       return {};
     }
   }
+
+  TimeType TimeOfIntersection(Line const& vLine, Plane const& vPlane) {
+    auto const delta = vPlane.GetCenter() - vLine.GetR0();
+    auto const v = vLine.GetV0();
+    auto const n = vPlane.GetNormal();
+    auto const c = n.dot(v);
+
+    if (c.magnitude() == 0) {
+      return std::numeric_limits<TimeType::value_type>::infinity() * 1_s;
+    } else {
+      return n.dot(delta) / c;
+    }
+  }
 } // namespace corsika::process::tracking_line

Processes/TrackingLine/TrackingLine.h

@@ -13,6 +13,7 @@
 #define _include_corsika_processes_TrackingLine_h_
 
 #include <corsika/geometry/Line.h>
+#include <corsika/geometry/Plane.h>
 #include <corsika/geometry/Sphere.h>
 #include <corsika/geometry/Trajectory.h>
 #include <corsika/geometry/Vector.h>
@@ -33,8 +34,10 @@ namespace corsika::process {
   namespace tracking_line {
 
     std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
-    TimeOfIntersection(corsika::geometry::Line const& line,
-                       geometry::Sphere const& sphere);
+    TimeOfIntersection(geometry::Line const&, geometry::Sphere const&);
+
+    corsika::units::si::TimeType TimeOfIntersection(geometry::Line const& line,
+                                                    geometry::Plane const&);
 
     class TrackingLine {