medium interface

Environment/HomogeneousMedium.h

+#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

Environment/HydrogenSphere/HydrogenSphere.h

-#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

Environment/IMediumModel.h

+#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

Environment/NuclearComposition.h

+#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

Environment/VolumeTreeNode.h

 #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

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