#10 IRBEM IGRF implementation producing incorrect values

[AaronH]

I have been doing some testing of the IGRF code included in IRBEM, and I believe that the values it is producing are incorrect. I have compared the IRBEM produced values with by my own code and the official IGRF code (provided by BGS) downloaded from https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html -- both my code and the official code agree, but the IRBEM code is often out by tens of nT.

For instance, using the geographic location lat=45 deg,lon=45 deg,alt=0 km, IRBEM produces the following field values:

BX: -35760.422335333125
BY: -31648.993464602798
BZ: -16675.949955399021
Per my code and the BGS code, the correct values are closer to:
Bx: -35780.284108
By: -31631.701466
Bz: -16636.035447

I have tried to determine exactly where this error is occurring, but to be honest the IRBEM IGRF code is rather opaque -- it is not a simple implementation of the series expansion formula from the IGRF website, and it is barely documented at all. Where it is documented, it appears to be out of date (igrf.f makes mention of conversion from geodetic coordinates, however no such conversion seems to be made, and the input coordinates are geographic cartesian). I do however know that the error is not due to the truncation of the Gauss coefficients at N=10, as I have tried artificially truncating my own code to N=10, and it does not explain the difference.

I realise that these errors are not hugely significant, but I feel that these should be addressed. I worry that if it is not simply due to a truncation of the coefficients, then it might be a coding error. I was unable to determine exactly what the code was doing, so I can't prove this. Part of the problem is that it is riddled with redundant code; for instance in the calculation of the Schmidt normalisation coefficients, there are the following lines of code:

        F0 = F0 * X * X / (4.D0 * X - 2.D0)
        F0 = F0 * (2.D0 * X - 1.D0) / X
        F = F0 * 0.5D0 * SQRT(2.D0)

It does not take much rearranging to realise that this is functionally equivalent to:

        F0 = F0 * X / 2.D0
        F = F0 / SQRT(2.D0)

There are several such redundancies, and it wouldn't surprise me if it was in one of these that the error has occurred.

My suggestion would be to scrap the entire code as it is and use the implementation from BGS (which is available under the MIT license, so can be included without issue), but I realise that could be quite a bit of work.

[AaronH]

I realise I forgot to specify that the above example was calculated for 2015/01/1, but the year used does not seem to matter -- the values are incorrect regardless.

[Paul OBrien]

I would guess this is a truncation issue. The ONERA group wrote this bit, I believe. If Sebastien weighs in, we might learn more about why it looks the way it does.

If you would like to modify the code to make it cleaner and to conform to the latest / official implementation that is welcome. Be mindful, however, that the truncation was made in order to maintain speed of execution. The library is intended to be a fast code and sometimes sacrifices accuracy to that end. If you do make a modification, please try to honor that intent. I.e., confirm that the new code runs approximately as fast as the original code.

Also, note that the OPTIONS array that we use as input to many functions is used to control numerical precision of the various routines. You might be able to shoehorn a truncation control parameter into that somehow. That would allow users working near the surface of the Earth to request higher order evaluations.

[AaronH]

Hi Paul,

On further investigation I think you're right. I had actually tested truncating the result to N=10 and was getting different results, but I now realise that's due to the way that IGRF interpolates the Gauss coefficients (if you don't specify the IGRF initialization frequency, it defaults to the middle of the year rather than the start of the year). I get values much closer if I emulate this.

I can try and modify the code to try and implement the full range of coefficients, although I don't have much of a timescale for this. I think it is possible to achieve performance parity or better while still implementing the full N=13, based on my own testing. With an improved algorithm I was able to get a 3-5x speed up compared to IRBEM when computing a large number of points, while still maintaining full N=13 accuracy.