Memory Error/Extrapolation bug

CHANGELOG.rst

@@ -8,6 +8,8 @@ Changelog
 * Added unit tests for todays date, ensuring that `apex.dat` is current
 * Added cron unit test to GitHub Actions CI
 * Added a logo
+* Correct indexing bug in Fortran source that was causing array overflow and
+  memory errors for extrapolated years beyond the latest formal IGRF fit.
 
 2.0.0 (2022-12-09)
 ------------------

apexpy/tests/test_Apex.py

@@ -1765,3 +1765,75 @@ def test_get_height_along_fieldline(self, lat, height):
         assert abs(height - fheight) < 1.0e-7, \
             "bad height calculation: {:.7f} != {:.7f}".format(height, fheight)
         return
+
+
+class TestApexMethodExtrapolateIGRF(object):
+    """Test the Apex methods on a year when IGRF must be extrapolated.
+
+    Notes
+    -----
+    Extrapolation should be done using a year within 5 years of the latest IGRF
+    model epoch.
+
+    """
+
+    def setup_method(self):
+        """Initialize all tests."""
+        self.apex_out = apexpy.Apex(date=2025, refh=300)
+        self.in_lat = 60
+        self.in_lon = 15
+        self.in_alt = 100
+        self.in_time = dt.datetime(2024, 2, 3, 4, 5, 6)
+        return
+
+    def teardown_method(self):
+        """Clean up after each test."""
+        del self.apex_out, self.in_lat, self.in_lon, self.in_alt
+        return
+
+    @pytest.mark.parametrize("method_name, out_comp",
+                             [("geo2apex",
+                               (56.25343704223633, 92.04932403564453)),
+                              ("apex2geo",
+                               (53.84184265136719, -66.93045806884766,
+                                3.6222547805664362e-06)),
+                              ("geo2qd",
+                               (56.82968521118164, 92.04932403564453)),
+                              ("apex2qd", (60.498401178276744, 15.0)),
+                              ("qd2apex", (59.49138097045895, 15.0))])
+    def test_method_scalar_input(self, method_name, out_comp):
+        """Test the user method against set values with scalars.
+
+        Parameters
+        ----------
+        method_name : str
+            Apex class method to be tested
+        out_comp : tuple of floats
+            Expected output values
+
+        """
+        # Get the desired methods
+        user_method = getattr(self.apex_out, method_name)
+
+        # Get the user output
+        user_out = user_method(self.in_lat, self.in_lon, self.in_alt)
+
+        # Evaluate the user output
+        np.testing.assert_allclose(user_out, out_comp, rtol=1e-5, atol=1e-5)
+
+        for out_val in user_out:
+            assert np.asarray(out_val).shape == (), "output is not a scalar"
+        return
+
+    def test_convert_to_mlt(self):
+        """Test conversion from mlon to mlt with scalars."""
+
+        # Get user output
+        user_out = self.apex_out.mlon2mlt(self.in_lon, self.in_time)
+
+        # Set comparison values
+        out_comp = 23.955474853515625
+
+        # Evaluate user output
+        np.testing.assert_allclose(user_out, out_comp, rtol=1e-5, atol=1e-5)
+        return

docs/installation.rst

@@ -131,6 +131,8 @@ To download an artifact:
 
 To install, use ``pip install .``
 
+.. _installation-build:
+
 Build from Source
 ^^^^^^^^^^^^^^^^^
 

docs/maintenance.rst

@@ -16,74 +16,48 @@ Updating IGRF
 The `International Geomagnetic Reference Field <https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html>`_
 is regularly updated to reflect the most recent changes to the Terrestrial
 magnetic field. apexpy currently uses IRGF-13 coefficients, which are provided
-in the ``apexpy/src/apexpy/igrf13coeff.txt`` file. To change or update the
-magnetic field coefficients used by apexpy, you need to update the python code.
-Assuming your new coefficient file has the same format, the process is simple:
-
-1. Clone the repository or your fork of the repository (see :ref:`contributing`)
-2. Update ``apexpy/src/apexpy/apex.py`` variable ``igrf_fn`` by setting
-   it equal to the new target filename
-3. Install the python package from the command line
-   (see :ref:`installation-cmd`)
-
-Updating ``apexsh.dat``
------------------------
+in the ``apexpy/apexpy/igrf13coeff.txt`` file. To change or update the
+magnetic field coefficients used by apexpy, you need to update the python code,
+then rerun the fortran program that builds ``apexpy/apexpy/apexsh.dat``. This 
+is what makes apexpy performant.   For more details, see Emmert et al. [2010] [1]_.
 
-After updating the IGRF coefficients file the ``apexsh.dat`` file needs to be
-rebuilt. This file is what makes apexpy performant. For more details, see
-Emmert et al. [2010] [1]_.
-
-Updating ``apexsh.dat`` is done by modifying and building the fortran source
-code in the ``apexpy/src/fortranapex`` directory. Working in that directory:
+Assuming your new coefficient file has the same format, the process is simple:
 
-1. Make sure a copy of the latest IGRF coefficient file is present in the
-   selfsame directory.
-2. Modify the ``igrffilein`` in ``checkapexsh.f90`` to the name of the IGRF
-   coefficient file (``igrf13coeff.txt``, for example).
-3. Modify ``checkapexsh.f90`` by adding the next 5 year epoch to the
+1. Clone the repository or your fork of the repository (see :ref:`contributing`).
+2. Update ``apexpy/apexpy/apex.py`` variable ``igrf_fn`` by setting
+   it equal to the new IGRF coefficient filename (``igrf13coeff.txt``, for example).
+3. In ``apexpy/fortranapex/checkapexsh.f90``, update the variable ``igrffilein``
+   to the new IGRF coefficent filename.  Relative paths are allowed.
+4. Modify ``checkapexsh.f90`` by adding the next 5 year epoch to the
    ``epochgrid`` variable and updating the ``nepochgrid`` variable as
    necessary. For example, if the newest IGRF coefficients are good up to 2025
    and ``epochgrid`` only has up to the year 2020, then add 2025 to
    ``epochgrid`` and then increment ``nepochgrid`` by 1.
-4. Build the ``apextest`` binary by running the ``make`` command.
-5. Copy the IGRF coefficient file (``cp ../apexpy/igrf13coeff.txt``) into the
-   current directory.
-6. Execute the ``apextest`` binary to generate the new ``apexsh.dat`` file.
-7. Copy the new ``apexsh.dat`` file to the ``apexpy/src/apexpy`` directory.
-
-After updating the ``apexsh.dat`` file, some of the unit tests and the
-documentation examples in the README and ``apexpy/docs/examples`` directory
-will need to be updated as well.
+5. Execute the ``apextest`` binary to generate the new ``apexsh.dat`` file.
+6. Update the unit tests in the class ``TestApexMethodExtrapolateIGRF`` in 
+   ``apexpy/apexpy/tests/test_Apex.py`` so that they check the methods are working
+   correctly with dates after the latest IGRF epoch (i.e., if the latest epoch is
+   2020, set the test to initialize with the year 2025).  You will have to update
+   the hard-coded confirmation values used by these tests.
+7. Commit all changes and create a pull request on GitHub to integrate your 
+   branch with updated IGRF into the main repository.
 
 Modifying Fortran Source
 ------------------------
 When modifying the fortran source code, it can be helpful to run a preliminary
-validation of the fortran output independent of the python wrapper.
+validation of the fortran output independent of the python wrapper. This should
+be done within the ``apexpy/fortranapex`` directory.
 
 1. Remove any existing binaries by running the ``make clean`` command.
 2. Build the ``apextest`` binary by running the ``make`` command.
-3. Copy the IGRF coefficient file (``cp ../apexpy/igrf13coeff.txt``) into the
-   current directory.
-4. Execute the ``apextext`` binary.
-5. Confirm the output printed to the screen matches the test output shown in
-   the comment blot at the bottom of ``checkapexsh.f90``.
-
-The output may not match the test output exactly due to floating point errors
-and improvements in the precision of the calculation.
-
-After updating the fortran source code, the signature file must be recreated so
-the python wrapper works correctly.  It is also a good idea to update
-``apexsh.dat`` following the instructions above.  Use `f2py <https://numpy.org/doc/stable/f2py/>`_
-to create a new signature file from the ``apexpy/src/fortranapex`` directory.
-::
-
-  f2py -m fortranapex apexsh.f90 igrf.f90 apex.f90 magfld.f90 makeapexsh.f90 -h fortranapex.pyf --overwrite-signature
-
-
-This will create the file ``fortranapex.pyf``.  Then reinstall the package with
-``pip`` from the root directory.  If the modifications involved adding or
-removing fortran source files, modify the list of extension sources in
-``setup.py``.
+3. Execute the ``apextext`` binary.
+4. Confirm the output printed to the screen matches the test output shown in
+   the comment block at the bottom of ``checkapexsh.f90``. The output may not match 
+   the test output exactly due to floating point errors and improvements in the 
+   precision of the calculation.
+5. If the modifications involved adding or removing fortran source files, modify the 
+   list of extension sources in ``setup.cfg``.
+6. Rebuild and install apexpy following the instructions in :ref:`installation-build`.
 
 Updating tests and style standards
 -----------------------------------

fortranapex/Makefile

@@ -6,7 +6,7 @@
 FC 	= gfortran  # mpif90, ifort
 LD      = $(FC)
 LDFLAGS = #-ipo
-FFLAGS 	= -Wall -O3 -fPIC #-ipo -check bounds -check pointer -check uninit
+FFLAGS 	= -Wall -O3 -fPIC -fcheck=all -fbounds-check -pedantic -Wextra #-ipo -check bounds -check pointer -check uninit
 SFLAGS  = -static #-libgfortran -static-libgcc (AGB, static should work?)
 STATIC  = 0
 

fortranapex/checkapexsh.f90

@@ -27,8 +27,8 @@ program checkapexsh
   integer(4), parameter :: nepochgrid=26
   ! integer(4), parameter :: nepochgrid=3       ! For testing/debugging
   integer(4)            :: lmax=3, nmmax=6
-  character(1000)       :: apexshfile='apexsh.dat'
-  character(len=1000)   :: igrffilein='igrf13coeffs.txt'
+  character(1000)       :: apexshfile='../apexpy/apexsh.dat'
+  character(len=1000)   :: igrffilein='../apexpy/igrf13coeffs.txt'
   real(8)               :: epochgrid(0:nepochgrid - 1)
   real(8)               :: epoch
   real(8)               :: glat, glon, alt, hr, prec, error

fortranapex/igrf.f90

@@ -48,16 +48,19 @@ subroutine read_igrf(filename_in)
     end do
 
     ! Read epochs
-    num_epochs = count([(s(i:i + 3), i = 1, len_trim(s))]== 'IGRF')
-    num_epochs = count([(s(i:i + 3), i = 1, len_trim(s))]== 'DGRF') + num_epochs
+    num_epochs = 0
+    do i=1, len_trim(s)
+      if ((s(i:i+3) == 'IGRF').or.(s(i:i+3) == 'DGRF')) then
+        num_epochs = num_epochs + 1
+      end if
+    end do
     allocate(epoch(1:num_epochs))
     allocate(nmxe(1:num_epochs))
 
     ! Read epochs
-    read(100,*, iostat = state) s
-    do i = 1, num_epochs
-       epoch(i) = 1900.0 + real(i - 1) * 5.0d0
-    end do
+    read(unit = 100, fmt = '(A)', iostat = state) s
+    read(s(8:), *) epoch
+
 
     ! Number of coefficients
     do i = 1, num_epochs
@@ -89,7 +92,7 @@ subroutine read_igrf(filename_in)
     call fseek(100, offset, 0, state)
 
     ! Assign the variables for Gauss coefficients
-    allocate(g(1:num_sh, 1:num_epochs))
+    allocate(g(1:num_sh, 1:num_epochs+1))
     allocate(nm(1:num_sh, 2))
 
     ! Read coefficients

fortranapex/magfld.f90

@@ -172,11 +172,11 @@ subroutine cofrm(date, filename)
       call exit(1)
     end if
 
-    ! Set the date and time
-    iy = 1
-    do while (date > epoch(iy))
-      iy = iy + 1
+    ! Set the date and time index
+    do iy = 1, nepo
+      if (date < epoch(iy)) exit
     end do
+    iy = iy - 1
 
     ngh = nght * nepo
     nmax1 = int(nmxe(iy))