new class

corsika/detail/media/WMM.inl → corsika/detail/media/WorldMagneticModel.inl

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#pragma once
+#include <corsika/framework/core/Logging.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
 
+#include <boost/filesystem.hpp>
 #include <boost/math/special_functions/factorials.hpp>
 #include <boost/math/special_functions/legendre.hpp>
 
+#include <stdexcept>
+#include <string>
 #include <cmath>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-#include <fstream>
 
 namespace corsika {
-  inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year,
-                                     const LengthType altitude, const double latitude,
-                                     const double longitude) {
-    if (year < 2020 || year > 2025) {
-      CORSIKA_LOG_WARN("Year has to be between 2020 and 2025.");
+
+  inline WorldMagneticModel::WorldMagneticModel(Point const& center,
+                                                std::string const& dataFile)
+      : center_(center) {
+
+    // Read in coefficients
+    boost::filesystem::path const path = corsika::corsika_data(dataFile);
+    boost::filesystem::ifstream file(path, std::ios::in);
+
+    // 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.");
+    }
+
+    // WorldMagneticModel 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.
+
+    std::string line;
+    while (getline(file >> std::ws, line)) {
+
+      double epoch;
+      std::string model_name;
+      std::string release_date; // just for WMM
+      int nMax = 12; // the spherical max n (l) shell (for IGRF), for WMM this is 12
+                     // Note that n=l=0 is the monopole and is not included in the model.
+      int dummyInt;
+      double dummyDouble;
+
+      std::istringstream firstLine(line);
+
+      // check comments and ignore:
+      if (firstLine.peek() == '#' ||                     // normal comment
+          line.size() == 0 ||                            // empty line
+          line.find("9999999999999999999999999") == 0) { // crazy WMM comment
+        continue;
+      }
+
+      // check IGRF format:
+      if (firstLine >> model_name >> epoch >> nMax >> dummyInt >> dummyInt >>
+          dummyDouble >> dummyDouble >> dummyDouble >> dummyDouble >> model_name >>
+          dummyInt) {
+        static bool info = false;
+        if (!info) {
+          CORSIKA_LOG_INFO("Reading IGRF input data format.");
+          info = true;
+        }
+      } else {
+        // check WMM format:
+        firstLine.clear();
+        firstLine.seekg(0, std::ios::beg);
+        if (firstLine >> epoch >> model_name >> release_date) {
+          CORSIKA_LOG_INFO("Reading WMM input data format.");
+        } else {
+          CORSIKA_LOG_ERROR("line: {}", line);
+          throw std::runtime_error("Incompatible input data for WorldMagneticModel");
+        }
+      }
+
+      int nPar = 0;
+      for (int i = 0; i < nMax; ++i) { nPar += i + 2; }
+      int iEpoch = int(epoch);
+
+      if (parameters_.count(iEpoch) != 0) {
+        throw std::runtime_error(
+            "WorldMagneticModel input file has duplicate Epoch. Fix.");
+      }
+      parameters_[iEpoch] = std::vector<ParameterLine>(nPar);
+
+      for (int i = 0; i < nPar; i++) {
+        file >> parameters_[iEpoch][i].n >> parameters_[iEpoch][i].m >>
+            parameters_[iEpoch][i].g >> parameters_[iEpoch][i].h >>
+            parameters_[iEpoch][i].dg >> parameters_[iEpoch][i].dh;
+        file.ignore(9999999, '\n');
+      }
+    }
+    file.close();
+  }
+
+  inline MagneticFieldVector WorldMagneticModel::getField(double const year,
+                                                          LengthType const altitude,
+                                                          double const latitude,
+                                                          double const longitude) {
+
+    int iYear = int(year);
+    int iEpoch = 0;
+    for (auto parIt = parameters_.rbegin(); parIt != parameters_.rend(); ++parIt) {
+      if (parIt->first <= iYear) {
+        iEpoch = parIt->first;
+        break;
+      }
+    }
+    double epoch = double(iEpoch);
+    CORSIKA_LOG_DEBUG("Found Epoch {} for year {}", iEpoch, year);
+    if (iEpoch == 0) {
+      CORSIKA_LOG_WARN("Year {} is before first EPOCH. Results unclear.", year);
     }
     if (altitude < -1_km || altitude > 850_km) {
       CORSIKA_LOG_WARN("Altitude should be between -1_km and 850_km.");
     }
-    if (latitude < -90 || latitude > 90) {
+    if (latitude <= -90 || latitude >= 90) {
       CORSIKA_LOG_ERROR("Latitude has to be between -90 and 90 degree.");
-      abort();
+      throw std::runtime_error("Latitude has to be between -90 and 90 degree.");
     } else if (latitude < -89.992 || latitude > 89.992) {
       CORSIKA_LOG_WARN("Latitude is close to the poles.");
     }
@@ -49,59 +134,42 @@ namespace corsika {
     LengthType r = sqrt(p * p + z * z);
     double lat_sph = asin(z / r);
 
-    const int length = 90;
-    double epoch, g[length], h[length], g_dot[length], h_dot[length];
-    std::string model_name;
-    std::string release_date;
-    int n[length], m[length];
+    double legendre, next_legendre, derivate_legendre;
+    double magneticfield[3] = {0, 0, 0};
 
-    // Read in coefficients
-    boost::filesystem::path const path = corsika::corsika_data("GeoMag/WMM.COF");
-    boost::filesystem::ifstream file(path, std::ios::in);
-    // Exit if file opening failed
-    if (!file.is_open()) {
-      CORSIKA_LOG_ERROR("Failed opening WMM.COF.");
-      abort();
-    }
+    for (size_t j = 0; j < parameters_[iEpoch].size(); j++) {
 
-    file >> epoch >> model_name >> release_date;
-    for (int i = 0; i < length; i++) {
-      file >> n[i] >> m[i] >> g[i] >> h[i] >> g_dot[i] >> h_dot[i];
-      // Time interpolation
-      g[i] = g[i] + (year - epoch) * g_dot[i];
-      h[i] = h[i] + (year - epoch) * h_dot[i];
-    }
-    file.close();
+      ParameterLine p = parameters_[iEpoch][j];
 
-    double legendre, next_legendre, derivate_legendre;
-    double magneticfield[3] = {0, 0, 0};
+      // Time interpolation
+      p.g = p.g + (year - epoch) * p.dg;
+      p.h = p.h + (year - epoch) * p.dh;
 
-    for (int j = 0; j < length; j++) {
-      legendre = boost::math::legendre_p(n[j], m[j], sin(lat_sph));
-      next_legendre = boost::math::legendre_p(n[j] + 1, m[j], sin(lat_sph));
+      legendre = boost::math::legendre_p(p.n, p.m, sin(lat_sph));
+      next_legendre = boost::math::legendre_p(p.n + 1, p.m, sin(lat_sph));
 
       // Schmidt semi-normalization and Condon-Shortley phase term
-      if (m[j] > 0) {
-        legendre *= sqrt(2 * boost::math::factorial<double>(n[j] - m[j]) /
-                         boost::math::factorial<double>(n[j] + m[j])) *
-                    pow(-1, m[j]);
-        next_legendre *= sqrt(2 * boost::math::factorial<double>(n[j] + 1 - m[j]) /
-                              boost::math::factorial<double>(n[j] + 1 + m[j])) *
-                         pow(-1, m[j]);
+      if (p.m > 0) {
+        legendre *= sqrt(2 * boost::math::factorial<double>(p.n - p.m) /
+                         boost::math::factorial<double>(p.n + p.m)) *
+                    pow(-1, p.m);
+        next_legendre *= sqrt(2 * boost::math::factorial<double>(p.n + 1 - p.m) /
+                              boost::math::factorial<double>(p.n + 1 + p.m)) *
+                         pow(-1, p.m);
       }
       derivate_legendre =
-          (n[j] + 1) * tan(lat_sph) * legendre -
-          sqrt(pow(n[j] + 1, 2) - pow(m[j], 2)) / cos(lat_sph) * next_legendre;
+          (p.n + 1) * tan(lat_sph) * legendre -
+          sqrt(pow(p.n + 1, 2) - pow(p.m, 2)) / cos(lat_sph) * next_legendre;
 
       magneticfield[0] +=
-          pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) *
-          (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * derivate_legendre;
+          pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+          (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * derivate_legendre;
       magneticfield[1] +=
-          pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) * m[j] *
-          (g[j] * sin(m[j] * lon) - h[j] * cos(m[j] * lon)) * legendre;
+          pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) * p.m *
+          (p.g * sin(p.m * lon) - p.h * cos(p.m * lon)) * legendre;
       magneticfield[2] +=
-          (n[j] + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, n[j] + 2) *
-          (g[j] * cos(m[j] * lon) + h[j] * sin(m[j] * lon)) * legendre;
+          (p.n + 1) * pow(constants::EarthRadius::Geomagnetic_reference / r, p.n + 2) *
+          (p.g * cos(p.m * lon) + p.h * sin(p.m * lon)) * legendre;
     }
     magneticfield[0] *= -1;
     magneticfield[1] /= cos(lat_sph);
@@ -115,7 +183,8 @@ namespace corsika {
     magneticfield_geo[2] = magneticfield[0] * sin(lat_sph - lat_geo) +
                            magneticfield[2] * cos(lat_sph - lat_geo);
 
-    return MagneticFieldVector{Cs, magneticfield_geo[0] * 1_nT,
+    return MagneticFieldVector{center_.getCoordinateSystem(), magneticfield_geo[0] * 1_nT,
                                magneticfield_geo[1] * 1_nT, magneticfield_geo[2] * -1_nT};
   }
+
 } // namespace corsika
\ No newline at end of file

corsika/media/WMM.hpp

-/*
- * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
- *
- * This software is distributed under the terms of the GNU General Public
- * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
- * the license.
- */
-
-#include <boost/filesystem.hpp>
-#include <boost/math/special_functions/factorials.hpp>
-#include <boost/math/special_functions/legendre.hpp>
-
-#include <cmath>
-#include <corsika/framework/core/Logging.hpp>
-#include <corsika/framework/utility/CorsikaData.hpp>
-#include <fstream>
-
-namespace corsika {
-  /**
-   * A magnetic field calculated with the WMM model
-   *
-   * @param  year        Year of the evaluation, between 2020 and 2025.
-   * @param  altitude    Height of the location to evaluate the field at,
-                         in km between -1 and 850.
-   * @param  latitude    Latitude of the location to evaluate the field at,
-                         in degrees between -90 and 90 (negative for southern hemisphere).
-   * @param  longitute   Longitude of the location to evaluate the field at,
-                         in degrees between -180 and 180 (negative for western
-   hemisphere).
-   *
-   * @returns    The magnetic field vector in nT.
-   *
-   */
-
-  inline MagneticFieldVector get_wmm(const CoordinateSystemPtr Cs, const double year,
-                                     const LengthType altitude, const double latitude,
-                                     const double longitude);
-} // namespace corsika
-
-#include <corsika/detail/media/WMM.inl>

corsika/media/WorldMagneticModel.hpp

+#pragma once
+
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/geometry/PhysicalGeometry.hpp>
+
+#include <fstream>
+#include <string>
+
+namespace corsika {
+
+  /**
+   * A magnetic field calculated with the WMM or IGRF model.
+   */
+
+  class WorldMagneticModel {
+
+    /**
+     * Internal data structure for a single shell of the spherical harmonic
+     * parameterization.
+     */
+    struct ParameterLine {
+      int n;
+      int m;
+      double g;
+      double h;
+      double dg; //< time dependence of g in "per year"
+      double dh; //< time dependence of h in "per year"
+    };
+
+  public:
+    /**
+     * Construct a new World Magnetic Model object.
+     *
+     * @param center Center of Earth.
+     * @param data Data table to read.
+     */
+    WorldMagneticModel(Point const& center, std::string const& data = "GeoMag/WMM.COF");
+
+    /**
+     * Calculates the value of the magnetic field.
+     *
+     * @param  year        Year of the evaluation, between 2020 and 2025.
+     * @param  altitude    Height of the location to evaluate the field at,
+     *                    in km between -1 and 850.
+     * @param  latitude    Latitude of the location to evaluate the field at,
+     *                   in degrees between -90 and 90 (negative for southern hemisphere).
+     * @param  longitute   Longitude of the location to evaluate the field at,
+     *                      in degrees between -180 and 180 (negative for western
+     *   hemisphere).
+     *
+     * @returns    The magnetic field vector in nT.
+     */
+    MagneticFieldVector getField(double const year, LengthType const altitude,
+                                 double const latitude, double const longitude);
+
+  private:
+    Point center_;                                         //< Center of Earth.
+    std::map<int, std::vector<ParameterLine>> parameters_; //< epoch to Parameters map
+  };
+} // namespace corsika
+
+#include <corsika/detail/media/WorldMagneticModel.inl>
\ No newline at end of file

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/WMM.hpp>
+#include <corsika/media/WorldMagneticModel.hpp>
 
 #include <SetupTestTrajectory.hpp>
 #include <corsika/setup/SetupTrajectory.hpp>
@@ -25,7 +25,7 @@ using namespace corsika;
 
 TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
 
-  logging::set_level(logging::level::info);
+  logging::set_level(logging::level::trace);
 
   CoordinateSystemPtr const& gCS = get_root_CoordinateSystem();
   Point const gOrigin(gCS, {0_m, 0_m, 0_m});
@@ -74,18 +74,42 @@ TEST_CASE("UniformMagneticField w/ Homogeneous Medium") {
   }
 
   SECTION("WorldMagneticModel") {
-    // create earth magnetic field vector
-    MagneticFieldVector Earth_B_1 = get_wmm(gCS, 2022.5, 100_km, -80, -120);
-    MagneticFieldVector Earth_B_2 = get_wmm(gCS, 2022.5, 0_km, 0, 120);
-    MagneticFieldVector Earth_B_3 = get_wmm(gCS, 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_THROWS(WorldMagneticModel(gOrigin, "GeoMag/IDoNotExist.COF"));
+
+    {
+      WorldMagneticModel 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));
+    }
+    {
+      WorldMagneticModel wmm(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);
+      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));
+    }
   }
 }