Adding IGRF13 model.

eoxmagmod/eoxmagmod/__init__.py

@@ -140,7 +140,7 @@
 except ImportError:
     pass
 
-__version__ = '0.9.6'
+__version__ = '0.9.7'
 __author__ = 'Martin Paces ([email protected])'
 __copyright__ = 'Copyright (C) 2014 EOX IT Services GmbH'
 __licence__ = 'EOX licence (MIT style)'

eoxmagmod/eoxmagmod/data/__init__.py

@@ -44,6 +44,7 @@
 CHAOS_CORE_LATEST = CHAOS7_CORE_LATEST
 IGRF11 = join(_DIRNAME, 'igrf11coeffs.txt')
 IGRF12 = join(_DIRNAME, 'IGRF12.shc')
+IGRF13 = join(_DIRNAME, 'igrf13coeffs.txt')
 SIFM = join(_DIRNAME, 'SIFM.shc')
 LCS1 = join(_DIRNAME, 'LCS-1.shc')
 MF7 = join(_DIRNAME, 'MF7.shc')

eoxmagmod/eoxmagmod/magnetic_model/loader_igrf.py

@@ -34,9 +34,6 @@
 
 def load_model_igrf(path):
     """ Load model from an IGRF coefficient file.
-
-    Note: IGRF file format was discontinued after IGRF11. Starting by IGRF12
-    the IRFF models uses SHC format which is loaded by the SHC file loader.
     """
     return SphericalHarmomicGeomagneticModel(load_coeff_igrf(path))
 

eoxmagmod/eoxmagmod/magnetic_model/tests/coefficient_loaders.py

@@ -36,7 +36,7 @@
 from eoxmagmod.data import (
     CHAOS6_CORE_LATEST, CHAOS6_STATIC,
     CHAOS7_CORE_LATEST, CHAOS7_STATIC,
-    IGRF11, IGRF12, SIFM, WMM_2015,
+    IGRF11, IGRF12, IGRF13, SIFM, WMM_2015,
     EMM_2010_STATIC, EMM_2010_SECVAR,
     LCS1, MF7,
 )
@@ -186,6 +186,17 @@ class TestCoeffSIFM(TestCase, ShcTestMixIn):
     validity = decimal_year_to_mjd2000((2013.4976, 2015.4962))
 
 
+class TestCoeffIGRF13(TestCase, CoefficietLoaderTestMixIn):
+    class_ = SparseSHCoefficientsTimeDependentDecimalYear
+    degree = 13
+    min_degree = 1
+    validity = decimal_year_to_mjd2000((1900.0, 2025.0))
+
+    @staticmethod
+    def load():
+        return load_coeff_igrf(IGRF13)
+
+
 class TestCoeffIGRF12(TestCase, ShcTestMixIn):
     class_ = SparseSHCoefficientsTimeDependentDecimalYear
     path = IGRF12

eoxmagmod/eoxmagmod/magnetic_model/tests/model_loaders.py

@@ -34,9 +34,7 @@
 from numpy.random import uniform
 from numpy.testing import assert_allclose
 from eoxmagmod.magnetic_time import mjd2000_to_magnetic_universal_time
-from eoxmagmod.time_util import (
-    decimal_year_to_mjd2000, decimal_year_to_mjd2000_simple,
-)
+from eoxmagmod.time_util import decimal_year_to_mjd2000
 from eoxmagmod.magnetic_model.loader_shc import (
     load_model_shc, load_model_shc_combined,
 )
@@ -59,7 +57,7 @@
     EMM_2010_STATIC, EMM_2010_SECVAR, WMM_2015,
     CHAOS6_CORE_LATEST, CHAOS6_STATIC,
     CHAOS7_CORE_LATEST, CHAOS7_STATIC,
-    IGRF11, IGRF12, SIFM, LCS1, MF7,
+    IGRF11, IGRF12, IGRF13, SIFM, LCS1, MF7,
 )
 from eoxmagmod.magnetic_model.tests.data import (
     SWARM_MMA_SHA_2C_TEST_DATA,
@@ -475,6 +473,19 @@ def load(self):
         return load_model_shc(IGRF12, interpolate_in_decimal_years=True)
 
 
+class TestIGRF13(TestCase, SHModelTestMixIn):
+    reference_values = (
+        -15522.5, (30.0, 40.0, 8000.0),
+        (15259.57386772841, -159.00767967612023, -14015.952721753336)
+    )
+    degree = 13
+    min_degree = 1
+    validity = decimal_year_to_mjd2000((1900.0, 2025.0))
+
+    def load(self):
+        return load_model_igrf(IGRF13)
+
+
 class TestSIFM(TestCase, SHModelTestMixIn):
     reference_values = (
         5295.36, (30.0, 40.0, 8000.0),

eoxmagmod/eoxmagmod/magnetic_model/tests/parser_igrf.py

@@ -31,7 +31,7 @@
 from unittest import TestCase, main
 from numpy import abs as aabs
 from eoxmagmod.magnetic_model.parser_igrf import parse_igrf_file
-from eoxmagmod.data import IGRF11
+from eoxmagmod.data import IGRF11, IGRF13
 
 
 class TestIGRFParser(TestCase):
@@ -60,6 +60,13 @@ def test_parse_igrf_file_igrf11(self):
             "degree_max": 13,
         })
 
+    def test_parse_igrf_file_igrf13(self):
+        data = self.parse(IGRF13)
+        self._assert_valid(data, {
+            "degree_min": 1,
+            "degree_max": 13,
+        })
+
 
 if __name__ == "__main__":
     main()

eoxmagmod/eoxmagmod/tests/data/update_quasi_dipole_test_data.py

@@ -0,0 +1,96 @@
+#-------------------------------------------------------------------------------
+#
+#  Quasi-Dipole coordinates - test
+#
+# Author: Martin Paces <[email protected]>
+#
+#-------------------------------------------------------------------------------
+# Copyright (C) 2018 EOX IT Services GmbH
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies of this Software or works derived from this Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+# THE SOFTWARE.
+#-------------------------------------------------------------------------------
+#pylint: disable=missing-docstring
+
+import sys
+try:
+    # Python 2
+    from itertools import izip as zip
+except ImportError:
+    pass
+from numpy import array
+from eoxmagmod import mjd2000_to_decimal_year
+from eoxmagmod.quasi_dipole_coordinates import (
+    eval_mlt, eval_subsol, eval_qdlatlon_with_base_vectors,
+)
+
+EARTH_RADIUS = 6371.2 # km
+START_TIME = 0.0    # MJD2000 / 2000-01-01T00:00:00Z
+END_TIME = 7305.0   # MJD2000 / 2020-01-01T00:00:00Z
+
+
+def load_coordinates(file_in):
+
+    _ = next(file_in)
+
+    times, lats, lons, rads = [], [], [], []
+
+    for line in file_in:
+        time, _, lat, lon, rad = [float(v) for v in line.split()[:5]]
+        times.append(time)
+        lats.append(lat)
+        lons.append(lon)
+        rads.append(rad)
+
+    return array(times), array(lats), array(lons), array(rads)
+
+
+def update_test_data(file_out, time_mjd2000, lat, lon, rad):
+    """ Generate test dataset. """
+    time_decimal_year = mjd2000_to_decimal_year(time_mjd2000)
+
+    qdlat, qdlon, f11, f12, f21, f22, f__ = eval_qdlatlon_with_base_vectors(
+        lat, lon, rad, time_decimal_year
+    )
+    mlt = eval_mlt(qdlon, time_mjd2000)
+    sol_lat, sol_lon = eval_subsol(time_mjd2000)
+
+    header = [
+        "MJD2000", "DecimalYear", "Latitude", "Longitude", "Radius",
+        "QDLatitude", "QDLongitude", "F11", "F12", "F21", "F22", "F",
+        "MagneticLocalTime", "SubsolarLatitude", "SubsolarLongitude",
+    ]
+    records = zip(
+        time_mjd2000, time_decimal_year, lat, lon, rad,
+        qdlat, qdlon, f11, f12, f21, f22, f__,
+        mlt, sol_lat, sol_lon
+    )
+
+    file_out.write("\t".join(header) + "\n")
+    for record in records:
+        file_out.write("\t".join("%.14e" % value for value in record) + "\n")
+
+
+def main(filename):
+    with open(filename) as file_in:
+        time, lat, lon, rad = load_coordinates(file_in)
+    update_test_data(sys.stdout, time, lat, lon, rad)
+
+
+if __name__ == "__main__":
+    main(sys.argv[1])