T01 IDL vs Python differences

[jameswilburlewis]

The IDL and Python outputs from the T01 model are differing by a few nanotesla for a large part of the 2007-03-23 THEMIS-A ephemeris which I'm testing with. The differences are occurring in regions with low field strength (sometimes 16x the reference field value), so I think this is more than just roundoff error.

See the attached plot: top two panels are the IDL and Python T01 field outputs; third panel is the difference, bottom panel shows the GSM positions in RE used as model inputs. The other parameters were: pdyn = 2.0, dsti = -30.0, yimf = 0.0, zimf = -5.0, g1 = 6.0, g2 = 10.0.

[madokatext]

The python T01 use different coefficients "c_sy" (in func "full_rc") and "a" (in func "t01") as the original FORTRAN77 T01. I think it's the main cause of the difference results between python and IDL T01. However, I don't know why they have different coefficients, since both of them have the comments ”LATEST MODIFICATIONS/BUGS REMOVED: JUNE 24, 2006“.

Another weird point is that the IDL and the FORTRAN77 version still have slightly differnent results, but much smaller difference than that between python and FORTRAN77. The IDL geopack should be compiled from FORTRAN77 codes, so there shouldn't be differences in theory.

The situation is so complicated in T01 that I finally gave up dealing with the bugs of it.

[jameswilburlewis]

Interesting! I wonder if Haje Korth has made a few changes in the Fortran T01 code in the process of developing the IDL version?

I'm going to try using ChatGPT-4 to translate the TA15 and TA16 Fortran models to Python. If that goes well, maybe I'll redo T01 the same way and see what happens.

[madokatext]

I did another test, if I keep python coefficients "c_sy" (in func "full_rc"), "a" (in func "t01"), "a1,a2,rrc1,dd1,rrc2,dd2,p1,r1,dr1,dla1,p2,r2,dr2,dla2,p3,r3,dr3" (in func "ap") the same with FORTRAN77, then they got close results (less than 1% difference). That means the python T01 code may not have serious bugs, the main difference of the results only comes from coefficients differences.

[jameswilburlewis]

@madokatext : After making the changes you suggested, I also see much improved results when comparing to IDL! The maximum error is now down to about 0.25 nT in the inner magnetosphere. I'll submit a PR, and hopefully we can get all the model corrections released to pypi soon.