diff --git a/corsika/detail/media/WMM.inl b/corsika/detail/media/WMM.inl
deleted file mode 100644
index 3eecac4800bd830af5fceadfd8ce43d476d384d9..0000000000000000000000000000000000000000
--- a/corsika/detail/media/WMM.inl
+++ /dev/null
@@ -1,121 +0,0 @@
-/*
- * (c) Copyright 2020 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.
- */
-
-#pragma once
-
-#include <boost/math/special_functions/factorials.hpp>
-#include <boost/math/special_functions/legendre.hpp>
-
-#include <cmath>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-#include <fstream>
-
-namespace corsika {
- inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year,
- const LengthType altitude, const double latitude,
- const double longitude) {
- if (year < 2020 || year > 2025) {
- CORSIKA_LOG_WARN("Year has to be between 2020 and 2025.");
- }
- if (altitude < -1_km || altitude > 850_km) {
- CORSIKA_LOG_WARN("Altitude should be between -1_km and 850_km.");
- }
- if (latitude < -90 || latitude > 90) {
- CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
- abort();
- } else if (latitude < -89.992 || latitude > 89.992) {
- CORSIKA_LOG_WARN("Latitude is close to the poles.");
- }
- if (longitude < -180 || longitude > 180) {
- CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
- }
-
- const double lat_geo = latitude * constants::pi / 180;
- const double lon = longitude * constants::pi / 180;
-
- // Transform into spherical coordinates
- const double f = 1 / 298.257223563;
- const double e_squared = f * (2 - f);
- LengthType R_c =
- constants::EarthRadius::Equatorial / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
- LengthType p = (R_c + altitude) * cos(lat_geo);
- LengthType z = sin(lat_geo) * (altitude + R_c * (1 - e_squared));
- LengthType r = sqrt(p * p + z * z);
- double lat_sph = asin(z / r);
-
- const int length = 90;
- double epoch, g[length], h[length], g_dot[length], h_dot[length];
- std::string model_name;
- std::string release_date;
- int n[length], m[length];
-
- // Read in coefficients
- boost::filesystem::path const path = corsika::corsika_data("GeoMag/WMM.COF");
- boost::filesystem::ifstream file(path, std::ios::in);
- // Exit if file opening failed
- if (!file.is_open()) {
- CORSIKA_LOG_ERROR("Failed opening WMM.COF.");
- abort();
- }
-
- file >> epoch >> model_name >> release_date;
- for (int i = 0; i < length; i++) {
- file >> n[i] >> m[i] >> g[i] >> h[i] >> g_dot[i] >> h_dot[i];
- // Time interpolation
- g[i] = g[i] + (year - epoch) * g_dot[i];
- h[i] = h[i] + (year - epoch) * h_dot[i];
- }
- file.close();
-
- double legendre, next_legendre, derivate_legendre;
- double magneticfield[3] = {0, 0, 0};
-
- for (int j = 0; j < length; j++) {
- legendre = boost::math::legendre_p(n[j], m[j], sin(lat_sph));
- next_legendre = boost::math::legendre_p(n[j] + 1, m[j], sin(lat_sph));
-
- // Schmidt semi-normalization and Condon-Shortley phase term
- if (m[j] > 0) {
- legendre *= sqrt(2 * boost::math::factorial<double>(n[j] - m[j]) /
- boost::math::factorial<double>(n[j] + m[j])) *
- pow(-1, m[j]);
- next_legendre *= sqrt(2 * boost::math::factorial<double>(n[j] + 1 - m[j]) /
- boost::math::factorial<double>(n[j] + 1 + m[j])) *
- pow(-1, m[j]);
- }
- derivate_legendre =
- (n[j] + 1) * tan(lat_sph) * legendre -
- sqrt(pow(n[j] + 1, 2) - pow(m[j], 2)) / cos(lat_sph) * next_legendre;
-
- magneticfield[0] +=
- pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) *
- (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * derivate_legendre;
- magneticfield[1] +=
- pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) * m[j] *
- (g[j] * sin(m[j] * lon) - h[j] * cos(m[j] * lon)) * legendre;
- magneticfield[2] +=
- (n[j] + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) *
- (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * legendre;
- }
- magneticfield[0] *= -1;
- magneticfield[1] /= cos(lat_sph);
- magneticfield[2] *= -1;
-
- // Transform back into geodetic coordinates
- double magneticfield_geo[3];
- magneticfield_geo[0] = magneticfield[0] * cos(lat_sph - lat_geo) -
- magneticfield[2] * sin(lat_sph - lat_geo);
- magneticfield_geo[1] = magneticfield[1];
- magneticfield_geo[2] = magneticfield[0] * sin(lat_sph - lat_geo) +
- magneticfield[2] * cos(lat_sph - lat_geo);
-
- return MagneticFieldVector{Cs, magneticfield_geo[0] * 1_nT,
- magneticfield_geo[1] * 1_nT, magneticfield_geo[2] * -1_nT};
- }
-} // namespace corsika
\ No newline at end of file
diff --git a/corsika/detail/media/WorldMagneticModel.inl b/corsika/detail/media/WorldMagneticModel.inl
new file mode 100644
index 0000000000000000000000000000000000000000..ce1365320bc23ade677f4b4b1023c926bbf182c9
--- /dev/null
+++ b/corsika/detail/media/WorldMagneticModel.inl
@@ -0,0 +1,190 @@
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
+
+#include <boost/filesystem.hpp>
+#include <boost/math/special_functions/factorials.hpp>
+#include <boost/math/special_functions/legendre.hpp>
+
+#include <stdexcept>
+#include <string>
+#include <cmath>
+
+namespace corsika {
+
+ inline WorldMagneticModel::WorldMagneticModel(Point const& center,
+ std::string const& dataFile)
+ : 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);
+ throw std::runtime_error("Cannot load WorldMagneticModel data.");
+ }
+
+ // WorldMagneticModel supports two types of input data: WMM.COF and IGRF.COF
+ // They have only slightly different format and content and can be easily
+ // differentiated here.
+
+ std::string line;
+ while (getline(file >> std::ws, line)) {
+
+ double epoch;
+ std::string model_name;
+ std::string release_date; // just for WMM
+ int nMax = 12; // the spherical max n (l) shell (for IGRF), for WMM this is 12
+ // Note that n=l=0 is the monopole and is not included in the model.
+ int dummyInt;
+ double dummyDouble;
+
+ std::istringstream firstLine(line);
+
+ // check comments and ignore:
+ if (firstLine.peek() == '#' || // normal comment
+ line.size() == 0 || // empty line
+ line.find("9999999999999999999999999") == 0) { // crazy WMM comment
+ continue;
+ }
+
+ // check IGRF format:
+ if (firstLine >> model_name >> epoch >> nMax >> dummyInt >> dummyInt >>
+ dummyDouble >> dummyDouble >> dummyDouble >> dummyDouble >> model_name >>
+ dummyInt) {
+ static bool info = false;
+ if (!info) {
+ CORSIKA_LOG_INFO("Reading IGRF input data format.");
+ info = true;
+ }
+ } else {
+ // check WMM format:
+ firstLine.clear();
+ firstLine.seekg(0, std::ios::beg);
+ if (firstLine >> epoch >> model_name >> release_date) {
+ CORSIKA_LOG_INFO("Reading WMM input data format.");
+ } else {
+ CORSIKA_LOG_ERROR("line: {}", line);
+ throw std::runtime_error("Incompatible input data for WorldMagneticModel");
+ }
+ }
+
+ int nPar = 0;
+ for (int i = 0; i < nMax; ++i) { nPar += i + 2; }
+ int iEpoch = int(epoch);
+
+ if (parameters_.count(iEpoch) != 0) {
+ throw std::runtime_error(
+ "WorldMagneticModel input file has duplicate Epoch. Fix.");
+ }
+ parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
+
+ for (int i = 0; i < nPar; i++) {
+ file >> parameters_[iEpoch][i].n >> parameters_[iEpoch][i].m >>
+ parameters_[iEpoch][i].g >> parameters_[iEpoch][i].h >>
+ parameters_[iEpoch][i].dg >> parameters_[iEpoch][i].dh;
+ file.ignore(9999999, '\n');
+ }
+ }
+ file.close();
+ }
+
+ inline MagneticFieldVector WorldMagneticModel::getField(double const year,
+ LengthType const altitude,
+ double const latitude,
+ 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;
+ }
+ }
+ double epoch = double(iEpoch);
+ CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch, year);
+ if (iEpoch == 0) {
+ CORSIKA_LOG_WARN("Year {} is before first EPOCH. Results unclear.", year);
+ }
+ if (altitude < -1_km || altitude > 850_km) {
+ CORSIKA_LOG_WARN("Altitude should be between -1_km and 850_km.");
+ }
+ if (latitude <= -90 || latitude >= 90) {
+ CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
+ throw std::runtime_error("Latitude has to be between -90 and 90 degree.");
+ } else if (latitude < -89.992 || latitude > 89.992) {
+ CORSIKA_LOG_WARN("Latitude is close to the poles.");
+ }
+ if (longitude < -180 || longitude > 180) {
+ CORSIKA_LOG_WARN("Longitude should be between -180 and 180 degree.");
+ }
+
+ const double lat_geo = latitude * constants::pi / 180;
+ const double lon = longitude * constants::pi / 180;
+
+ // Transform into spherical coordinates
+ const double f = 1 / 298.257223563;
+ const double e_squared = f * (2 - f);
+ LengthType R_c =
+ constants::EarthRadius::Equatorial / sqrt(1 - e_squared * pow(sin(lat_geo), 2));
+ LengthType p = (R_c + altitude) * cos(lat_geo);
+ LengthType z = sin(lat_geo) * (altitude + R_c * (1 - e_squared));
+ LengthType r = sqrt(p * p + z * z);
+ double lat_sph = asin(z / r);
+
+ double legendre, next_legendre, derivate_legendre;
+ double magneticfield[3] = {0, 0, 0};
+
+ for (size_t j = 0; j < parameters_[iEpoch].size(); j++) {
+
+ ParameterLine p = parameters_[iEpoch][j];
+
+ // Time interpolation
+ p.g = p.g + (year - epoch) * p.dg;
+ p.h = p.h + (year - epoch) * p.dh;
+
+ legendre = boost::math::legendre_p(p.n, p.m, sin(lat_sph));
+ next_legendre = boost::math::legendre_p(p.n + 1, p.m, sin(lat_sph));
+
+ // Schmidt semi-normalization and Condon-Shortley phase term
+ if (p.m > 0) {
+ legendre *= sqrt(2 * boost::math::factorial<double>(p.n - p.m) /
+ boost::math::factorial<double>(p.n + p.m)) *
+ pow(-1, p.m);
+ next_legendre *= sqrt(2 * boost::math::factorial<double>(p.n + 1 - p.m) /
+ boost::math::factorial<double>(p.n + 1 + p.m)) *
+ pow(-1, p.m);
+ }
+ derivate_legendre =
+ (p.n + 1) * tan(lat_sph) * legendre -
+ sqrt(pow(p.n + 1, 2) - pow(p.m, 2)) / cos(lat_sph) * next_legendre;
+
+ magneticfield[0] +=
+ pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+ (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * derivate_legendre;
+ magneticfield[1] +=
+ pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) * p.m *
+ (p.g * sin(p.m * lon) - p.h * cos(p.m * lon)) * legendre;
+ magneticfield[2] +=
+ (p.n + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+ (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * legendre;
+ }
+ magneticfield[0] *= -1;
+ magneticfield[1] /= cos(lat_sph);
+ magneticfield[2] *= -1;
+
+ // Transform back into geodetic coordinates
+ double magneticfield_geo[3];
+ magneticfield_geo[0] = magneticfield[0] * cos(lat_sph - lat_geo) -
+ magneticfield[2] * sin(lat_sph - lat_geo);
+ magneticfield_geo[1] = magneticfield[1];
+ magneticfield_geo[2] = magneticfield[0] * sin(lat_sph - lat_geo) +
+ magneticfield[2] * cos(lat_sph - lat_geo);
+
+ return MagneticFieldVector{center_.getCoordinateSystem(), magneticfield_geo[0] * 1_nT,
+ magneticfield_geo[1] * 1_nT, magneticfield_geo[2] * -1_nT};
+ }
+
+} // namespace corsika
\ No newline at end of file
diff --git a/corsika/media/WMM.hpp b/corsika/media/WMM.hpp
deleted file mode 100644
index 9eb93b9d55854a88e059a3e0623a5509aadf0868..0000000000000000000000000000000000000000
--- a/corsika/media/WMM.hpp
+++ /dev/null
@@ -1,40 +0,0 @@
-/*
- * (c) Copyright 2020 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 <boost/filesystem.hpp>
-#include <boost/math/special_functions/factorials.hpp>
-#include <boost/math/special_functions/legendre.hpp>
-
-#include <cmath>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-#include <fstream>
-
-namespace corsika {
- /**
- * A magnetic field calculated with the WMM model
- *
- * @param year Year of the evaluation, between 2020 and 2025.
- * @param altitude Height of the location to evaluate the field at,
- in km between -1 and 850.
- * @param latitude Latitude of the location to evaluate the field at,
- in degrees between -90 and 90 (negative for southern hemisphere).
- * @param longitute Longitude of the location to evaluate the field at,
- in degrees between -180 and 180 (negative for western
- hemisphere).
- *
- * @returns The magnetic field vector in nT.
- *
- */
-
- inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year,
- const LengthType altitude, const double latitude,
- const double longitude);
-} // namespace corsika
-
-#include <corsika/detail/media/WMM.inl>
diff --git a/corsika/media/WorldMagneticModel.hpp b/corsika/media/WorldMagneticModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..642e25e73c34f64ddd7ac02e04f988957a1f4418
--- /dev/null
+++ b/corsika/media/WorldMagneticModel.hpp
@@ -0,0 +1,62 @@
+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+#include <fstream>
+#include <string>
+
+namespace corsika {
+
+ /**
+ * A magnetic field calculated with the WMM or IGRF model.
+ */
+
+ class WorldMagneticModel {
+
+ /**
+ * Internal data structure for a single shell of the spherical harmonic
+ * parameterization.
+ */
+ struct ParameterLine {
+ int n;
+ int m;
+ double g;
+ double h;
+ double dg; //< time dependence of g in "per year"
+ double dh; //< time dependence of h in "per year"
+ };
+
+ public:
+ /**
+ * Construct a new World Magnetic Model object.
+ *
+ * @param center Center of Earth.
+ * @param data Data table to read.
+ */
+ WorldMagneticModel(Point const& center, std::string const& data = "GeoMag/WMM.COF");
+
+ /**
+ * Calculates the value of the magnetic field.
+ *
+ * @param year Year of the evaluation, between 2020 and 2025.
+ * @param altitude Height of the location to evaluate the field at,
+ * in km between -1 and 850.
+ * @param latitude Latitude of the location to evaluate the field at,
+ * in degrees between -90 and 90 (negative for southern hemisphere).
+ * @param longitute Longitude of the location to evaluate the field at,
+ * in degrees between -180 and 180 (negative for western
+ * hemisphere).
+ *
+ * @returns The magnetic field vector in nT.
+ */
+ MagneticFieldVector getField(double const year, LengthType const altitude,
+ double const latitude, double const longitude);
+
+ private:
+ Point center_; //< Center of Earth.
+ std::map<int, std::vector<ParameterLine>> parameters_; //< epoch to Parameters map
+ };
+} // namespace corsika
+
+#include <corsika/detail/media/WorldMagneticModel.inl>
\ No newline at end of file
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index b792b7b16f4da697d096d6fae21361c8114e771a..27ee637bcde3bc27acabc9d9094a0217d51ae193 100644
--- a/tests/media/testMagneticField.cpp
+++ b/tests/media/testMagneticField.cpp
@@ -14,7 +14,7 @@
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/VolumeTreeNode.hpp>
-#include <corsika/media/WMM.hpp>
+#include <corsika/media/WorldMagneticModel.hpp>
#include <SetupTestTrajectory.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -25,7 +25,7 @@ using namespace corsika;
TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::trace);
CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
@@ -74,18 +74,42 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
}
SECTION("WorldMagneticModel") {
- // create earth magnetic field vector
- MagneticFieldVector Earth_B_1 = get_wmm(gCS, 2022.5, 100_km, -80, -120);
- MagneticFieldVector Earth_B_2 = get_wmm(gCS, 2022.5, 0_km, 0, 120);
- MagneticFieldVector Earth_B_3 = get_wmm(gCS, 2020, 100_km, 80, 0);
- CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5815).margin(0.03));
- CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(14803).margin(0.03));
- CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.03));
- CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(39684.7).margin(0.03));
- CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-42.2).margin(0.03));
- CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(10809.5).margin(0.03));
- CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6261.8).margin(0.03));
- CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(-185.5).margin(0.03));
- CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
+
+ CHECK_THROWS(WorldMagneticModel(gOrigin, "GeoMag/IDoNotExist.COF"));
+
+ {
+ WorldMagneticModel wmm(gOrigin, "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);
+ CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5815).margin(0.03));
+ CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(14803).margin(0.03));
+ CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.03));
+ CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(39684.7).margin(0.03));
+ CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-42.2).margin(0.03));
+ CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(10809.5).margin(0.03));
+ CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6261.8).margin(0.03));
+ CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(-185.5).margin(0.03));
+ CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
+ }
+ {
+ WorldMagneticModel wmm(gOrigin, "GeoMag/IGRF13.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);
+ CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5815).margin(0.03));
+ CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(14803).margin(0.03));
+ CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.03));
+ CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(39684.7).margin(0.03));
+ CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-42.2).margin(0.03));
+ CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(10809.5).margin(0.03));
+ CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6261.8).margin(0.03));
+ CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(-185.5).margin(0.03));
+ CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
+ }
}
}