diff --git a/corsika/detail/media/GeomagneticModel.inl b/corsika/detail/media/GeomagneticModel.inl index ed92768381b953204a48bd5197ecd8af99f2cdd0..1747f95132082c26a06f58ccd2ef5d9e482d1395 100644 --- a/corsika/detail/media/GeomagneticModel.inl +++ b/corsika/detail/media/GeomagneticModel.inl @@ -95,11 +95,16 @@ namespace corsika { int iYear = int(year); int iEpoch = 0; + int next_Epoch = 0; for (auto parIt = parameters_.rbegin(); parIt != parameters_.rend(); ++parIt) { if (parIt->first <= iYear) { iEpoch = parIt->first; break; } + if (parIt->first >= iYear) { + next_Epoch = parIt->first; + break; + } } double epoch = double(iEpoch); CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch, year); @@ -142,6 +147,11 @@ namespace corsika { // Time interpolation p.g = p.g + (year - epoch) * p.dg; p.h = p.h + (year - epoch) * p.dh; + if (next_Epoch != 0) { + ParameterLine next_p = parameters_[next_Epoch][j]; + p.g = p.g + (next_p.g - p.g) * (year - epoch) / (next_Epoch - epoch); + p.h = p.h + (next_p.h - p.h) * (year - epoch) / (next_Epoch - epoch); + } legendre = pow(-1, p.m) * std::assoc_legendre(p.n, p.m, sin(lat_sph)); next_legendre = pow(-1, p.m) * std::assoc_legendre(p.n + 1, p.m, sin(lat_sph)); diff --git a/tests/media/testMagneticField.cpp b/tests/media/testMagneticField.cpp index d6b0c9adb05e6d431f2354255da448658a2e12de..de8b12fabb1522e8031c86cc0661131b9fae4e94 100644 --- a/tests/media/testMagneticField.cpp +++ b/tests/media/testMagneticField.cpp @@ -99,21 +99,21 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") { // create earth magnetic field vector MagneticFieldVector Earth_B_1 = igrf.getField(2022.5, 100_km, -80, -120); - MagneticFieldVector Earth_B_2 = igrf.getField(1935, 10_km, 80, 120); + MagneticFieldVector Earth_B_2 = igrf.getField(1937.5, 10_km, 80, 120); MagneticFieldVector Earth_B_3 = igrf.getField(1990, 50_km, 80, 0); - MagneticFieldVector Earth_B_4 = igrf.getField(2010, 0_km, 0, -120); + MagneticFieldVector Earth_B_4 = igrf.getField(2012, 0_km, 0, -120); CHECK(Earth_B_1.getX(gCS) / 1_nT == Approx(5821).margin(0.5)); CHECK(Earth_B_1.getY(gCS) / 1_nT == Approx(-14796).margin(0.5)); - CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(49776).margin(0.5)); - CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(3014).margin(0.5)); - CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-1072).margin(0.5)); - CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(-57594).margin(0.5)); + CHECK(Earth_B_1.getZ(gCS) / 1_nT == Approx(497756).margin(0.5)); + CHECK(Earth_B_2.getX(gCS) / 1_nT == Approx(2957).margin(0.5)); + CHECK(Earth_B_2.getY(gCS) / 1_nT == Approx(-1014).margin(0.5)); + CHECK(Earth_B_2.getZ(gCS) / 1_nT == Approx(-57610).margin(0.5)); CHECK(Earth_B_3.getX(gCS) / 1_nT == Approx(6545).margin(0.5)); CHECK(Earth_B_3.getY(gCS) / 1_nT == Approx(1473).margin(0.5)); CHECK(Earth_B_3.getZ(gCS) / 1_nT == Approx(-52884).margin(0.5)); - CHECK(Earth_B_4.getX(gCS) / 1_nT == Approx(30137).margin(0.5)); - CHECK(Earth_B_4.getY(gCS) / 1_nT == Approx(-4805).margin(0.5)); - CHECK(Earth_B_4.getZ(gCS) / 1_nT == Approx(-5161).margin(0.5)); + CHECK(Earth_B_4.getX(gCS) / 1_nT == Approx(30052).margin(0.5)); + CHECK(Earth_B_4.getY(gCS) / 1_nT == Approx(-4738).margin(0.5)); + CHECK(Earth_B_4.getZ(gCS) / 1_nT == Approx(-5190).margin(0.5)); } } }