From 145f56b49ee6b35802bd91a06819361c83f21f9f Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Tue, 4 Dec 2018 20:42:01 +0100
Subject: [PATCH] intersection of line and sphere
---
Documentation/Examples/cascade_example.cc | 20 +++++---
Environment/Environment.h | 19 ++++++-
Environment/IMediumModel.h | 1 +
Environment/VolumeTreeNode.h | 11 ----
Framework/Cascade/Cascade.h | 2 -
Framework/Cascade/testCascade.cc | 11 ++--
Framework/Geometry/Line.h | 5 +-
Framework/Geometry/Sphere.h | 3 ++
Framework/Geometry/Vector.h | 11 ++++
Framework/Geometry/testGeometry.cc | 16 ++++--
Processes/TrackingLine/CMakeLists.txt | 27 +++++-----
Processes/TrackingLine/TrackingLine.h | 53 +++++++++++++++++---
Processes/TrackingLine/testTrackingLine.cc | 58 ++++++++++++++++++++++
13 files changed, 184 insertions(+), 53 deletions(-)
create mode 100644 Processes/TrackingLine/testTrackingLine.cc
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 5caff6b79..9b664046b 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -23,6 +23,7 @@
#include <corsika/cascade/SibStack.h>
#include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/geometry/Sphere.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/units/PhysicalUnits.h>
@@ -39,6 +40,7 @@ using namespace corsika::geometry;
using namespace corsika::environment;
using namespace std;
+using namespace corsika::units::si;
static int fCount = 0;
@@ -47,7 +49,7 @@ public:
ProcessSplit() {}
template <typename Particle>
- double MinStepLength(Particle& p, setup::Trajectory&) const {
+ double MinStepLength(Particle& p, corsika::setup::Trajectory&) const {
// beam particles for sibyll : 1, 2, 3 for p, pi, k
// read from cross section code table
@@ -99,7 +101,7 @@ public:
}
template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, Trajectory&, Stack&) const {
+ EProcessReturn DoContinuous(Particle&, corsika::setup::Trajectory&, Stack&) const {
// corsika::utls::ignore(p);
return EProcessReturn::eOk;
}
@@ -242,12 +244,14 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
-int main() {
+int main() {
Environment env;
-
- auto theMedium = env.GetUniverse()::CreateNode<Sphere>({Point{env.GetCS(), 0_m, 0_m, 0_m}, 100_km});
-
- tracking_line::TrackingLine<setup::Stack> tracking(environment);
+ auto& universe = env.GetUniverse();
+
+ auto theMedium = Environment::CreateNode<Sphere>(
+ Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
+
+ tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
const auto sequence = p0 + p1;
@@ -257,7 +261,7 @@ int main() {
stack.Clear();
auto particle = stack.NewParticle();
- EnergyType E0 = 100_GeV;
+ EnergyType E0 = 1_EeV;
particle.SetEnergy(E0);
particle.SetPID(Code::Proton);
EAS.Init();
diff --git a/Environment/Environment.h b/Environment/Environment.h
index 3c20ba22a..0bc2fa112 100644
--- a/Environment/Environment.h
+++ b/Environment/Environment.h
@@ -20,11 +20,26 @@ namespace corsika::environment {
class Environment {
public:
+ using IEnvironmentModel = corsika::setup::IEnvironmentModel;
+
auto& GetUniverse() { return universe; }
- auto const& GetCS() const { return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();}
+ auto const& GetCoordinateSystem() const {
+ return corsika::geometry::RootCoordinateSystem::GetInstance().GetRootCS();
+ }
+
+ // factory method for creation of VolumeTreeNodes
+ template <class VolumeType, typename... Args>
+ static auto CreateNode(Args&&... args) {
+ static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
+ "unusable type provided, needs to be derived from "
+ "\"corsika::geometry::Volume\"");
+
+ return std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+ std::make_unique<VolumeType>(std::forward<Args>(args)...));
+ }
private:
- VolumeTreeNode<corsika::setup::IEnvironmentModel>::VTNUPtr universe;
+ VolumeTreeNode<IEnvironmentModel>::VTNUPtr universe;
};
} // namespace corsika::environment
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
index a80950831..4a1eedd86 100644
--- a/Environment/IMediumModel.h
+++ b/Environment/IMediumModel.h
@@ -2,6 +2,7 @@
#define _include_IMediumModel_h
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/Line.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Trajectory.h>
#include <corsika/units/PhysicalUnits.h>
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 523aebcd2..470bb8aa0 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -102,17 +102,6 @@ namespace corsika::environment {
return std::make_shared<MediumType>(std::forward<Args>(args)...);
}
- // factory method for creation of nodes
- template <class VolumeType, typename... Args>
- static auto CreateNode(Args&&... args) {
- static_assert(std::is_base_of_v<corsika::geometry::Volume, VolumeType>,
- "unusable type provided, needs to be derived from "
- "\"corsika::geometry::Volume\"");
-
- return std::make_unique<VolumeTreeNode<IModelProperties>>(
- std::make_unique<VolumeType>(std::forward<Args>(args)...));
- }
-
private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index e46960dc9..8207a036f 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -19,8 +19,6 @@
#include <corsika/setup/SetupTrajectory.h>
-using namespace corsika::units::si;
-
namespace corsika::cascade {
template <typename Tracking, typename ProcessList, typename Stack>
diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc
index ad4daa2e4..5017d5acf 100644
--- a/Framework/Cascade/testCascade.cc
+++ b/Framework/Cascade/testCascade.cc
@@ -9,6 +9,8 @@
* the license.
*/
+#include <corsika/environment/Environment.h>
+
#include <corsika/cascade/Cascade.h>
#include <corsika/process/ProcessSequence.h>
@@ -35,6 +37,7 @@ using namespace corsika::geometry;
#include <iostream>
using namespace std;
+using namespace corsika::units::si;
static int fCount = 0;
@@ -73,10 +76,10 @@ public:
private:
};
-
TEST_CASE("Cascade", "[Cascade]") {
+ corsika::environment::Environment env; // dummy environment
- tracking_line::TrackingLine<setup::Stack> tracking;
+ tracking_line::TrackingLine<setup::Stack> tracking(env);
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
@@ -92,8 +95,8 @@ TEST_CASE("Cascade", "[Cascade]") {
EnergyType E0 = 100_GeV;
particle.SetEnergy(E0);
particle.SetPosition(Point(rootCS, {0_m, 0_m, 10_km}));
- particle.SetMomentum(
- corsika::stack::super_stupid::MomentumVector(rootCS, {0*newton*second, 0*newton*second, -1*newton*second}));
+ particle.SetMomentum(corsika::stack::super_stupid::MomentumVector(
+ rootCS, {0 * newton * second, 0 * newton * second, -1 * newton * second}));
EAS.Init();
EAS.Run();
diff --git a/Framework/Geometry/Line.h b/Framework/Geometry/Line.h
index 80f6c16bf..a79ba9ce8 100644
--- a/Framework/Geometry/Line.h
+++ b/Framework/Geometry/Line.h
@@ -33,7 +33,7 @@ namespace corsika::geometry {
Point GetPosition(corsika::units::si::TimeType t) const { return r0 + v0 * t; }
LengthType ArcLength(corsika::units::si::TimeType t1,
- corsika::units::si::TimeType t2) const {
+ corsika::units::si::TimeType t2) const {
return v0.norm() * (t2 - t1);
}
@@ -41,6 +41,9 @@ namespace corsika::geometry {
corsika::units::si::LengthType t) const {
return t / v0.norm();
}
+
+ auto GetR0() const { return r0; }
+ auto GetV0() const { return v0; }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 3b6c52860..3e5458b7e 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -31,6 +31,9 @@ namespace corsika::geometry {
bool Contains(Point const& p) const override {
return fRadius * fRadius > (fCenter - p).squaredNorm();
}
+
+ auto& GetCenter() const { return fCenter; }
+ auto GetRadius() const { return fRadius; }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index cb4e49005..904b70398 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -181,6 +181,17 @@ namespace corsika::geometry {
bareResult);
}
}
+
+ template <typename dim2>
+ auto dot(Vector<dim2> pV) const {
+ auto const c1 = GetComponents().eVector;
+ auto const c2 = pV.GetComponents(*BaseVector<dim>::cs).eVector;
+ auto const bareResult = c1.dot(c2);
+
+ using ProdQuantity = phys::units::detail::Product<dim, dim2, double, double>;
+
+ return ProdQuantity(bareResult);
+ }
};
} // namespace corsika::geometry
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 63f119fc0..8130114aa 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/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})));
}
@@ -154,8 +162,7 @@ TEST_CASE("Trajectories") {
auto const t = 1_s;
Trajectory<Line> base(line, t);
- CHECK(line.GetPosition(t).GetCoordinates() ==
- base.GetPosition(1.).GetCoordinates());
+ CHECK(line.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(1));
}
@@ -183,8 +190,7 @@ TEST_CASE("Trajectories") {
auto const t = 1234_s;
Trajectory<Helix> const base(helix, t);
- CHECK(helix.GetPosition(t).GetCoordinates() ==
- base.GetPosition(1.).GetCoordinates());
+ CHECK(helix.GetPosition(t).GetCoordinates() == base.GetPosition(1.).GetCoordinates());
CHECK(base.ArcLength(1_s, 2_s) / 1_m == Approx(5));
}
diff --git a/Processes/TrackingLine/CMakeLists.txt b/Processes/TrackingLine/CMakeLists.txt
index 2911c9b78..83ed893c1 100644
--- a/Processes/TrackingLine/CMakeLists.txt
+++ b/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)
diff --git a/Processes/TrackingLine/TrackingLine.h b/Processes/TrackingLine/TrackingLine.h
index 602d54732..386215801 100644
--- a/Processes/TrackingLine/TrackingLine.h
+++ b/Processes/TrackingLine/TrackingLine.h
@@ -1,14 +1,20 @@
#ifndef _include_corsika_processes_TrackinLine_h_
#define _include_corsika_processes_TrackinLine_h_
-#include <corsika/geometry/Vector.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 <optional>
+
+#include <iostream>
+
using namespace corsika;
namespace corsika::process {
@@ -16,19 +22,54 @@ 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<corsika::units::si::TimeType> TimeOfIntersection(
+ geometry::Trajectory<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) {
+ auto const t = (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
+ return t;
+ } 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);
+ 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;
+ geometry::Line traj(p.GetPosition(), velocity);
return geometry::Trajectory<corsika::geometry::Line>(traj, 100_ns);
}
};
- } // namespace stack_inspector
+ } // namespace tracking_line
} // namespace corsika::process
diff --git a/Processes/TrackingLine/testTrackingLine.cc b/Processes/TrackingLine/testTrackingLine.cc
new file mode 100644
index 000000000..a07101096
--- /dev/null
+++ b/Processes/TrackingLine/testTrackingLine.cc
@@ -0,0 +1,58 @@
+/**
+ * (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;
+
+
+TEST_CASE("TrackingLine") {
+ corsika::environment::Environment env; // dummy environment
+ auto const& cs = env.GetCoordinateSystem();
+
+ tracking_line::TrackingLine<setup::Stack> 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());
+ REQUIRE(opt.value() / 9_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());
+ }
+}
--
GitLab