interpolation between epochs

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));

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));
     }
   }
 }