search

spedas / pyspedas

Python-based Space Physics Environment Data Analysis Software
https://pyspedas.readthedocs.io/
MIT License
143 stars 58 forks source link

Runtime warnings in geopack t01 and t96 models #830

Closed jameswilburlewis closed 2 months ago

jameswilburlewis commented 2 months ago

I haven't noticed these before, but under python 3.9 and 3.12, I'm seeing some runtime warnings from the t01 and t96 models:

t96:

15-Apr-24 09:05:44: /Users/jwl/PycharmProjects/pyspedas/venv/lib/python3.9/site-packages/geopack/geopack.py:167: RuntimeWarning: invalid value encountered in cast
  irp3 = np.int64(r+2)

15-Apr-24 09:05:44: /Users/jwl/PycharmProjects/pyspedas/venv/lib/python3.9/site-packages/geopack/geopack.py:168: RuntimeWarning: divide by zero encountered in scalar divide
  nm = np.int64(3+30/irp3)

15-Apr-24 09:05:44: /Users/jwl/PycharmProjects/pyspedas/venv/lib/python3.9/site-packages/geopack/geopack.py:168: RuntimeWarning: invalid value encountered in cast
  nm = np.int64(3+30/irp3)

t01:

15-Apr-24 09:05:24: /Users/jwl/PycharmProjects/pyspedas/venv/lib/python3.9/site-packages/geopack/geopack.py:167: RuntimeWarning: invalid value encountered in cast
  irp3 = np.int64(r+2)

15-Apr-24 09:05:24: /Users/jwl/PycharmProjects/pyspedas/venv/lib/python3.9/site-packages/geopack/geopack.py:168: RuntimeWarning: divide by zero encountered in scalar divide
  nm = np.int64(3+30/irp3)

15-Apr-24 09:05:24: /Users/jwl/PycharmProjects/pyspedas/venv/lib/python3.9/site-packages/geopack/geopack.py:168: RuntimeWarning: invalid value encountered in cast
  nm = np.int64(3+30/irp3)

The relevant code in geopack.py is in the igrf_geo routine:

    # In this new version, the optimal value of the parameter nm (maximal order of the spherical
    # harmonic expansion) is not user-prescribed, but calculated inside the subroutine, based
    # on the value of the radial distance r:
    irp3 = np.int64(r+2)
    nm = np.int64(3+30/irp3)
    if nm > 13: nm = 13
    k = nm+1

So I guess something is wrong with an R value being passed in (maybe NaN or inf?), the np.int64 cast fails returning 0, then you get the divide-by-zero error in the next line.

jameswilburlewis commented 2 months ago

Actually, the problem is with the pyspedas geopack test suite, not geopack itself. After interpolating the positions onto the proton_speed support data, there are NaNs in the position variable getting passed to geopack.