update to magnetic field model testing

corsika/detail/media/GeomagneticModel.inl

@@ -7,9 +7,6 @@
  */
 
 #include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-
-#include <boost/filesystem.hpp>
 
 #include <stdexcept>
 #include <string>
@@ -18,16 +15,15 @@
 namespace corsika {
 
   inline GeomagneticModel::GeomagneticModel(Point const& center,
-                                            std::string const& dataFile)
+                                            boost::filesystem::path const path)
       : center_(center) {
 
     // Read in coefficients
-    boost::filesystem::path const path = corsika::corsika_data(dataFile);
     boost::filesystem::ifstream file(path, std::ios::in);
 
     // Exit if file opening failed
     if (!file.is_open()) {
-      CORSIKA_LOG_ERROR("Failed opening data file {}", dataFile);
+      CORSIKA_LOG_ERROR("Failed opening data file {}", path);
       throw std::runtime_error("Cannot load GeomagneticModel data.");
     }
 
@@ -93,6 +89,11 @@ namespace corsika {
       }
     }
     file.close();
+
+    if (parameters_.size() == 0) {
+      CORSIKA_LOG_ERROR("No input data read!");
+      throw std::runtime_error("No input data read");
+    }
   }
 
   inline MagneticFieldVector GeomagneticModel::getField(double const year,
@@ -101,17 +102,14 @@ namespace corsika {
                                                         double const longitude) {
 
     int iYear = int(year);
-    int iEpoch = 0;
-    for (auto parIt = parameters_.rbegin(); parIt != parameters_.rend(); ++parIt) {
-      if (parIt->first <= iYear) {
-        iEpoch = parIt->first;
-        break;
-      }
+    auto iEpoch = parameters_.rbegin();
+    for (; iEpoch != parameters_.rend(); ++iEpoch) {
+      if (iEpoch->first <= iYear) { break; }
     }
-    double epoch = double(iEpoch);
-    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch, year);
-    if (iEpoch == 0) {
+    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch->first, year);
+    if (iEpoch == parameters_.rend()) {
       CORSIKA_LOG_WARN("Year {} is before first EPOCH. Results unclear.", year);
+      iEpoch--; // move one epoch back
     }
     if (altitude < -1_km || altitude > 600_km) {
       CORSIKA_LOG_WARN("Altitude should be between -1_km and 600_km.");
@@ -125,6 +123,7 @@ namespace corsika {
     if (longitude < -180 || longitude > 180) {
       CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
     }
+    double epoch = double(iEpoch->first);
 
     const double lat_geo = latitude * constants::pi / 180;
     const double lon = longitude * constants::pi / 180;
@@ -142,17 +141,23 @@ namespace corsika {
     double legendre, next_legendre, derivate_legendre;
     double magneticfield[3] = {0, 0, 0};
 
-    for (size_t j = 0; j < parameters_[iEpoch].size(); j++) {
+    for (size_t j = 0; j < iEpoch->second.size(); j++) {
 
-      ParameterLine p = parameters_[iEpoch][j];
+      ParameterLine p = iEpoch->second[j];
 
       // Time interpolation
-      p.g = p.g + (year - epoch) * p.dg;
-      p.h = p.h + (year - epoch) * p.dh;
-      if (iEpoch < 2020) {
-        ParameterLine next_p = parameters_[iEpoch + 5][j];
-        p.g = p.g + (next_p.g - p.g) * (year - epoch) / 5;
-        p.h = p.h + (next_p.h - p.h) * (year - epoch) / 5;
+      if (iEpoch == parameters_.rbegin() || p.dg != 0 || p.dh != 0) {
+        // this is the latest epoch in time, or time-dependence (dg/dh) was specified
+        // we use the extrapolation factors dg/dh:
+        p.g = p.g + (year - epoch) * p.dg;
+        p.h = p.h + (year - epoch) * p.dh;
+      } else {
+        // we linearly interpolate between two epochs
+        auto const nextEpoch = --iEpoch; // next epoch
+        ParameterLine const next_p = nextEpoch->second[j];
+        const double length = nextEpoch->first - epoch;
+        p.g = p.g + (next_p.g - p.g) * (year - epoch) / length;
+        p.h = p.h + (next_p.h - p.h) * (year - epoch) / length;
       }
 
       legendre = pow(-1, p.m) * std::assoc_legendre(p.n, p.m, sin(lat_sph));

corsika/media/GeomagneticModel.hpp

@@ -10,7 +10,9 @@
 
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/PhysicalGeometry.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
 
+#include <boost/filesystem.hpp>
 #include <fstream>
 #include <string>
 
@@ -44,7 +46,8 @@ namespace corsika {
      * @param center Center of Earth.
      * @param data Data table to read.
      */
-    GeomagneticModel(Point const& center, std::string const& data = "GeoMag/WMM.COF");
+    GeomagneticModel(Point const& center, boost::filesystem::path const path =
+                                              corsika::corsika_data("GeoMag/WMM.COF"));
 
     /**
      * Calculates the value of the magnetic field.

tests/media/CMakeLists.txt

@@ -12,3 +12,9 @@ set (test_media_sources
 CORSIKA_ADD_TEST (testMedia SOURCES ${test_media_sources})
 
 target_link_libraries (testMedia CONAN_PKG::boost) # for math/quadrature/gauss_kronrod.hpp
+
+target_compile_definitions (
+  testMedia
+  PRIVATE
+  REFDATADIR="${CMAKE_CURRENT_SOURCE_DIR}"
+)

tests/media/EmptyGeomagneticData.COF

+# no data here...
\ No newline at end of file

tests/media/testMagneticField.cpp

@@ -9,6 +9,7 @@
 #include <corsika/framework/core/ParticleProperties.hpp>
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
 #include <corsika/media/HomogeneousMedium.hpp>
 #include <corsika/media/IMediumModel.hpp>
 #include <corsika/media/UniformMagneticField.hpp>
@@ -19,10 +20,13 @@
 #include <SetupTestTrajectory.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>
 
+#include <boost/filesystem.hpp>
 #include <catch2/catch.hpp>
 
 using namespace corsika;
 
+const std::string refDataDir = std::string(REFDATADIR); // from cmake
+
 TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
 
   logging::set_level(logging::level::trace);
@@ -75,15 +79,19 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
 
   SECTION("GeomagneticModel") {
 
-    CHECK_THROWS(GeomagneticModel(gOrigin, "GeoMag/IDoNotExist.COF"));
+    // non existing data file
+    CHECK_THROWS(GeomagneticModel(gOrigin, corsika_data("GeoMag/IDoNotExist.COF")));
+    // and empty data file
+    CHECK_THROWS(GeomagneticModel(gOrigin, refDataDir + "/EmptyGeomagneticData.COF"));
 
     {
-      GeomagneticModel wmm(gOrigin, "GeoMag/WMM.COF");
+      GeomagneticModel wmm(gOrigin, corsika_data("GeoMag/WMM.COF"));
 
       // create earth magnetic field vector
       MagneticFieldVector Earth_B_1 = wmm.getField(2022.5, 100_km, -80, -120);
       MagneticFieldVector Earth_B_2 = wmm.getField(2022.5, 0_km, 0, 120);
       MagneticFieldVector Earth_B_3 = wmm.getField(2020, 100_km, 80, 0);
+      MagneticFieldVector Earth_B_4 = wmm.getField(2019.999, 100_km, 80, 0);
       CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5815).margin(0.05));
       CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(-14803).margin(0.05));
       CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.05));
@@ -93,9 +101,12 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
       CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6261.8).margin(0.05));
       CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(185.5).margin(0.05));
       CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.05));
+      CHECK(Earth_B_4.getX(gCS) / 1_nT == Approx(6261.8).margin(0.05));
+      CHECK(Earth_B_4.getY(gCS) / 1_nT == Approx(185.5).margin(0.05));
+      CHECK(Earth_B_4.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.05));
     }
     {
-      GeomagneticModel igrf(gOrigin, "GeoMag/IGRF13.COF");
+      GeomagneticModel igrf(gOrigin, corsika_data("GeoMag/IGRF13.COF"));
 
       // create earth magnetic field vector
       MagneticFieldVector Earth_B_1 = igrf.getField(2022.5, 100_km, -80, -120);