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

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

@@ -58,6 +58,8 @@ class Apex(object):
     -----
     The calculations use IGRF-12 with coefficients from 1900 to 2020 [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
 
@@ -800,7 +802,6 @@ 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
@@ -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