IAP GITLAB

Skip to content
Snippets Groups Projects
Commit 1c42ced0 authored by André Schmidt's avatar André Schmidt Committed by ralfulrich
Browse files

named GeomagneticModel

parent b7d0a10a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
corsika/detail/media/WorldMagneticModel.inl → corsika/detail/media/GeomagneticModel.inl
+ 7
− 7
View file @ 1c42ced0
...@@ -9,7 +9,7 @@ ...@@ -9,7 +9,7 @@
namespace corsika { namespace corsika {
inline WorldMagneticModel::WorldMagneticModel(Point const& center, inline GeomagneticModel::GeomagneticModel(Point const& center,
std::string const& dataFile) std::string const& dataFile)
: center_(center) { : center_(center) {
...@@ -20,10 +20,10 @@ namespace corsika { ...@@ -20,10 +20,10 @@ namespace corsika {
// Exit if file opening failed // Exit if file opening failed
if (!file.is_open()) { if (!file.is_open()) {
CORSIKA_LOG_ERROR("Failed opening data file {}", dataFile); 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 // They have only slightly different format and content and can be easily
// differentiated here. // differentiated here.
...@@ -64,7 +64,7 @@ namespace corsika { ...@@ -64,7 +64,7 @@ namespace corsika {
CORSIKA_LOG_INFO("Reading WMM input data format."); CORSIKA_LOG_INFO("Reading WMM input data format.");
} else { } else {
CORSIKA_LOG_ERROR("line: {}", line); 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 { ...@@ -74,7 +74,7 @@ namespace corsika {
if (parameters_.count(iEpoch) != 0) { if (parameters_.count(iEpoch) != 0) {
throw std::runtime_error( throw std::runtime_error(
"WorldMagneticModel input file has duplicate Epoch. Fix."); "GeomagneticModel input file has duplicate Epoch. Fix.");
} }
parameters_[iEpoch] = std::vector<ParameterLine>(nPar); parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
...@@ -88,7 +88,7 @@ namespace corsika { ...@@ -88,7 +88,7 @@ namespace corsika {
file.close(); file.close();
} }
inline MagneticFieldVector WorldMagneticModel::getField(double const year, inline MagneticFieldVector GeomagneticModel::getField(double const year,
LengthType const altitude, LengthType const altitude,
double const latitude, double const latitude,
double const longitude) { double const longitude) {
...@@ -182,7 +182,7 @@ namespace corsika { ...@@ -182,7 +182,7 @@ namespace corsika {
magneticfield[2] * cos(lat_sph - lat_geo); magneticfield[2] * cos(lat_sph - lat_geo);
return MagneticFieldVector{center_.getCoordinateSystem(), magneticfield_geo[0] * 1_nT, 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 } // namespace corsika
\ No newline at end of file
corsika/media/WorldMagneticModel.hpp → corsika/media/GeomagneticModel.hpp
+ 5
− 3
View file @ 1c42ced0
...@@ -10,9 +10,11 @@ namespace corsika { ...@@ -10,9 +10,11 @@ namespace corsika {
/** /**
* A magnetic field calculated with the WMM or IGRF model. * 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 * Internal data structure for a single shell of the spherical harmonic
...@@ -34,7 +36,7 @@ namespace corsika { ...@@ -34,7 +36,7 @@ namespace corsika {
* @param center Center of Earth. * @param center Center of Earth.
* @param data Data table to read. * @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. * Calculates the value of the magnetic field.
...@@ -59,4 +61,4 @@ namespace corsika { ...@@ -59,4 +61,4 @@ namespace corsika {
}; };
} // namespace corsika } // namespace corsika
#include <corsika/detail/media/WorldMagneticModel.inl> #include <corsika/detail/media/GeomagneticModel.inl>
\ No newline at end of file \ No newline at end of file
examples/corsika.cpp
+ 3
− 2
View file @ 1c42ced0
...@@ -33,6 +33,7 @@ ...@@ -33,6 +33,7 @@
#include <corsika/media/Environment.hpp> #include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp> #include <corsika/media/FlatExponential.hpp>
#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp> #include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp> #include <corsika/media/LayeredSphericalAtmosphereBuilder.hpp>
...@@ -41,7 +42,6 @@ ...@@ -41,7 +42,6 @@
#include <corsika/media/UniformMagneticField.hpp> #include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp> #include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp> #include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/WMM.hpp>
#include <corsika/modules/BetheBlochPDG.hpp> #include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp> #include <corsika/modules/LongitudinalProfile.hpp>
...@@ -200,11 +200,12 @@ int main(int argc, char** argv) { ...@@ -200,11 +200,12 @@ int main(int argc, char** argv) {
EnvType env; EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem(); CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m}; Point const center{rootCS, 0_m, 0_m, 0_m};
GeomagneticModel wmm(center, "GeoMag/WMM.COF");
// build a Linsley US Standard atmosphere into `env` // build a Linsley US Standard atmosphere into `env`
create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>( create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, 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 === */ /* === END: SETUP ENVIRONMENT AND ROOT COORDINATE SYSTEM === */
......
examples/vertical_EAS.cpp
+ 3
− 2
View file @ 1c42ced0
...@@ -35,6 +35,7 @@ ...@@ -35,6 +35,7 @@
#include <corsika/media/Environment.hpp> #include <corsika/media/Environment.hpp>
#include <corsika/media/FlatExponential.hpp> #include <corsika/media/FlatExponential.hpp>
#include <corsika/media/GeomagneticModel.hpp>
#include <corsika/media/HomogeneousMedium.hpp> #include <corsika/media/HomogeneousMedium.hpp>
#include <corsika/media/IMagneticFieldModel.hpp> #include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/NuclearComposition.hpp> #include <corsika/media/NuclearComposition.hpp>
...@@ -42,7 +43,6 @@ ...@@ -42,7 +43,6 @@
#include <corsika/media/UniformMagneticField.hpp> #include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/ShowerAxis.hpp> #include <corsika/media/ShowerAxis.hpp>
#include <corsika/media/CORSIKA7Atmospheres.hpp> #include <corsika/media/CORSIKA7Atmospheres.hpp>
#include <corsika/media/WMM.hpp>
#include <corsika/modules/BetheBlochPDG.hpp> #include <corsika/modules/BetheBlochPDG.hpp>
#include <corsika/modules/LongitudinalProfile.hpp> #include <corsika/modules/LongitudinalProfile.hpp>
...@@ -126,11 +126,12 @@ int main(int argc, char** argv) { ...@@ -126,11 +126,12 @@ int main(int argc, char** argv) {
EnvType env; EnvType env;
CoordinateSystemPtr const& rootCS = env.getCoordinateSystem(); CoordinateSystemPtr const& rootCS = env.getCoordinateSystem();
Point const center{rootCS, 0_m, 0_m, 0_m}; Point const center{rootCS, 0_m, 0_m, 0_m};
GeomagneticModel wmm(center, "GeoMag/WMM.COF");
// build a Linsley US Standard atmosphere into `env` // build a Linsley US Standard atmosphere into `env`
create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>( create_5layer_atmosphere<setup::EnvironmentInterface, MyExtraEnv>(
env, AtmosphereId::LinsleyUSStd, center, Medium::AirDry1Atm, 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 // pre-setup particle stack
unsigned short const A = std::stoi(std::string(argv[1])); unsigned short const A = std::stoi(std::string(argv[1]));
......
tests/media/testMagneticField.cpp
+ 20
− 20
View file @ 1c42ced0
...@@ -14,7 +14,7 @@ ...@@ -14,7 +14,7 @@
#include <corsika/media/UniformMagneticField.hpp> #include <corsika/media/UniformMagneticField.hpp>
#include <corsika/media/IMagneticFieldModel.hpp> #include <corsika/media/IMagneticFieldModel.hpp>
#include <corsika/media/VolumeTreeNode.hpp> #include <corsika/media/VolumeTreeNode.hpp>
#include <corsika/media/WorldMagneticModel.hpp> #include <corsika/media/GeomagneticModel.hpp>
#include <SetupTestTrajectory.hpp> #include <SetupTestTrajectory.hpp>
#include <corsika/setup/SetupTrajectory.hpp> #include <corsika/setup/SetupTrajectory.hpp>
...@@ -73,42 +73,42 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { ...@@ -73,42 +73,42 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
medium.getMagneticField(Point(gCS, 0_m, 0_m, 0_m)).getComponents(gCS)); 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 // create earth magnetic field vector
MagneticFieldVector Earth_B_1 = wmm.getField(2022.5, 100_km, -80, -120); 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_2 = wmm.getField(2022.5, 0_km, 0, 120);
MagneticFieldVector Earth_B_3 = wmm.getField(2020, 100_km, 80, 0); 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.getX(gCS) / 1_nT == Approx(5815).margin(0.05));
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.05));
CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49755.3).margin(0.03)); 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.03)); 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.03)); 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.03)); 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.03)); 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.03)); 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.03)); 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 // create earth magnetic field vector
MagneticFieldVector Earth_B_1 = wmm.getField(2022.5, 100_km, -80, -120); MagneticFieldVector Earth_B_1 = igrf.getField(2022.5, 100_km, -80, -120);
MagneticFieldVector Earth_B_2 = wmm.getField(2022.5, 0_km, 0, 120); MagneticFieldVector Earth_B_2 = igrf.getField(2022.5, 0_km, 0, 120);
MagneticFieldVector Earth_B_3 = wmm.getField(2020, 100_km, 80, 0); 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.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_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.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_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.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)); CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52429.1).margin(0.03));
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

data privacy - impressum