From 48e18eaec5221fa1b001a2138a873a0555f3a5a2 Mon Sep 17 00:00:00 2001
From: Philipp Windischhofer <windischhofer@uchicago.edu>
Date: Tue, 2 Jul 2024 09:13:08 -0500
Subject: [PATCH] test pattern for environment dump

---
 corsika/media/HomogeneousMedium.hpp      |  5 ++
 corsika/media/IMagneticFieldModel.hpp    |  9 +++
 corsika/media/IMediumModel.hpp           |  7 ++
 corsika/media/IMediumPropertyModel.hpp   |  9 +++
 corsika/media/IRefractiveIndexModel.hpp  |  8 ++
 corsika/media/MediumPropertyModel.hpp    |  5 ++
 corsika/media/UniformMagneticField.hpp   |  5 ++
 corsika/media/UniformRefractiveIndex.hpp |  5 ++
 tests/media/CMakeLists.txt               |  3 +
 tests/media/testEnvironmentDump.cpp      | 95 ++++++++++++++++++++++++
 10 files changed, 151 insertions(+)
 create mode 100644 tests/media/testEnvironmentDump.cpp

diff --git a/corsika/media/HomogeneousMedium.hpp b/corsika/media/HomogeneousMedium.hpp
index d41c0bffa..7a8bd598b 100644
--- a/corsika/media/HomogeneousMedium.hpp
+++ b/corsika/media/HomogeneousMedium.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <iostream>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/Line.hpp>
 #include <corsika/framework/geometry/Point.hpp>
@@ -35,6 +36,10 @@ namespace corsika {
     LengthType getArclengthFromGrammage(BaseTrajectory const&,
                                         GrammageType grammage) const override;
 
+    void DescribeMassDensity(std::ostream& stream) const override {
+      stream << "mass_density = " << density_.magnitude() << std::endl;
+    }
+    
   private:
     MassDensityType const density_;
     NuclearComposition const nuclComp_;
diff --git a/corsika/media/IMagneticFieldModel.hpp b/corsika/media/IMagneticFieldModel.hpp
index 3efd34546..9b19140d4 100644
--- a/corsika/media/IMagneticFieldModel.hpp
+++ b/corsika/media/IMagneticFieldModel.hpp
@@ -8,6 +8,8 @@
 
 #pragma once
 
+#include <iostream>
+
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
@@ -40,6 +42,13 @@ namespace corsika {
      */
     virtual ~IMagneticFieldModel() = default; // LCOV_EXCL_LINE
 
+    virtual void DescribeMagneticField(std::ostream&) const { }; // = 0;   // make pure-virtual once implemented everywhere
+    
+    virtual void AddProperties(std::ostream& stream) const override {
+      Model::AddProperties(stream);
+      DescribeMagneticField(stream);
+    };
+
   }; // END: class MagneticField
 
 } // namespace corsika
diff --git a/corsika/media/IMediumModel.hpp b/corsika/media/IMediumModel.hpp
index e3e030392..e77e4d4eb 100644
--- a/corsika/media/IMediumModel.hpp
+++ b/corsika/media/IMediumModel.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <iostream>
 #include <corsika/media/NuclearComposition.hpp>
 #include <corsika/framework/geometry/Point.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
@@ -44,6 +45,12 @@ namespace corsika {
                                                 GrammageType) const = 0;
 
     virtual NuclearComposition const& getNuclearComposition() const = 0;
+
+    virtual void DescribeMassDensity(std::ostream&) const { }; // = 0;  // make pure-virtual once implemented everywhere
+    
+    virtual void AddProperties(std::ostream& stream) const {
+      DescribeMassDensity(stream);
+    };
   };
 
 } // namespace corsika
diff --git a/corsika/media/IMediumPropertyModel.hpp b/corsika/media/IMediumPropertyModel.hpp
index 94e80fb51..7e6519230 100644
--- a/corsika/media/IMediumPropertyModel.hpp
+++ b/corsika/media/IMediumPropertyModel.hpp
@@ -8,6 +8,8 @@
 
 #pragma once
 
+#include <iostream>
+
 #include <corsika/media/MediumProperties.hpp>
 
 #include <corsika/framework/geometry/Point.hpp>
@@ -37,6 +39,13 @@ namespace corsika {
      */
     virtual ~IMediumPropertyModel() = default;
 
+    virtual void DescribeMedium(std::ostream&) const { }; // = 0;   // make pure-virtual once implemented everywhere
+    
+    virtual void AddProperties(std::ostream& stream) const override {
+      TModel::AddProperties(stream);
+      DescribeMedium(stream);
+    };
+
   }; // END: class IMediumTypeModel
 
 } // namespace corsika
diff --git a/corsika/media/IRefractiveIndexModel.hpp b/corsika/media/IRefractiveIndexModel.hpp
index b042540fb..f14095959 100644
--- a/corsika/media/IRefractiveIndexModel.hpp
+++ b/corsika/media/IRefractiveIndexModel.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <iostream>
 #include <corsika/framework/geometry/Point.hpp>
 
 namespace corsika {
@@ -35,6 +36,13 @@ namespace corsika {
      */
     virtual ~IRefractiveIndexModel() = default;
 
+    virtual void DescribeRefractiveIndex(std::ostream&) const { }; // = 0;   // make pure-virtual once implemented everywhere
+
+    virtual void AddProperties(std::ostream& stream) const override {
+      TModel::AddProperties(stream);
+      DescribeRefractiveIndex(stream);
+    };
+    
   }; // END: class IRefractiveIndexModel
 
 } // namespace corsika
diff --git a/corsika/media/MediumPropertyModel.hpp b/corsika/media/MediumPropertyModel.hpp
index 97bce6ff8..81f704b4a 100644
--- a/corsika/media/MediumPropertyModel.hpp
+++ b/corsika/media/MediumPropertyModel.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <iostream>
 #include <corsika/media/IMediumPropertyModel.hpp>
 
 namespace corsika {
@@ -45,6 +46,10 @@ namespace corsika {
      */
     void setMedium(Medium const medium);
 
+    void DescribeMedium(std::ostream& stream) const override {
+      stream << "A MediumPropertyModel with some attributes " << std::endl;
+    };
+    
   }; // END: class MediumPropertyModel
 
 } // namespace corsika
diff --git a/corsika/media/UniformMagneticField.hpp b/corsika/media/UniformMagneticField.hpp
index 02f45304a..2e0873eee 100644
--- a/corsika/media/UniformMagneticField.hpp
+++ b/corsika/media/UniformMagneticField.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <iostream>
 #include <corsika/media/IMagneticFieldModel.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>
@@ -55,6 +56,10 @@ namespace corsika {
      */
     void setMagneticField(MagneticFieldVector const& Bfield) { B_ = Bfield; }
 
+    void DescribeMagneticField(std::ostream& stream) const override {
+      stream << "B-field = " << B_.getNorm().magnitude() << std::endl;
+    }
+    
   private:
     MagneticFieldVector B_; ///< The constant magnetic field we use.
 
diff --git a/corsika/media/UniformRefractiveIndex.hpp b/corsika/media/UniformRefractiveIndex.hpp
index c5bcebd6c..b38d4c16f 100644
--- a/corsika/media/UniformRefractiveIndex.hpp
+++ b/corsika/media/UniformRefractiveIndex.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include <iostream>
 #include <corsika/media/IRefractiveIndexModel.hpp>
 
 namespace corsika {
@@ -51,6 +52,10 @@ namespace corsika {
      */
     void setRefractiveIndex(double const n);
 
+    void DescribeRefractiveIndex(std::ostream& stream) const override {
+      stream << "IOR = " << n_ << std::endl;
+    } 
+
   }; // END: class RefractiveIndex
 
 } // namespace corsika
diff --git a/tests/media/CMakeLists.txt b/tests/media/CMakeLists.txt
index 5a8442123..e886d98a2 100644
--- a/tests/media/CMakeLists.txt
+++ b/tests/media/CMakeLists.txt
@@ -18,3 +18,6 @@ target_compile_definitions (
   PRIVATE
   REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
 )
+
+add_executable (testEnvironmentDump testEnvironmentDump.cpp)
+target_link_libraries (testEnvironmentDump CORSIKA8)
diff --git a/tests/media/testEnvironmentDump.cpp b/tests/media/testEnvironmentDump.cpp
new file mode 100644
index 000000000..46cf7e3b1
--- /dev/null
+++ b/tests/media/testEnvironmentDump.cpp
@@ -0,0 +1,95 @@
+#include <corsika/framework/geometry/SeparationPlane.hpp>
+#include <corsika/media/Environment.hpp>
+#include <corsika/media/HomogeneousMedium.hpp>
+#include <corsika/media/IMagneticFieldModel.hpp>
+#include <corsika/media/MediumPropertyModel.hpp>
+#include <corsika/media/NuclearComposition.hpp>
+#include <corsika/media/ShowerAxis.hpp>
+#include <corsika/media/UniformMagneticField.hpp>
+#include <corsika/media/UniformRefractiveIndex.hpp>
+
+using namespace corsika;
+using namespace std;
+
+using EnvironmentInterface =
+    IRefractiveIndexModel<IMediumPropertyModel<IMagneticFieldModel<IMediumModel>>>;
+using EnvType = Environment<EnvironmentInterface>;
+
+int main(int argc, char** argv) {
+
+  EnvType env;
+  CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
+
+  Point const surfaceCenter{rootCS, 0_m, 0_m, 0_m};
+  auto iceBlock = EnvType::createNode<SeparationPlane>(
+      Plane(surfaceCenter, DirectionVector(rootCS, {0, 0, 1})));
+
+  auto airBlock = EnvType::createNode<SeparationPlane>(
+      Plane(surfaceCenter, DirectionVector(rootCS, {0, 0, -1})));
+  
+  // see https://www.ngdc.noaa.gov/geomag/calculators/magcalc.shtml#igrfwmm
+  // for 90 deg south, 1.95km alt.
+  MagneticFieldVector const southPoleBField(rootCS, -8.734_uT, 14.357_uT, 51.694_uT);
+
+  // Creating the universe where all mediums exist.
+  using NodeT = VolumeTreeNode<EnvironmentInterface>;
+  NodeT* const universe = env.getUniverse().get();
+  
+  // Ice properties
+  NuclearComposition const nuclearCompositionIce{{Code::Hydrogen, Code::Oxygen},
+                                                 {2.0 / 3.0, 1.0 / 3.0}};
+  MassDensityType const iceDensity = 0.919_g / (1_cm * 1_cm * 1_cm);
+  double const iceRefractiveIndex = 1.42;
+
+  auto const mediumIce = std::make_shared<UniformRefractiveIndex<MediumPropertyModel<
+      UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>>>(
+      iceRefractiveIndex, Medium::WaterIce, southPoleBField, iceDensity,
+      nuclearCompositionIce);
+
+  iceBlock->setModelProperties(mediumIce);
+  universe->addChild(std::move(iceBlock));
+
+  // Air properties
+  // TODO: The following is just a guess, need to use actual air
+  NuclearComposition const nuclearCompositionAir{{Code::Nitrogen, Code::Oxygen},
+						 {3.0 / 4.0, 1.0 / 4.0}};
+  MassDensityType const airDensity = 1.2e-3_g / (1_cm * 1_cm * 1_cm);
+  double const airRefractiveIndex = 1.00001;
+
+  auto const mediumAir = std::make_shared<UniformRefractiveIndex<MediumPropertyModel<
+    UniformMagneticField<HomogeneousMedium<EnvironmentInterface>>>>>(
+      airRefractiveIndex, Medium::AirDry1Atm, southPoleBField, airDensity,
+      nuclearCompositionAir);
+
+  airBlock->setModelProperties(mediumAir);
+  universe->addChild(std::move(airBlock));
+
+  auto dumper = [](const NodeT& node){
+    const EnvironmentInterface* props = &node.getModelProperties();
+    if(props != nullptr) {
+      std::cout << "------" << std::endl;      
+      props -> AddProperties(std::cout);      
+      std::cout << "------" << std::endl;
+    }
+  };
+  universe->walk(dumper);
+
+  // Output:
+  
+  // ------
+  // mass_density = 919
+  // B-field = 5.43569e-05
+  // A MediumPropertyModel with some attributes
+  // IOR = 1.42
+  // ------
+  // ------
+  // mass_density = 1.2
+  // B-field = 5.43569e-05
+  // A MediumPropertyModel with some attributes
+  // IOR = 1.00001
+  // ------
+  
+  std::cout << "done" << std::endl;
+  
+  return 0;
+}
-- 
GitLab