diff --git a/Environment/HomogeneousMedium.h b/Environment/HomogeneousMedium.h
new file mode 100644
index 0000000000000000000000000000000000000000..b1a465959741daedc3a2d27c9cf44986be47da2c
--- /dev/null
+++ b/Environment/HomogeneousMedium.h
@@ -0,0 +1,29 @@
+#ifndef _include_HomogeneousMedium_h_
+#define _include_HomogeneousMedium_h_
+
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/units/PhysicalUnits.h>
+
+/**
+ * a homogeneous medium
+ */
+
+namespace corsika::environment {
+
+ template <class T>
+ class HomogeneousMedium : T {
+ MassDensityType const fDensity;
+ NuclearComposition const fNuclComp;
+
+ public:
+ HomogeneousMedium(MassDensityType pDensity, NuclearComposition pNuclComp)
+ : fDensity(pDensity)
+ , fNuclComp(pNuclComp){};
+
+ MassDensityType GetMassDensity(Point const& p) const override { return fDensity; }
+ NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }
+ };
+
+} // namespace corsika::environment
+#endif
diff --git a/Environment/HydrogenSphere/HydrogenSphere.h b/Environment/HydrogenSphere/HydrogenSphere.h
deleted file mode 100644
index ca5b0f31439f0d35c82bc21e64b283f866df3758..0000000000000000000000000000000000000000
--- a/Environment/HydrogenSphere/HydrogenSphere.h
+++ /dev/null
@@ -1,41 +0,0 @@
-#ifndef _include_HydrogenSphere_h_
-#define _include_HydrogenSphere_h_
-
-#include <corsika/geometry/CoordinateSystem.h>
-#include <corsika/particles/ParticleProperties.h>
-#include <corsika/units/PhysicalUnits.h>
-
-/**
- * a fSphere homogeneously filled with hydrogen
- */
-
-namespace corsika::environment {
-
-class HydrogenSphere {
- CoordinateSystem const& fCS;
- corsika::geometry::Sphere const fSphere;
- MassDensityType const fDensity;
-
-public:
- HydrogenSphere(corsika::geometry::CoordinateSystem const& pEnvCS, LengthType pRadius,
- MassDensityType pDensity)
- : fCS(pEnvCS) fSphere(corsika::geometry::Point(pEnvCS, {0_m, 0_m, 0_m}), radius)
- , fDensity(pDensity) {}
-
- auto GetTargetParticle(Point const& p) const {
- return fSphere.isInside(p) ? corsika::particles::Code::Proton
- : corsika::particles::Code::Unknown;
- }
-
- MassDensityType GetDensity(Point const& p) const { return fSphere.isInside(p) ? density : 0_kg / (meter*meter*meter); };
-
- GetMagneticField(Point const& p) {
- QuantityVector<magnetic_flux_density_d> components{0 * corsika::units::tesla,
- 0. * corsika::units::tesla,
- 0. * corsika::units::tesla};
- return Vector<magnetic_flux_density_d>(fCS, components);
- }
-};
-
-}
-#endif
diff --git a/Environment/IMediumModel.h b/Environment/IMediumModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..f291246cb6da9ffe4fa139c72cc1af61a48fe40a
--- /dev/null
+++ b/Environment/IMediumModel.h
@@ -0,0 +1,25 @@
+#ifndef _include_IMediumModel_h
+#define _include_IMediumModels_h
+
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/geometry/BaseTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <tuple>
+#include <vector>
+
+namespace corsika::environment {
+
+ class IMediumModel {
+ public:
+ virtual ~IMediumModel() = default;
+
+ virtual corsika::units::si::MassDensityType GetMassDensity(
+ corsika::geometry::Point const&) const = 0;
+ virtual corsika::units::si::GrammageType IntegratedGrammage(BaseTrajectory const&,
+ double, double) const = 0;
+ virtual NuclearComposition const& GetNuclearComposition() const = 0;
+ };
+
+} // namespace corsika::environment
+
+#endif
diff --git a/Environment/NuclearComposition.h b/Environment/NuclearComposition.h
new file mode 100644
index 0000000000000000000000000000000000000000..0fda0ee44ee31acebc8d8f760a2d4a48ba8750b4
--- /dev/null
+++ b/Environment/NuclearComposition.h
@@ -0,0 +1,34 @@
+#ifndef _include_NuclearComposition_h
+#define _include_NuclearComposition_h
+
+#include <algorithm>
+#include <vector>
+
+namespace corsika::environment {
+ class NuclearComposition {
+ std::vector<float> const fNumberFractions; //<! relative fractions of number density
+ std::vector<corsika::particles::Code> const
+ fComponents; //!< particle codes of consitutents
+
+ public:
+ NuclearComposition(std::vector<corsika::particles::Code> pComponents,
+ std::vector<float> pFractions)
+ : fNumberFractions(pFractions)
+ , fComponents(pComponents) {
+ auto const sumFractions =
+ std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
+
+ if (!(0.999f < sum && sum < 1.001f)) {
+ throw std::string("element fractions do not add up to 1");
+ }
+ }
+
+ auto size() const { return fNumberFractions.size(); }
+
+ auto const& GetFractions() const { return fNumberFractions; }
+ auto const& GetComponents() const { return fComponents; }
+ };
+
+} // namespace corsika::environment
+
+#endif
diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h
index 2dcc0caf60c8456229ed29aafd02d633680a7037..7026ff06e33d5cfff661b440731bae70387da543 100644
--- a/Environment/VolumeTreeNode.h
+++ b/Environment/VolumeTreeNode.h
@@ -1,113 +1,114 @@
#ifndef _include_VolumeTreeNode_H
#define _include_VolumeTreeNode_H
+#include <corsika/geometry/Volume.h>
+#include <memory>
+#include <vector>
+
namespace corsika::environment {
-class Empty {}; //<! intended for usage as default template arguments
+ class Empty {}; //<! intended for usage as default template argument
-template <typename IModelProperties=Empty>
-class VolumeTreeNode {
-public:
+ 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;
+ using VolUPtr = std::unique_ptr<corsika::geometry::Volume>;
+
+ VolumeTreeNode(VolUPtr pVolume = nullptr)
+ : fGeoVolume(std::move(pVolume)) {}
+
+ bool Contains(corsika::geometry::Point const& p) const {
+ return fGeoVolume->Contains(p);
+ }
+
+ VolumeTreeNode<IModelProperties> const* Excludes(
+ corsika::geometry::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
+ VolumeTreeNode<IModelProperties> const* GetContainingNode(
+ corsika::geometry::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 (*childContainsIter)->GetContainingNode(p);
+ 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;
+ 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.
}
-
- auto const& GetModelProperties() const {
- return *fModelProperties;
+
+ 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;
+ 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)...);
+ 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)...));
+ 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:
+
+ private:
std::vector<VTNUPtr> fChildNodes;
std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
VolUPtr fGeoVolume;
- std::shared_ptr<IModelProperties> fModelProperties;
-};
+ std::shared_ptr<IModelProperties> fModelProperties;
+ };
-}
+} // namespace corsika::environment
#endif
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index 5e7368dc6e62f0ce14c190e4f62144e3f4305b82..515db9e08415de1b91ec4a65d018dc50b173f902 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -38,6 +38,7 @@ namespace corsika::units::si {
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>, double>;
} // end namespace corsika::units::si