From 956d5ae4cb7e6364f27ca1947cb0be8e17dba2b8 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Tue, 21 Dec 2021 11:17:00 +0100
Subject: [PATCH] update to magnetic field model testing
---
corsika/detail/media/GeomagneticModel.inl | 51 +++++++++++++----------
corsika/media/GeomagneticModel.hpp | 5 ++-
tests/media/CMakeLists.txt | 6 +++
tests/media/EmptyGeomagneticData.COF | 1 +
tests/media/testMagneticField.cpp | 17 ++++++--
5 files changed, 53 insertions(+), 27 deletions(-)
create mode 100644 tests/media/EmptyGeomagneticData.COF
diff --git a/corsika/detail/media/GeomagneticModel.inl b/corsika/detail/media/GeomagneticModel.inl
index 910792b7f..9e15b454d 100644
--- a/corsika/detail/media/GeomagneticModel.inl
+++ b/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));
diff --git a/corsika/media/GeomagneticModel.hpp b/corsika/media/GeomagneticModel.hpp
index 0c3d76096..bb9b3b923 100644
--- a/corsika/media/GeomagneticModel.hpp
+++ b/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.
diff --git a/tests/media/CMakeLists.txt b/tests/media/CMakeLists.txt
index c1265193c..5a8442123 100644
--- a/tests/media/CMakeLists.txt
+++ b/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}"
+)
diff --git a/tests/media/EmptyGeomagneticData.COF b/tests/media/EmptyGeomagneticData.COF
new file mode 100644
index 000000000..ba3fabdf8
--- /dev/null
+++ b/tests/media/EmptyGeomagneticData.COF
@@ -0,0 +1 @@
+# no data here...
\ No newline at end of file
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index f46297fa6..3d311e64b 100644
--- a/tests/media/testMagneticField.cpp
+++ b/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);
--
GitLab