diff --git a/corsika/detail/media/WMM.inl b/corsika/detail/media/WMM.inl
new file mode 100644
index 0000000000000000000000000000000000000000..3eecac4800bd830af5fceadfd8ce43d476d384d9
--- /dev/null
+++ b/corsika/detail/media/WMM.inl
@@ -0,0 +1,121 @@
+/*
+ * (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/media/WMM.hpp b/corsika/media/WMM.hpp
index a1f84351a4b51e679691c92a14aca1814a9d1798..9eb93b9d55854a88e059a3e0623a5509aadf0868 100644
--- a/corsika/media/WMM.hpp
+++ b/corsika/media/WMM.hpp
@@ -25,14 +25,16 @@ namespace corsika {
* @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).
+ 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);
+
+ 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/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index 1bbe1b718a2c9ed10dc6bc256ad00a53ab70f5c9..b792b7b16f4da697d096d6fae21361c8114e771a 100644
--- a/tests/media/testMagneticField.cpp
+++ b/tests/media/testMagneticField.cpp
@@ -30,72 +30,62 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
Point const gOrigin(gCS, {0_m, 0_m, 0_m});
- // setup our interface types
- using IModelInterface = IMagneticFieldModel<IMediumModel>;
- using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;
-
- // the composition we use for the homogenous medium
- NuclearComposition const protonComposition({Code::Proton}, {1.});
-
- // create a magnetic field vector
- Vector B0(gCS, 0_T, 0_T, 0_T);
-
- // the constant density
- const auto density{19.2_g / cube(1_cm)};
-
- // create our atmospheric model
- AtmModel medium(B0, density, protonComposition);
-
- // and test at several locations
- CHECK(B0.getComponents(gCS) ==
- medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
- CHECK(
- B0.getComponents(gCS) ==
- medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).getComponents(gCS));
- CHECK(B0.getComponents(gCS) ==
- medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
-
- // create a new magnetic field vector
- Vector B1(gCS, 23_T, 57_T, -4_T);
-
- // and update this atmospheric model
- medium.setMagneticField(B1);
-
- // and test at several locations
- CHECK(B1.getComponents(gCS) ==
- medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
- CHECK(
- B1.getComponents(gCS) ==
- medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km)).getComponents(gCS));
- CHECK(B1.getComponents(gCS) ==
- medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
-
- // check the density and nuclear composition
- CHECK(density == medium.getMassDensity(Point(gCS, 0_m, 0_m, 0_m)));
- medium.getNuclearComposition();
-
- // create a line of length 1 m
- Line const line(gOrigin, Vector<SpeedType::dimension_type>(
- gCS, {1_m / second, 0_m / second, 0_m / second}));
-
- // the end time of our line
- auto const tEnd = 1_s;
-
- // and the associated trajectory
- setup::Trajectory const track =
- setup::testing::make_track<setup::Trajectory>(line, tEnd);
-
- // 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));
+ SECTION("UniformMagneticField interface") {
+
+ // setup our interface types
+ using IModelInterface = IMagneticFieldModel<IMediumModel>;
+ using AtmModel = UniformMagneticField<HomogeneousMedium<IModelInterface>>;
+
+ // the composition we use for the homogenous medium
+ NuclearComposition const protonComposition({Code::Proton}, {1.});
+
+ // create a magnetic field vector
+ Vector B0(gCS, 0_T, 0_T, 0_T);
+
+ // the constant density
+ const auto density{19.2_g / cube(1_cm)};
+
+ // create our atmospheric model
+ AtmModel medium(B0, density, protonComposition);
+
+ // and test at several locations
+ CHECK(B0.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
+ CHECK(B0.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))
+ .getComponents(gCS));
+ CHECK(B0.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
+
+ // create a new magnetic field vector
+ Vector B1(gCS, 23_T, 57_T, -4_T);
+
+ // and update this atmospheric model
+ medium.setMagneticField(B1);
+
+ // and test at several locations
+ CHECK(B1.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, -10_m, 4_m, 35_km)).getComponents(gCS));
+ CHECK(B1.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 1000_km, -1000_km, 1000_km))
+ .getComponents(gCS));
+ CHECK(B1.getComponents(gCS) ==
+ medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
+ }
+
+ 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));
+ }
}