Merge pull request #45 from gregstarr/develop

PR for issues #12, #25, #47, #28

README.rst

@@ -7,8 +7,9 @@ Overview
 This is a Python wrapper for the Apex fortran library by
 Emmert et al. [2010] [1]_, which allows converting between geodetic, modified
 apex, and quasi-dipole coordinates as well as getting modified apex and
-quasi-dipole base vectors (Richmond [1995] [2]_). MLT calculations are also
-included. The package is free software (MIT license).
+quasi-dipole base vectors (Richmond [1995] [2]_). The geodetic system used here
+is WGS84. MLT calculations are also included. The package is free software
+(MIT license).
 
 Quick start
 ===========

setup.cfg

@@ -52,7 +52,6 @@ addopts =
   --doctest-modules
   --doctest-glob='*.rst'
   -rxEfsw
-  --strict
   --ignore=docs/conf.py
   --ignore=setup.py
   --ignore=.eggs

setup.py

@@ -49,7 +49,7 @@ def read(*names, **kwargs):
         package_dir={'': 'src'},
         py_modules=[path.splitext(path.basename(pp))[0]
                     for pp in glob('src/*.py')],
-        package_data={'apexpy': ['apexsh.dat']},
+        package_data={'apexpy': ['apexsh.dat', 'igrf13coeffs.txt']},
         zip_safe=False,
         classifiers=[
             # complete classifier list:

src/apexpy/apex.py

@@ -56,7 +56,9 @@ class Apex(object):
 
     Notes
     -----
-    The calculations use IGRF-12 with coefficients from 1900 to 2020 [1]_.
+    The calculations use IGRF-13 with coefficients from 1900 to 2025 [1]_.
+
+    The geodetic reference ellipsoid is WGS84.
 
     References
     ----------
@@ -443,9 +445,9 @@ def _qd2apex_nonvectorized(self, qlat, qlon, height):
                 estr += '({:.3g}) for qlat {:.3g}'.format(self.refh, qlat)
                 raise ApexHeightError(estr)
 
-        alat = np.sign(qlat) * np.degrees(np.arccos(np.sqrt((self.RE +
-                                                             self.refh) /
-                                                            (self.RE + hA))))
+        sqlat = np.sign(qlat) if qlat != 0 else 1
+        alat = sqlat * np.degrees(np.arccos(np.sqrt((self.RE + self.refh) /
+                                                    (self.RE + hA))))
 
         return alat, alon
 
@@ -779,9 +781,9 @@ def basevectors_apex(self, lat, lon, height, coords='geo', precision=1e-10):
 
         Parameters
         ----------
-        lat, lon : (N,) array_like or float
-            Latitude
         lat : (N,) array_like or float
+            Latitude
+        lon : (N,) array_like or float
             Longitude
         height : (N,) array_like or float
             Altitude in km
@@ -800,11 +802,10 @@ def basevectors_apex(self, lat, lon, height, coords='geo', precision=1e-10):
 
         Returns
         -------
-        f1, f2 : (2, N) or (2,) ndarray
         f3, g1, g2, g3, d1, d2, d3, e1, e2, e3 : (3, N) or (3,) ndarray
 
-        Note
-        ----
+        Notes
+        -----
         `f3`, `g1`, `g2`, and `g3` are not part of the Fortran code
         by Emmert et al. [2010] [5]_. They are calculated by this
         Python library according to the following equations in
@@ -930,7 +931,8 @@ def set_epoch(self, year):
         # f2py
         self.year = np.float64(year)
         fa.loadapxsh(self.datafile, self.year)
-        fa.cofrm(self.year)
+        igrf_fn = os.path.join(os.path.dirname(__file__), 'igrf13coeffs.txt')
+        fa.cofrm(self.year, igrf_fn)
 
     def set_refh(self, refh):
         """Updates the apex reference height for all subsequent conversions.

src/fortranapex/fortranapex.pyf

@@ -3,54 +3,6 @@
 
 python module fortranapex ! in 
     interface  ! in :fortranapex
-        subroutine cofrm(date) ! in :fortranapex:magfld.f
-            real :: date
-            integer :: nmax
-            real dimension(255) :: gb
-            real dimension(225) :: gv
-            integer, optional :: ichg=-99999
-            common /magcof/ nmax,gb,gv,ichg
-        end subroutine cofrm
-        subroutine dypol(colat,elon,vp) ! in :fortranapex:magfld.f
-            real :: colat
-            real :: elon
-            real :: vp
-            integer :: nmax
-            real dimension(255) :: gb
-            real dimension(225) :: gv
-            integer :: ichg
-            common /magcof/ nmax,gb,gv,ichg
-        end subroutine dypol
-        subroutine feldg(ienty,glat,glon,alt,bnrth,beast,bdown,babs) ! in :fortranapex:magfld.f
-            integer intent(in) :: ienty
-            real(4) intent(in) :: glat
-            real(4) intent(in) :: glon
-            real(4) intent(in) :: alt
-            real(4) intent(out) :: bnrth
-            real(4) intent(out) :: beast
-            real(4) intent(out) :: bdown
-            real(4) intent(out) :: babs
-            integer :: nmax
-            real dimension(255) :: gb
-            real dimension(225) :: gv
-            integer :: ichg
-            common /magcof/ nmax,gb,gv,ichg
-        end subroutine feldg
-        subroutine gd2cart(gdlat,glon,alt,x,y,z) ! in :fortranapex:magfld.f
-            real :: gdlat
-            real :: glon
-            real :: alt
-            real :: x
-            real :: y
-            real :: z
-        end subroutine gd2cart
-        subroutine convrt(i,gdlat,alt,x1,x2) ! in :fortranapex:magfld.f
-            integer :: i
-            real :: gdlat
-            real :: alt
-            real :: x1
-            real :: x2
-        end subroutine convrt
         subroutine apex(date,dlat,dlon,alt,a,alat,alon,bmag,xmag,ymag,zmag,v) ! in :fortranapex:apex.f
             real :: date
             real :: dlat
@@ -223,8 +175,8 @@ python module fortranapex ! in
         end module apxshmodule
         subroutine loadapxsh(datafilenew,epochnew) ! in :fortranapex:apexsh.f90
             use apxshmodule
-            character*1000 :: datafilenew
-            real(kind=4) :: epochnew
+            character*1000 intent(in) :: datafilenew
+            real(kind=4) intent(in) :: epochnew
         end subroutine loadapxsh
         subroutine allocatearrays ! in :fortranapex:apexsh.f90
             use apxshmodule
@@ -303,6 +255,56 @@ python module fortranapex ! in
             real(kind=8) dimension(nmax + 1,mmax + 1),intent(out),depend(nmax,mmax) :: v
             real(kind=8) dimension(nmax + 1,mmax + 1),intent(out),depend(nmax,mmax) :: w
         end subroutine alfbasis
+        subroutine cofrm(date,filename) ! in :fortranapex:magfld.f
+            use igrf
+            real(kind=4) intent(in) :: date
+            character*1000 intent(in) :: filename
+            integer :: nmax
+            real dimension(255) :: gb
+            real dimension(225) :: gv
+            integer, optional :: ichg=-99999
+            common /magcof/ nmax,gb,gv,ichg
+        end subroutine cofrm
+        subroutine dypol(colat,elon,vp) ! in :fortranapex:magfld.f
+            real :: colat
+            real :: elon
+            real :: vp
+            integer :: nmax
+            real dimension(255) :: gb
+            real dimension(225) :: gv
+            integer :: ichg
+            common /magcof/ nmax,gb,gv,ichg
+        end subroutine dypol
+        subroutine feldg(ienty,glat,glon,alt,bnrth,beast,bdown,babs) ! in :fortranapex:magfld.f
+            integer intent(in) :: ienty
+            real intent(in) :: glat
+            real intent(in) :: glon
+            real intent(in) :: alt
+            real intent(out) :: bnrth
+            real intent(out) :: beast
+            real intent(out) :: bdown
+            real intent(out) :: babs
+            integer :: nmax
+            real dimension(255) :: gb
+            real dimension(225) :: gv
+            integer :: ichg
+            common /magcof/ nmax,gb,gv,ichg
+        end subroutine feldg
+        subroutine gd2cart(gdlat,glon,alt,x,y,z) ! in :fortranapex:magfld.f
+            real :: gdlat
+            real :: glon
+            real :: alt
+            real :: x
+            real :: y
+            real :: z
+        end subroutine gd2cart
+        subroutine convrt(i,gdlat,alt,x1,x2) ! in :fortranapex:magfld.f
+            integer :: i
+            real :: gdlat
+            real :: alt
+            real :: x1
+            real :: x2
+        end subroutine convrt
         subroutine makeapxsh(datafilein,epochgridin,nepochin,lmaxin,mmaxin,nmaxin) ! in :fortranapex:makeapexsh.f90
             use apxshmodule
             character*128 intent(in) :: datafilein

src/fortranapex/igrf.f90

@@ -27,10 +27,6 @@ subroutine read_igrf(filename_in,GYR,HYR,GT,HT,NEPO,NGHT,EPOCH,NMXE)
     !num_sh = NGHT-2*sqrt(real(NGHT))
     !write(*,*) num_sh
 
-    write(*,*) '----------------------------'
-    write(*,*) '    Read IGRF coefficeints  '
-    write(*,*) '----------------------------'
-
     ! Open IGRF file
     open(unit=100, file=filename_in,status='old',iostat=state)
     if (state /= 0) then
@@ -47,7 +43,6 @@ subroutine read_igrf(filename_in,GYR,HYR,GT,HT,NEPO,NGHT,EPOCH,NMXE)
     num_epochs=count([(s(i:i+3),i=1,len_trim(s))].eq.'IGRF')
     num_epochs=count([(s(i:i+3),i=1,len_trim(s))].eq.'DGRF')+num_epochs
     allocate(EPOCH(1:num_epochs))
-    write(*,*) 'Number of epochs read: ',num_epochs
     allocate(NMXE(1:num_epochs))
 
     ! Read epochs
@@ -107,7 +102,6 @@ subroutine read_igrf(filename_in,GYR,HYR,GT,HT,NEPO,NGHT,EPOCH,NMXE)
     GT =0.0d0
     HYR=0.0d0
     HT =0.0d0
-    write(*,*) 'ASSIGN COEFFICIENTS: ',NEPO,NGHT
     do e=1,NEPO+1
        do l=1,L_max
           if (e .le. NEPO) then

src/fortranapex/magfld.f

@@ -1,6 +1,6 @@
 C  FILE NAME: magfld.f
 
-      SUBROUTINE COFRM (DATE)
+      SUBROUTINE COFRM (DATE, FILENAME)
 C          Define the International Geomagnetic Reference Field (IGRF) as a
 C          scalar potential field using a truncated series expansion with
 C          Schmidt semi-normalized associated Legendre functions of degree n and
@@ -9,7 +9,8 @@ SUBROUTINE COFRM (DATE)
 C          rate after the last epoch.
 C
 C          INPUTS:
-C            DATE = yyyy.fraction (UT)
+C            DATE = yyyy.fraction (UT)\
+C            FILENAME = filename
 C          OUTPUTS (in comnon block MAGCOF):
 C            NMAX = Maximum order of spherical harmonic coefficients used
 C            GB   = Coefficients for magnetic field calculation
@@ -129,9 +130,8 @@ SUBROUTINE COFRM (DATE)
       ICHG = 1
 
 c          Load coefficients
-      FILENAME='igrf13coeffs.txt'
       if (.not. allocated(GYR)) then
-        call read_igrf(filename,GYR,HYR,GT,HT,NEPO,NGHT,EPOCH,NMXE)
+        call read_igrf(FILENAME,GYR,HYR,GT,HT,NEPO,NGHT,EPOCH,NMXE)
       endif
       NGH=NGHT*NEPO
 

tests/test_Apex.py

@@ -310,6 +310,16 @@ def test_convert_qd2apex():
                     apex_out.qd2apex(60, 15, height=100))
 
 
+def test_convert_qd2apex_at_equator():
+    """Test the quasi-dipole to apex conversion at the magnetic equator"""
+    apex_out = Apex(date=2000, refh=80)
+    elat, elon = apex_out.convert(lat=0.0, lon=0, source='qd', dest='apex',
+                                  height=120.0)
+    clat, clon = apex_out.convert(lat=0.001, lon=0, source='qd', dest='apex',
+                                  height=120.0)
+    assert_allclose([elat, elon], [clat, clon], atol=1e-4)
+
+
 def test_convert_qd2mlt():
     datetime = dt.datetime(2000, 3, 9, 14, 25, 58)
     apex_out = Apex(date=2000, refh=300)

tests/test_cmd.py

@@ -8,7 +8,9 @@
 import subprocess
 
 os.chdir(os.path.join(os.path.dirname(os.path.abspath(__file__)), '..'))
-outfile = 'tests/output.txt'
+outfile = os.path.join('tests', 'output.txt')
+infile = os.path.join('tests', 'test_convert.txt')
+singlefile = os.path.join('tests', 'test_convert_single_line.txt')
 
 
 def setup_function(function):
@@ -23,10 +25,10 @@ def setup_function(function):
 
 def test_module_invocation():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'apex', '2015',
-                             '--height', '300', '-i', 'tests/test_convert.txt',
-                             '-o', outfile])
+                             '--height', '300', '-i', infile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
+
     assert os.path.isfile(outfile)
     data = np.loadtxt(outfile)
     np.testing.assert_allclose(data, [[57.469547, 93.639816],
@@ -36,8 +38,7 @@ def test_module_invocation():
 
 def test_convert_YYYY():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'apex', '2015',
-                             '--height', '300', '-i', 'tests/test_convert.txt',
-                             '-o', outfile])
+                             '--height', '300', '-i', infile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
     assert os.path.isfile(outfile)
@@ -49,8 +50,7 @@ def test_convert_YYYY():
 
 def test_convert_YYYYMM():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'apex', '201501',
-                             '--height', '300', '-i', 'tests/test_convert.txt',
-                             '-o', outfile])
+                             '--height', '300', '-i', infile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
     assert os.path.isfile(outfile)
@@ -63,7 +63,7 @@ def test_convert_YYYYMM():
 def test_convert_YYYYMMDD():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'apex',
                              '20150101', '--height', '300', '-i',
-                             'tests/test_convert.txt', '-o', outfile])
+                             infile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
     assert os.path.isfile(outfile)
@@ -76,7 +76,7 @@ def test_convert_YYYYMMDD():
 def test_convert_YYYYMMDDHHMMSS():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'apex',
                              '20150101000000', '--height', '300', '-i',
-                             'tests/test_convert.txt', '-o', outfile])
+                             infile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
     assert os.path.isfile(outfile)
@@ -89,8 +89,7 @@ def test_convert_YYYYMMDDHHMMSS():
 def test_convert_single_line():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'apex',
                              '20150101000000', '--height', '300', '-i',
-                             'tests/test_convert_single_line.txt', '-o',
-                             outfile])
+                             singlefile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
     assert os.path.isfile(outfile)
@@ -122,8 +121,7 @@ def test_convert_refh():
 def test_convert_mlt():
     pipe = subprocess.Popen(['python', '-m', 'apexpy', 'geo', 'mlt',
                              '20150101000000', '--height', '300', '-i',
-                             'tests/test_convert_single_line.txt', '-o',
-                             outfile])
+                             singlefile, '-o', outfile])
     pipe.communicate()
     pipe.wait()
     assert os.path.isfile(outfile)