Bug #1153

IGRF calculation different in sniffer than in mgd77list/mgd77magref

Added by Mike 4 months ago. Updated 2 months ago.

Status:NewStart date:2017-09-20
Priority:NormalDue date:
Assignee:-% Done:


Target version:Candidate for next bugfix release
Affected version:6.x-svn Platform:


Because I'm creating both a new mgd77 file and testing it along-track, I know that the sniffer should detect that reported and recomputed magnetic anomalies are identical. This is expected when you create a new archive file the same day as you test it with the sniffer.

This is not the case, however. Recomputed magnetic anomalies computed internally for regression analysis are different than reported so the output e77 file contains an error message (?-E-km1609-mag-07: Anomaly differs from mtf1-IGRF (m: 1 b: -10 rms: 10 r: 1 sig: 1 dec: 1). [Recompute])

I added some debug lines in the sniffer just beneath the old_anom[n] = ... line, like this:

old_anom[n] = D[i].number[MGD77_MAG];
fprintf (fpout,"date %lf old_anom %f mtf1 %f igrf %f new_anom %f old-new %f\n",date,D[i].number[MGD77_MAG],D[i].number[MGD77_MTF2-(int)mtf1],IGRF[MGD77_IGRF_F],new_anom[n],old_anom[n]-new_anom[n]);
fprintf (fpout,"%lf %lf 0 %.9f\n",D[i].number[MGD77_LONGITUDE],D[i].number[MGD77_LATITUDE],date);

After recompiling sniffer I run it on the test file km1609.mgd77 as follows:

mgd77sniffer km1609.mgd77 -R-177.4/-176.85/-6.4/-6.15 -De

The output looks like this

  1. Cruise km1609 ID km1609 MGD77 FILE VERSION: 20170920 N_RECS: 976
  2. Examined: Wed Sep 20 22:55:31 2017 by mtchandl
  3. Arguments: km1609.mgd77 -R-177.4/-176.85/-6.4/-6.15 -De
    N Errata table verification status
  4. mgd77manage applies corrections if the errata table is verified (toggle 'N' above to 'Y' after review)
  5. For instructions on E77 format and usage, see http://gmt.soest.hawaii.edu/mgd77/errata.php
  6. Verified by:
  8. Errata: Header
    ?-E-km1609-H01-04: Invalid Data Center File Number: () [km1609]
    ?-E-km1609-H11-06: Invalid Rightmost Longitude: (-176) [-177]
    Y-W-km1609-nav-00: Flagged 0.10 % of records with bad navigation
    Y-I-km1609-faa-06: Anomaly equivalent to gobs-IGF80 (m: 1 b: 0 rms: 0 r: 1 sig: 1 dec: 1)
    Y-I-km1609-faa-15: Anomaly differs from gobs-IGF80+eot (m: -0.989896 b: 2.96857 rms: 42.1107 r: 0.960914 sig: 1 dec: 0)
    ?-E-km1609-faa-16: Free-air anomalies may have been computed using IGF 1930 (m: 0.936866 b: 8.58291 rms: 9.92853 r: 0.999936 sig: 1 dec: 0). [Adjust to IGF 1980]
    date 2017.880462 old_anom -286.600000 mtf1 35669.600000 igrf 35945.102616 new_anom -275.600000 old-new -11.000000
    -177.351860 -6.368800 0 2017.880461635
    ?-E-km1609-mag-07: Anomaly differs from mtf1-IGRF (m: 1 b: -10 rms: 10 r: 1 sig: 1 dec: 1). [Recompute]
  9. Errata: Data
    ? km1609 2016-12-14T22:13:00.00 254 0-0-V GRAD: eot excessive
    ? km1609 2016-12-14T22:13:15.00 255 0-0-UVW GRAD: gobs, eot, faa excessive
    ? km1609 2016-12-14T22:16:15.00 273 0-0-OQ GRAD: mtf1, mag excessive
    ? km1609 2016-12-14T22:55:22.80 528 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T22:55:46.50 531 0-0-L GRAD: depth excessive
    Y km1609 2016-12-14T23:03:45.00 583 C-0-0 NAV: excessive speed
    ? km1609 2016-12-14T23:32:10.74 778 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:32:30.12 781 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:34:05.16 789 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:42:07.62 845 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:42:47.58 850 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:43:25.98 853 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:43:44.28 855 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:54:59.16 941 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:55:15.24 944 0-0-L GRAD: depth excessive
    ? km1609 2016-12-14T23:59:55.38 976 0-0-L GRAD: depth excessive

The two debug lines show position time and computed IGRF values of
-177.351860 -6.368800 0 2017.880461635 35945.102616

In mgd77magref the output is different:
echo "-177.351860 -6.368800 0 2017.880461635" | mgd77magref -Ft/0

In mgd77list the output is almost the same as magref:
mgd77list km1609 -Imct -Ftime,lon,lat,'mag!=NaN',igrf | head -n 1
2016-12-14T22:14:45 -177.352 -6.3688 35956.2

The difference between IGRF calculations is almost the same as the difference between reported and recomputed magnetic anomalies (e.g., old-new = -11.000000)

I wonder why sniffer is giving different values. Perhaps it is not reading the correct IGRF coefficients?

km1609.mgd77 - mgd77 file truncated to only 1000 rows (117 KB) Mike, 2017-09-20 23:07


#1 Updated by Joaquim 4 months ago

Hmm, they both call the same function, MGD77_igrf10syn(), so the difference has to lie on the input parameters. Can you debug the sniffer and find out all params send in to MGD77_igrf10syn()?
(Debuggers are MUCH better than fprintf's)

#2 Updated by Mike 4 months ago

I have't kept up with the debugging times, obviously, but I'll keep plugging away.

Looks like there's a problem with the date sniffer is passing to the MGD77_igrf10syn function, in this case. mgd77list shows the correct date of 2016-12-14, which in decimal is 2016.9536 (349 / 366).

mgd77list km1609 -Imct -Fytime,lon,lat,'mag!=NaN',igrf --FORMAT_FLOAT_OUT=%.9f -Imct | head -n 1
2016.953352212 -177.351860000 -6.368800000 35956.170229904

But in the sniffer output I posted yesterday, the date is nearly a year off at 2017.880461635.
I passed the same erroneous date to mgd77magref yesterday, though, and yielded a result that differs from the sniffer.

#3 Updated by Mike 4 months ago

The date calculation in sniffer appears to be off:

date = cal.year + cal.day_y / n_days + (cal.hour * GMT_HR2SEC_I + cal.min * GMT_MIN2SEC_I + cal.sec) * GMT_SEC2DAY;

This portion, (cal.hour * GMT_HR2SEC_I + cal.min * GMT_MIN2SEC_I + cal.sec) * GMT_SEC2DAY, gives days, but it needs to be divided by n_days before adding to cal.year and cal.day_y/n_days.

It gives a more reasonable result this way:

date = cal.year + cal.day_y / n_days + ((cal.hour * GMT_HR2SEC_I + cal.min * GMT_MIN2SEC_I + cal.sec) * GMT_SEC2DAY) / n_days;

Now sniffer debug lines output this:

date 2016.956282 cal.year 2016 cal.day_y 349 cal.hour 23 GMT_HR2SEC_I 3600 cal.min 59 GMT_MIN2SEC_I 60 cal.sec 0.000000 GMT_SEC2DAY 0.000012 old_anom -194.900000 mtf1 35670.500000 igrf 35865.406443 new_anom -195.000000 old-new 0.100000
-177.187970 -6.192660 0 2016.956282256

This date looks correct since it just a minute before the day changes to 350, and 350/366 is ~ .9563

For this same record, mgd77list outputs this:
mgd77list km1609 -Imct -Ftime,lon,lat,'mag!=NaN',igrf --FORMAT_FLOAT_OUT=%.9f -Imct | tail -n 1
2016-12-14T23:59:00 -177.187970000 -6.192660000 35865.439059134
mgd77list km1609 -Imct -Fytime,lon,lat,'mag!=NaN',igrf --FORMAT_FLOAT_OUT=%.9f -Imct | tail -n 1
2016.953550015 -177.187970000 -6.192660000 35865.439059134

Although most of the IGRF discrepancy is gone, mgd77list decimal time looks wrong. It's actually less than we'd expect for 2016-12-14T00:00:00, which would be 2016.9536. So apparently mgd77list also has a problem with calculating decimal year.

Interestingly, mgd77magref gives very different output for the same time, depending on whether decimal year is given:

% echo "-177.187970 -6.192660 0 2016.956282256" | mgd77magref -Ft/0
% echo "-177.187970 -6.192660 0 2016-12-14T23:59:00" | mgd77magref -Ft/0

Although decimal year looks correct, its IGRF value is quite different from mgd77list and the debugged sniffer.

#4 Updated by Mike 4 months ago

mgd77list uses the MGD77_cal_to_fyear function to calculate decimal year.

It is calculated slightly differently than in the sniffer:

return (cal→year + ((cal→day_y - 1.0) + (cal→hour * GMT_HR2SEC_I + cal→min * GMT_MIN2SEC_I + cal→sec) * GMT_SEC2DAY) / n_days);

I assume this formula, which subtracts a day from cal→day_y is the correct one. So I've replaced sniffer's date calculation with this function, i.e.,

date = MGD77_cal_to_fyear (GMT, &cal); /* Get date as decimal year */

The two programs now give identical igrf output

#5 Updated by Paul 2 months ago

  • Target version set to Candidate for next bugfix release
  • Affected version changed from 5.x-svn to 6.x-svn

Where did you update and what version?

Also available in: Atom PDF