diff --git a/README.rst b/README.rst index f289c686..a2b8a896 100644 --- a/README.rst +++ b/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 =========== diff --git a/setup.cfg b/setup.cfg index b3f75ce4..ed005def 100644 --- a/setup.cfg +++ b/setup.cfg @@ -52,7 +52,6 @@ addopts = --doctest-modules --doctest-glob='*.rst' -rxEfsw - --strict --ignore=docs/conf.py --ignore=setup.py --ignore=.eggs diff --git a/setup.py b/setup.py index d6ade3f1..97d88aca 100644 --- a/setup.py +++ b/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: diff --git a/src/apexpy/apex.py b/src/apexpy/apex.py index 7388120e..4f9e2e97 100644 --- a/src/apexpy/apex.py +++ b/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. diff --git a/src/fortranapex/igrf13coeffs.txt b/src/apexpy/igrf13coeffs.txt similarity index 100% rename from src/fortranapex/igrf13coeffs.txt rename to src/apexpy/igrf13coeffs.txt diff --git a/src/fortranapex/fortranapex.pyf b/src/fortranapex/fortranapex.pyf index fa4a6bd0..e639655f 100644 --- a/src/fortranapex/fortranapex.pyf +++ b/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 diff --git a/src/fortranapex/igrf.f90 b/src/fortranapex/igrf.f90 index 16cd7399..e997a6c8 100644 --- a/src/fortranapex/igrf.f90 +++ b/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 diff --git a/src/fortranapex/magfld.f b/src/fortranapex/magfld.f index d8311bcf..1240170a 100644 --- a/src/fortranapex/magfld.f +++ b/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 diff --git a/tests/test_Apex.py b/tests/test_Apex.py index 09092762..130094a8 100644 --- a/tests/test_Apex.py +++ b/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) diff --git a/tests/test_cmd.py b/tests/test_cmd.py index d795445d..20faa1da 100644 --- a/tests/test_cmd.py +++ b/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)