IAP GITLAB

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

named GeomagneticModel

parent 7021a3f9
No related branches found
No related tags found
1 merge request!404adding Geomagnetic Models
Hide whitespace changes
Inline Side-by-side
corsika/detail/media/WorldMagneticModel.inl → corsika/detail/media/GeomagneticModel.inl
+ 7
− 7
View file @ e441a86c
......@@ -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
corsika/media/WorldMagneticModel.hpp → corsika/media/GeomagneticModel.hpp
+ 5
− 3
View file @ e441a86c
......@@ -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
examples/corsika.cpp
+ 3
− 2
View file @ e441a86c
......@@ -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 === */
......
examples/vertical_EAS.cpp
+ 3
− 2
View file @ e441a86c
......@@ -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]));
......
tests/media/testMagneticField.cpp
+ 20
− 20
View file @ e441a86c
......@@ -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));
}
}
......
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