homogeneous environment

Documentation/Examples/cascade_example.cc

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -18,6 +17,8 @@
 #include <corsika/setup/SetupTrajectory.h>
 
 #include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
 
 #include <corsika/random/RNGManager.h>
 
@@ -29,6 +30,7 @@
 #include <corsika/units/PhysicalUnits.h>
 
 #include <iostream>
+#include <limits>
 #include <typeinfo>
 
 using namespace corsika;
@@ -192,9 +194,8 @@ public:
     //    rnd_ini_();
 
     corsika::random::RNGManager& rmng = corsika::random::RNGManager::GetInstance();
-    ;
-    const std::string str_name = "s_rndm";
-    rmng.RegisterRandomStream(str_name);
+
+    rmng.RegisterRandomStream("s_rndm");
 
     // //    corsika::random::RNG srng;
     // auto srng = rmng.GetRandomStream("s_rndm");
@@ -242,13 +243,25 @@ double s_rndm_(int&) {
   return rmng() / (double)rmng.max();
 }
 
+class A {};
+class B : public A {};
+
 int main() {
-  corsika::environment::Environment env; // dummy environment  
+  corsika::environment::Environment env; // dummy environment
   auto& universe = *(env.GetUniverse());
 
   auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
-      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 100_km);
-      
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m},
+      1_km * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel =
+      corsika::environment::HomogeneousMedium<corsika::environment::IMediumModel>;
+  theMedium->SetModelProperties<MyHomogeneousModel>(
+      1_g / (1_m * 1_m * 1_m),
+      corsika::environment::NuclearComposition(
+          std::vector<corsika::particles::Code>{corsika::particles::Code::Proton},
+          std::vector<float>{1.}));
+
   universe.AddChild(std::move(theMedium));
 
   tracking_line::TrackingLine<setup::Stack> tracking(env);
@@ -261,7 +274,7 @@ int main() {
 
   stack.Clear();
   auto particle = stack.NewParticle();
-  EnergyType E0 = 1_EeV;
+  EnergyType E0 = 100_TeV;
   particle.SetEnergy(E0);
   particle.SetPID(Code::Proton);
   EAS.Init();

Environment/Environment.h

@@ -18,14 +18,16 @@
 
 namespace corsika::environment {
   struct Universe : public corsika::geometry::Volume {
-      bool Contains(corsika::geometry::Point const&) const override {return true;}
+    bool Contains(corsika::geometry::Point const&) const override { return true; }
   };
 
+  // template <typename IEnvironmentModel>
   class Environment {
   public:
-    Environment() : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
-          std::make_unique<Universe>())) {}
-  
+    Environment()
+        : fUniverse(std::make_unique<VolumeTreeNode<IEnvironmentModel>>(
+              std::make_unique<Universe>())) {}
+
     using IEnvironmentModel = corsika::setup::IEnvironmentModel;
 
     auto& GetUniverse() { return fUniverse; }

Environment/HomogeneousMedium.h

@@ -25,7 +25,7 @@
 namespace corsika::environment {
 
   template <class T>
-  class HomogeneousMedium : T {
+  class HomogeneousMedium : public T {
     corsika::units::si::MassDensityType const fDensity;
     NuclearComposition const fNuclComp;
 
@@ -35,8 +35,8 @@ namespace corsika::environment {
         : fDensity(pDensity)
         , fNuclComp(pNuclComp){};
 
-    corsika::units::si::MassDensityType GetMassDensity([
-        [maybe_unused]] corsika::geometry::Point const& p) const override {
+    corsika::units::si::MassDensityType GetMassDensity(
+        corsika::geometry::Point const&) const override {
       return fDensity;
     }
     NuclearComposition const& GetNuclearComposition() const override { return fNuclComp; }

Environment/NuclearComposition.h

@@ -3,6 +3,7 @@
 
 #include <corsika/particles/ParticleProperties.h>
 #include <numeric>
+#include <stdexcept>
 #include <vector>
 
 namespace corsika::environment {
@@ -20,7 +21,7 @@ namespace corsika::environment {
           std::accumulate(pFractions.cbegin(), pFractions.cend(), 0.f);
 
       if (!(0.999f < sumFractions && sumFractions < 1.001f)) {
-        throw std::string("element fractions do not add up to 1");
+        throw std::runtime_error("element fractions do not add up to 1");
       }
     }
 

Environment/VolumeTreeNode.h

@@ -87,12 +87,12 @@ namespace corsika::environment {
 
     auto const& GetModelProperties() const { return *fModelProperties; }
 
-    template <typename ModelProperties, typename... Args>
+    template <typename TModelProperties, typename... Args>
     auto SetModelProperties(Args&&... args) {
-      static_assert(std::is_base_of_v<IModelProperties, ModelProperties>,
+      static_assert(std::is_base_of_v<IModelProperties, TModelProperties>,
                     "unusable type provided");
 
-      fModelProperties = std::make_shared<ModelProperties>(std::forward<Args>(args)...);
+      fModelProperties = std::make_shared<TModelProperties>(std::forward<Args>(args)...);
       return fModelProperties;
     }
 
@@ -111,7 +111,7 @@ namespace corsika::environment {
     std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes;
     VolumeTreeNode<IModelProperties> const* fParentNode = nullptr;
     VolUPtr fGeoVolume;
-    std::shared_ptr<IModelProperties> fModelProperties;
+    IMPSharedPtr fModelProperties;
   };
 
 } // namespace corsika::environment

Environment/testEnvironment.cc

+/**
+ * (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.
+ */
+
 #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
                           // cpp file
 

Framework/Cascade/testCascade.cc

-
 /**
  * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
  *
@@ -9,6 +8,8 @@
  * the license.
  */
 
+#include <limits>
+
 #include <corsika/environment/Environment.h>
 
 #include <corsika/cascade/Cascade.h>
@@ -80,12 +81,13 @@ private:
 TEST_CASE("Cascade", "[Cascade]") {
   corsika::environment::Environment env; // dummy environment
   auto& universe = *(env.GetUniverse());
-  auto const radius = 20_km;
+  auto const radius = 1_m * std::numeric_limits<double>::infinity();
+  ;
   auto theMedium = corsika::environment::Environment::CreateNode<Sphere>(
-    Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);      
+      Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, radius);
   universe.AddChild(std::move(theMedium));
 
-  tracking_line::TrackingLine<setup::Stack> tracking(env);  
+  tracking_line::TrackingLine<setup::Stack> tracking(env);
 
   stack_inspector::StackInspector<setup::Stack> p0(true);
   ProcessSplit p1;

Framework/Geometry/CoordinateSystem.cc

@@ -10,6 +10,7 @@
  */
 
 #include <corsika/geometry/CoordinateSystem.h>
+#include <stdexcept>
 
 using namespace corsika::geometry;
 
@@ -40,7 +41,7 @@ EigenTransform CoordinateSystem::GetTransformation(CoordinateSystem const& pFrom
     commonBase = a;
 
   } else {
-    throw std::string("no connection between coordinate systems found!");
+    throw std::runtime_error("no connection between coordinate systems found!");
   }
 
   EigenTransform t = EigenTransform::Identity();

Framework/Geometry/CoordinateSystem.h

@@ -15,6 +15,7 @@
 #include <corsika/geometry/QuantityVector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <Eigen/Dense>
+#include <stdexcept>
 
 typedef Eigen::Transform<double, 3, Eigen::Affine> EigenTransform;
 typedef Eigen::Translation<double, 3> EigenTranslation;
@@ -60,7 +61,7 @@ namespace corsika::geometry {
 
     auto rotate(QuantityVector<phys::units::length_d> axis, double angle) const {
       if (axis.eVector.isZero()) {
-        throw std::string("null-vector given as axis parameter");
+        throw std::runtime_error("null-vector given as axis parameter");
       }
 
       EigenTransform const rotation{Eigen::AngleAxisd(angle, axis.eVector.normalized())};
@@ -71,7 +72,7 @@ namespace corsika::geometry {
     auto translateAndRotate(QuantityVector<phys::units::length_d> translation,
                             QuantityVector<phys::units::length_d> axis, double angle) {
       if (axis.eVector.isZero()) {
-        throw std::string("null-vector given as axis parameter");
+        throw std::runtime_error("null-vector given as axis parameter");
       }
 
       EigenTransform const transf{Eigen::AngleAxisd(angle, axis.eVector.normalized()) *

Processes/TrackingLine/TrackingLine.h

@@ -24,6 +24,7 @@
 #include <algorithm>
 #include <iostream>
 #include <optional>
+#include <stdexcept>
 #include <utility>
 
 using namespace corsika;
@@ -39,8 +40,9 @@ namespace corsika::process {
       corsika::environment::Environment const& fEnvironment;
 
     public:
-      std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>> TimeOfIntersection(
-          corsika::geometry::Line const& traj, geometry::Sphere const& sphere) {
+      std::optional<std::pair<corsika::units::si::TimeType, corsika::units::si::TimeType>>
+      TimeOfIntersection(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);
@@ -61,7 +63,8 @@ namespace corsika::process {
 
         if (discriminant.magnitude() > 0) {
           (-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm());
-          return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()), (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
+          return std::make_pair((-alpha - sqrt(discriminant)) / (2 * v0.squaredNorm()),
+                                (-alpha + sqrt(discriminant)) / (2 * v0.squaredNorm()));
         } else {
           return {};
         }
@@ -92,13 +95,12 @@ namespace corsika::process {
               volume); // for the moment we are a bit bold here and assume
                        // everything is a sphere, crashes with exception if not
 
-          if (auto opt = TimeOfIntersection(line, sphere); opt.has_value())
-          {
+          if (auto opt = TimeOfIntersection(line, sphere); opt.has_value()) {
             auto const [t1, t2] = *opt;
             if (t1.magnitude() >= 0)
-                intersectionTimes.push_back(t1);
+              intersectionTimes.push_back(t1);
             else if (t2.magnitude() >= 0)
-                intersectionTimes.push_back(t2);
+              intersectionTimes.push_back(t2);
           }
         };
 
@@ -111,11 +113,17 @@ namespace corsika::process {
         auto const minIter =
             std::min_element(intersectionTimes.cbegin(), intersectionTimes.cend());
 
+        TimeType min;
+
         if (minIter == intersectionTimes.cend()) {
-          throw std::string("no intersection with anything!");
+          min = 1_s; // todo: do sth. more reasonable as soon as tracking is able
+                     // to handle the numerics properly
+          //~ throw std::runtime_error("no intersection with anything!");
+        } else {
+          min = *minIter;
         }
 
-        return geometry::Trajectory<corsika::geometry::Line>(line, *minIter);
+        return geometry::Trajectory<corsika::geometry::Line>(line, min);
       }
     };