diff --git a/Documentation/Examples/geometry_example.cc b/Documentation/Examples/geometry_example.cc
index 6a05974fd6d5e0224d809072b99979e29578f7d2..4091ab904c9dc7a8b233b21af78f2f45b0874ad9 100644
--- a/Documentation/Examples/geometry_example.cc
+++ b/Documentation/Examples/geometry_example.cc
@@ -39,10 +39,10 @@ int main() {
std::cout << "p2-p1 norm^2: " << norm << std::endl;
Sphere s(p1, 10_m); // define a sphere around a point with a radius
- std::cout << "p1 inside s: " << s.isInside(p2) << std::endl;
+ std::cout << "p1 inside s: " << s.Contains(p2) << std::endl;
Sphere s2(p1, 3_um); // another sphere
- std::cout << "p1 inside s2: " << s2.isInside(p2) << std::endl;
+ std::cout << "p1 inside s2: " << s2.Contains(p2) << std::endl;
// let's try parallel projections:
auto const v1 = Vector<length_d>(root, {1_m, 1_m, 0_m});
diff --git a/Environment/BaseEnvironment.h b/Environment/BaseEnvironment.h
deleted file mode 100644
index e85febecad782321c189145fb57e0fc9534e671f..0000000000000000000000000000000000000000
--- a/Environment/BaseEnvironment.h
+++ /dev/null
@@ -1,57 +0,0 @@
-#ifndef _include_EnvironmentBase_h
-#define _include_EnvironmentBase_h
-
-#include <vector>
-#include <memory>
-#include <algorithm>
-#include <corsika/geometry/Sphere.h>
-#include <corsika/geometry/Point.h>
-
-namespace corsika::environment {
-
-class BaseEnvironment {
- std::vector<std::unique_ptr<BaseEnvironment>> childNodes, excludedOverlapNodes;
- corsika::geometry::Sphere shape;
-
-public:
- bool ShapeContains(corsika::geometry::Point const& p) const {
- return shape.Contains(p);
- }
-
- BaseEnvironment* GetContainingNode(corsika::geometry::Point const& p) const {
- if (!shape.Contains(p))
- {
- return nullptr;
- }
-
- auto const predicate = [&] (auto const& s) {return bool(s->GetContainingNode(p));};
-
- auto const childFountIt = std::find_if(childNodes.cbegin(), childNodes.cend(), predicate);
-
- if (childFoundIt == childNodes.cend())
- {
- if (auto const excludedIt = std::find_if(excludedOverlapNodes.cbegin(), excludedOverlapNodes.cend(), predicate); excluded == excludedOverlapNodes.cend())
- {
- return *this;
- }
- else
- {
- return excludedIt;
- }
- }
-
- return (*childFoundIt)->GetContainingNode(p);
- }
-
- /**
- * Get mass density at given point \a p. It is the caller's responsibility
- * to call this method from the appropriate node in whose domain \a p is located.
- */
- virtual MassDensityType GetDensity(corsika::geometry::Point const& p) = 0;
-
- virtual std::pair<corsika::particles::Code const*, corsika::particles::Code const*> GetTargetComposition(Point const& p)
-};
-
-}
-
-#endif
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
new file mode 100644
index 0000000000000000000000000000000000000000..2dcc0caf60c8456229ed29aafd02d633680a7037
--- /dev/null
+++ b/Environment/VolumeTreeNode.h
@@ -0,0 +1,113 @@
+#ifndef _include_VolumeTreeNode_H
+#define _include_VolumeTreeNode_H
+
+namespace corsika::environment {
+
+class Empty {}; //<! intended for usage as default template arguments
+
+template <typename IModelProperties=Empty>
+class VolumeTreeNode {
+public:
+ using VTNUPtr = std::unique_ptr<VolumeTreeNode<IModelProperties>>;
+ using IMPSharedPtr = std::shared_ptr<IModelProperties>;
+
+ VolumeTreeNode(VolUPtr pVolume = nullptr) : fGeoVolume(std::move(pVolume)) { }
+
+ bool Contains(Point const& p) const { return fGeoVolume->Contains(p); }
+
+ VolumeTreeNode<IModelProperties> const* Excludes(Point const& p) const {
+ auto exclContainsIter = std::find_if(fExcludedNodes.cbegin(), fExcludedNodes.cend(),
+ [&] (auto const& s) {return bool(s->Contains(p));});
+
+ return exclContainsIter != fExcludedNodes.cend() ? *exclContainsIter : nullptr;
+ }
+
+ /** returns a pointer to the sub-VolumeTreeNode which is "responsible" for the given
+ * \class Point \arg p, or nullptr iff \arg p is not contained in this volume.
+ */
+ VolumeTreeNode<IModelProperties> const* GetContainingNode(Point const& p) const {
+ if (!Contains(p))
+ {
+ return nullptr;
+ }
+
+ if (auto const childContainsIter = std::find_if(fChildNodes.cbegin(), fChildNodes.cend(),
+ [&] (auto const& s) {return bool(s->Contains(p));});
+ childContainsIter == fChildNodes.cend()) // not contained in any of the children
+ {
+ if (auto const exclContainsIter = Excludes(p)) // contained in any excluded nodes
+ {
+ return exclContainsIter->GetContainingNode(p);
+ }
+ else
+ {
+ return this;
+ }
+ }
+ else
+ {
+ return (*childContainsIter)->GetContainingNode(p);
+ }
+ }
+
+ void addChild(VTNUPtr pChild) {
+ pChild->fParentNode = this;
+ fChildNodes.push_back(std::move(pChild));
+ // It is a bad idea to return an iterator to the inserted element
+ // because it might get invalidated when the vector needs to grow
+ // later and the caller won't notice.
+ }
+
+ void excludeOverlapWith(VTNUPtr const& pNode) {
+ fExcludedNodes.push_back(pNode.get());
+ }
+
+ auto GetParent() const {
+ return fParentNode;
+ };
+
+ auto const& GetVolume() const {
+ return *fGeoVolume;
+ }
+
+ auto const& GetModelProperties() const {
+ return *fModelProperties;
+ }
+
+ template <typename ModelProperties, typename... Args>
+ auto SetModelProperties(Args&&... args) {
+ static_assert(std::is_base_of_v<IModelProperties, ModelProperties>, "unusable type provided");
+
+ fModelProperties = std::make_unique<ModelProperties>(std::forward<Args>(args)...);
+ }
+
+ auto SetModelProperties(IMPSharedPtr ptr) {
+ fModelProperties = IMPSharedPtr;
+ }
+
+ template <class MediumType, typename... Args>
+ static auto CreateMedium(Args&&... args) {
+ static_assert(std::is_base_of_v<IMediumModel, MediumType>, "unusable type provided, needs to be derived from \"IMediumModel\"");
+
+ return std::make_shared<MediumType>(std::forward<Args>(args)...);
+ }
+
+ // factory methods for creation of nodes
+ template <class VolumeType, typename... Args>
+ static auto CreateNode(Args&&...args) {
+ static_assert(std::is_base_of_v<Volume, VolumeType>, "unusable type provided, needs to be derived from \"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;
+ VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
+ VolUPtr fGeoVolume;
+ std::shared_ptr<IModelProperties> fModelProperties;
+};
+
+}
+
+#endif
diff --git a/Framework/Geometry/CMakeLists.txt b/Framework/Geometry/CMakeLists.txt
index a4fc43cb3364881d7f3e42939b0985b69df4afab..d4aef5a42682b482596a2a347c247e5311bf6d37 100644
--- a/Framework/Geometry/CMakeLists.txt
+++ b/Framework/Geometry/CMakeLists.txt
@@ -9,6 +9,7 @@ set (
Vector.h
Point.h
Sphere.h
+ Volume.h
CoordinateSystem.h
Helix.h
BaseVector.h
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
index 28b4e7dc63d5ded1ecd6d047df442af4557e4f20..76699690c5d49b40a7e25d4f44e95cbead64b5e0 100644
--- a/Framework/Geometry/Sphere.h
+++ b/Framework/Geometry/Sphere.h
@@ -1,23 +1,24 @@
#ifndef _include_SPHERE_H_
#define _include_SPHERE_H_
+#include <corsika/geometry/Volume.h>
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
namespace corsika::geometry {
- class Sphere {
- Point center;
- LengthType const radius;
+ class Sphere : public Volume {
+ Point const fCenter;
+ LengthType const fRadius;
public:
Sphere(Point const& pCenter, LengthType const pRadius)
- : center(pCenter)
- , radius(pRadius) {}
+ : fCenter(pCenter)
+ , fRadius(pRadius) {}
//! returns true if the Point p is within the sphere
- auto isInside(Point const& p) const {
- return radius * radius > (center - p).squaredNorm();
+ bool Contains(Point const& p) const override {
+ return fRadius * fRadius > (fCenter - p).squaredNorm();
}
};
diff --git a/Framework/Geometry/Volume.h b/Framework/Geometry/Volume.h
new file mode 100644
index 0000000000000000000000000000000000000000..badbeb9df00ee72d77894ddd175b699e2fac19ae
--- /dev/null
+++ b/Framework/Geometry/Volume.h
@@ -0,0 +1,19 @@
+#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
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 3c7f63af445c847541fd2c3fc665aa8ae4d3198d..440930718e9e1e6c5a92599ed1d1680988beab8c 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -116,8 +116,8 @@ TEST_CASE("Sphere") {
Sphere sphere(center, 5_m);
SECTION("isInside") {
- REQUIRE_FALSE(sphere.isInside(Point(rootCS, {100_m, 0_m, 0_m})));
- REQUIRE(sphere.isInside(Point(rootCS, {2_m, 3_m, 4_m})));
+ REQUIRE_FALSE(sphere.Contains(Point(rootCS, {100_m, 0_m, 0_m})));
+ REQUIRE(sphere.Contains(Point(rootCS, {2_m, 3_m, 4_m})));
}
}