diff --git a/corsika/detail/media/WorldMagneticModel.inl b/corsika/detail/media/GeomagneticModel.inl
similarity index 93%
rename from corsika/detail/media/WorldMagneticModel.inl
rename to corsika/detail/media/GeomagneticModel.inl
index 225ef0d5cceb5e82f19043da312cb27b9ef0b17f..6bf0ea20edd43b45c52c7fd3069abbb87491a4af 100644
--- a/corsika/detail/media/WorldMagneticModel.inl
+++ b/corsika/detail/media/GeomagneticModel.inl
@@ -9,7 +9,7 @@
namespace corsika {
- inline WorldMagneticModel::WorldMagneticModel(Point const& center,
+ inline GeomagneticModel::GeomagneticModel(Point const& center,
std::string const& dataFile)
: center_(center) {
@@ -20,10 +20,10 @@ namespace corsika {
// 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.");
+ throw std::runtime_error("Cannot load GeomagneticModel data.");
}
- // WorldMagneticModel supports two types of input data: WMM.COF and IGRF.COF
+ // GeomagneticModel 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.
@@ -64,7 +64,7 @@ namespace corsika {
CORSIKA_LOG_INFO("Reading WMM input data format.");
} else {
CORSIKA_LOG_ERROR("line: {}", line);
- throw std::runtime_error("Incompatible input data for WorldMagneticModel");
+ throw std::runtime_error("Incompatible input data for GeomagneticModel");
}
}
@@ -74,7 +74,7 @@ namespace corsika {
if (parameters_.count(iEpoch) != 0) {
throw std::runtime_error(
- "WorldMagneticModel input file has duplicate Epoch. Fix.");
+ "GeomagneticModel input file has duplicate Epoch. Fix.");
}
parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
@@ -88,7 +88,7 @@ namespace corsika {
file.close();
}
- inline MagneticFieldVector WorldMagneticModel::getField(double const year,
+ inline MagneticFieldVector GeomagneticModel::getField(double const year,
LengthType const altitude,
double const latitude,
double const longitude) {
@@ -182,7 +182,7 @@ namespace corsika {
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};
+ magneticfield_geo[1] * -1_nT, magneticfield_geo[2] * -1_nT};
}
} // namespace corsika
\ No newline at end of file
diff --git a/corsika/media/WorldMagneticModel.hpp b/corsika/media/GeomagneticModel.hpp
similarity index 86%
rename from corsika/media/WorldMagneticModel.hpp
rename to corsika/media/GeomagneticModel.hpp
index 642e25e73c34f64ddd7ac02e04f988957a1f4418..27a8f59c84d84b29fddb8a31b6f580ffe95404a4 100644
--- a/corsika/media/WorldMagneticModel.hpp
+++ b/corsika/media/GeomagneticModel.hpp
@@ -10,9 +10,11 @@ namespace corsika {
/**
* A magnetic field calculated with the WMM or IGRF model.
+ * WMM: https://geomag.bgs.ac.uk/documents/WMM2020_Report.pdf
+ * IGRF: https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html
*/
- class WorldMagneticModel {
+ class GeomagneticModel {
/**
* Internal data structure for a single shell of the spherical harmonic
@@ -34,7 +36,7 @@ namespace corsika {
* @param center Center of Earth.
* @param data Data table to read.
*/
- WorldMagneticModel(Point const& center, std::string const& data = "GeoMag/WMM.COF");
+ GeomagneticModel(Point const& center, std::string const& data = "GeoMag/WMM.COF");
/**
* Calculates the value of the magnetic field.
@@ -59,4 +61,4 @@ namespace corsika {
};
} // namespace corsika
-#include <corsika/detail/media/WorldMagneticModel.inl>
\ No newline at end of file
+#include <corsika/detail/media/GeomagneticModel.inl>
\ No newline at end of file
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index cb30870f4d3c6fbdc69ff5c612f77702c12e1d4c..3bb04e08c5462eb4967e3356aa39945835c68957 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -33,6 +33,7 @@
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
@@ -41,7 +42,6 @@
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
-#include <corsika/media/WMM.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
@@ -200,11 +200,12 @@ int main(int argc, char** argv) {
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
+ GeomagneticModel wmm(center, "GeoMag/WMM.COF");
// build a Linsley US Standard atmosphere into `env`
create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
- get_wmm(rootCS, 2022.5, 10_km, 49, 8.4));
+ wmm.getField(2022.5, 10_km, 49, 8.4));
/* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 6e12bb9c404ea434ffb071d30cc134faa8b70b79..dba4f18fe09c9869fee217af830f7bb1e2e51342 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -35,6 +35,7 @@
#include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp>
+#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/NuclearComposition.hpp>
@@ -42,7 +43,6 @@
#include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp>
-#include <corsika/media/WMM.hpp>
#include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp>
@@ -126,11 +126,12 @@ int main(int argc, char** argv) {
EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m};
+ GeomagneticModel wmm(center, "GeoMag/WMM.COF");
// build a Linsley US Standard atmosphere into `env`
create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm,
- get_wmm(rootCS, 2022.5, 10_km, 49, 8.4));
+ wmm.getField(2022.5, 10_km, 49, 8.4));
// pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1]));
diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp
index 27ee637bcde3bc27acabc9d9094a0217d51ae193..1d9b8e6751ed0746e4acf7b057e3e7dd3e880642 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/WorldMagneticModel.hpp>
+#include <corsika/media/GeomagneticModel.hpp>
#include <SetupTestTrajectory.hpp>
#include <corsika/setup/SetupTrajectory.hpp>
@@ -73,42 +73,42 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS));
}
- SECTION("WorldMagneticModel") {
+ SECTION("GeomagneticModel") {
- CHECK_THROWS(WorldMagneticModel(gOrigin, "GeoMag/IDoNotExist.COF"));
+ CHECK_THROWS(GeomagneticModel(gOrigin, "GeoMag/IDoNotExist.COF"));
{
- WorldMagneticModel wmm(gOrigin, "GeoMag/WMM.COF");
+ GeomagneticModel 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));
+ 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));
+ CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(39684.7).margin(0.05));
+ CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(42.2).margin(0.05));
+ CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(10809.5).margin(0.05));
+ 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));
}
{
- WorldMagneticModel wmm(gOrigin, "GeoMag/IGRF13.COF");
+ GeomagneticModel igrf(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);
+ MagneticFieldVector Earth_B_1 = igrf.getField(2022.5, 100_km, -80, -120);
+ MagneticFieldVector Earth_B_2 = igrf.getField(2022.5, 0_km, 0, 120);
+ MagneticFieldVector Earth_B_3 = igrf.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.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.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.getY(gCS) / 1_nT == Approx(185.5).margin(0.03));
CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
}
}