Bug report on t96.py

[madokatext]

Hello! I found that the func t96 of geopack python version has very different results with the IDL version at some places, such as dayside magnetosphere. But at other places, they get exactly the same results. I began to doubt some parts of the codes have some problems.

After long time of debugging, I found two bugs:

1. t96.py: func r2_birk: line 1410: "r2outer" should be "r2sheet".
2. t96.py: func birk1tot_02: line 844: "[x1,y1,z1,ps]" should be "[x2,y2,z2,ps]".

The IDL version of geopack is directly compiled from the original FORTRAN77 codes, so I think its results are correct. I diecrtly compared the FORTRAN77 codes with the python codes, and found the two bugs above.

[madokatext]

3. t96.py: func birk1tot_02: line 731: "&" should be "|".

[tsssss]

Thanks for catching the bugs. That's really helpful since I mostly use it in the nightside. Should I wait until you are done with debugging t96 then update the package?

---

From: madokatext @.> Sent: Friday, September 15, 2023 6:57 AM To: tsssss/geopack @.> Cc: Subscribed @.***> Subject: Re: [tsssss/geopack] Bug report on t96.py (Issue #18)

1. t96.py: func birk1tot_02: line 731: "&" should be "|".

— Reply to this email directly, view it on GitHubhttps://github.com/tsssss/geopack/issues/18#issuecomment-1721327737, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADBMLBTXC55PY4W7VBIRDPDX2RNEZANCNFSM6AAAAAA4YMDRFI. You are receiving this because you are subscribed to this thread.Message ID: @.***>

[madokatext]

I think you can update it right now. I haven't met problems using the package after that.

[jameswilburlewis]

I forked the repo and tried implementing the fixes @madokatext suggested. The agreement between IDL and Python is better now, but there are still some significant differences (a couple of nanotesla, but in a region with a much weaker field), and I think it's not yet completely fixed. For test data, I'm using the THEMIS-A ephemeris for 2007-03-23 (one minute cadence), with pdyn=2.0, dsti=-30.0, yimf=0.0, zimf=-5.0. (Probably not at all close to the actual conditions...) In IDL SPEDAS, I used the 'exact_tilt_times' keyword to match the Python geopack behavior.

Here's a plot of what I'm seeing after including madokatext's patches:

Top panel: IDL SPEDAS T96 output Next panel: Python T96 output Third panel: Difference between IDL and Python model outputs Bottom panel: GSM coordinates in Re (X,Y,Z, total)

[madokatext]

That's another bug, I've just corrected it. @jameswilburlewis

4. t96.py: func ringcurr96: line 425: "zs" should be "zsww", and add "zs=zsww" between line 440 and 441.

[jameswilburlewis]

Thanks @madokatext! With that additional fix, the results match to within 0.25 nT near perigee, and a relative tolerance of .001 everywhere.