Merge branch '453-simplified-media-constructors-for-common-media' into 'master'

Resolve "Simplified media constructors for common media"

Closes #453

See merge request AirShowerPhysics/corsika!425

examples/CMakeLists.txt

@@ -64,3 +64,5 @@ CORSIKA_REGISTER_EXAMPLE (corsika RUN_OPTIONS --energy=100 --zenith=10 --nevent=
 add_executable (mars mars.cpp)
 target_link_libraries (mars CORSIKA8::CORSIKA8)
 
+add_executable (environment environment.cpp)
+target_link_libraries (environment CORSIKA8::CORSIKA8)

examples/environment.cpp

+/*
+ * (c) Copyright 2022 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * 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/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMediumModel.hpp>
+#include <corsika/media/MediumProperties.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+
+using namespace corsika;
+
+// Example to show how to construct an environment in CORSIKA
+
+// MyMediumInterface is constructed via "mixin inheritance"
+// Here, we construct our base class with IMediumModel and IMediumPropertyModel:
+// the former allow us define a density profile and nuclear composite,
+// the later is used to determine energy losses parameters.
+using MyMediumInterface = IMediumPropertyModel<IMediumModel>;
+
+using MyEnvType = Environment<MyMediumInterface>;
+
+int main() {
+  // create environment and universe, empty so far
+  MyEnvType env;
+
+  // create a geometry object: a sphere with radius of 1000 m
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+  Point const center{rootCS, 0_m, 0_m, 0_m};
+  LengthType const radius = 1000_m;
+  auto sphere = std::make_unique<Sphere>(center, radius);
+
+  // create node from geometry object
+  auto node = std::make_unique<VolumeTreeNode<MyMediumInterface>>(std::move(sphere));
+
+  // set medium properties to our node, say it is water
+  NuclearComposition const nucl_comp{{Code::Hydrogen, Code::Oxygen}, {0.11, 0.89}};
+  MassDensityType const density = 1_g / (1_cm * 1_cm * 1_cm);
+
+  // create concrete implementation of MyMediumInterface by combining parts that
+  // implement the corresponding interfaces. HomogeneousMedium implements IMediumModel,
+  // MediumPropertyModel implements IMediumPropertyModel.
+  auto const medium =
+      std::make_shared<MediumPropertyModel<HomogeneousMedium<MyMediumInterface>>>(
+          Medium::WaterLiquid, density, nucl_comp);
+  node->setModelProperties(medium);
+
+  // put our node into universe
+  // node: this has to be down after setting node model properties, since
+  // std::move will make our previous defined node, which is a unique pointer
+  // un-referencable in the context
+  VolumeTreeNode<MyMediumInterface>* const universe = env.getUniverse().get();
+  universe->addChild(std::move(node));
+
+  // example to explore the media properties of the node
+  for (auto h = 0_m; h < radius; h += 100_m) {
+    Point const ptest{rootCS, 0_m, 0_m, h};
+    MyMediumInterface const& media_prop =
+        env.getUniverse()->getContainingNode(ptest)->getModelProperties();
+    MassDensityType const rho = media_prop.getMassDensity(ptest);
+    NuclearComposition const& nuc_comp = media_prop.getNuclearComposition();
+    std::vector<Code> const& nuclei = nuc_comp.getComponents();
+    std::vector<double> const& fractions = nuc_comp.getFractions();
+    CORSIKA_LOG_INFO(
+        "radius: {:.2f} m, density: {:.2f} g/cm^3, nuclei: ({:d}, {:d}), fractions: "
+        "({:.2f}, "
+        "{:.2f})",
+        h / 1_m, rho / (1_g / static_pow<3>(1_cm)), get_PDG(nuclei.at(0)),
+        get_PDG(nuclei.at(1)), fractions.at(0), fractions.at(1));
+  }
+
+  return 0;
+}