d2dtf(), dtf2d(), utctai(), dat(): Inexact table and method of detecting jump

[mrhso]

The pre-1972 _changes are approximated to $10^{-7}\kern{3 mu}\text{s}$, leading to strange results.

Example (After #91 fixed)

>>> from astropy.time import Time, TimeDelta
>>> t = Time('1971-12-31 23:59:60.107757999', format = 'iso', scale = 'utc', precision = 9)
>>> t
<Time object: scale='utc' format='iso' value=1971-12-31 23:59:60.107757999>
>>> t.tai
<Time object: scale='tai' format='iso' value=1972-01-01 00:00:10.000000002>

In theory, it should be earlier than 1972-01-01 00:00:10, for 1972-01-01 00:00:00 UTC = 1972-01-01 00:00:10 TAI and the step was -0.1077580 s.

And the method of detecting jump is also inexact.

/* Separate TAI-UTC change into per-day (DLOD) and any jump (DLEAP). */
   dlod = 2.0 * (dat12 - dat0);
   dleap = dat24 - (dat0 + dlod);

The end of a day may not be 24:00:00 but 24:00:00+dleap.

So $\text{dleap}=\text{dat24}-\left(\text{dat0}+\dfrac{24\kern{3 mu}\text{h}+\text{dleap}}{24\kern{3 mu}\text{h}}\text{dlod}\right)$, that is, $\text{dleap}=\dfrac{\text{ERFA\_DAYSEC}}{\text{ERFA\_DAYSEC}+\text{dlod}}(\text{dat24}-(\text{dat0}+\text{dlod}))$.

Interestingly, this inexact method just matches $10^{-7}\kern{3 mu}\text{s}$ precision, dleap is calculated correctly (just a coincidence).

So the table _changes and the method of detecting jump may need to be modified at the same time.

Date	TAI-UTC
1960-01-01	$\dfrac{70890899507113}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-37300)$
1961-01-01	$\dfrac{71140899510863}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-37300)$
1961-08-01	$\dfrac{68640899473363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-37300)$
1962-01-01	$\dfrac{92292899473363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{351}{312500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-37665)$
1963-11-01	$\dfrac{97292899538363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{351}{312500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-37665)$
1964-01-01	$\dfrac{162006499538363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1964-04-01	$\dfrac{167006499613363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1964-09-01	$\dfrac{172006499688363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1965-01-01	$\dfrac{177006499763363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1965-03-01	$\dfrac{182006499838363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1965-07-01	$\dfrac{187006499913363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1965-09-01	$\dfrac{192006499988363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{62500}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-38761)$
1966-01-01	$\dfrac{215658499988363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{31250}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-39126)$
1968-02-01	$\dfrac{210658499838363}{50000000000000}\kern{3 mu}\text{s}+\dfrac{81}{31250}\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-39126)$

/* Dates and Delta(AT)s */
   static const eraLEAPSECOND _changes[] = {
      { 1960,  1,  1.41781799014226 },
      { 1961,  1,  1.42281799021726 },
      { 1961,  8,  1.37281798946726 },
      { 1962,  1,  1.84585798946726 },
      { 1963, 11,  1.94585799076726 },
      { 1964,  1,  3.24012999076726 },
      { 1964,  4,  3.34012999226726 },
      { 1964,  9,  3.44012999376726 },
      { 1965,  1,  3.54012999526726 },
      { 1965,  3,  3.64012999676726 },
      { 1965,  7,  3.74012999826726 },
      { 1965,  9,  3.84012999976726 },
      { 1966,  1,  4.31316999976726 },
      { 1968,  2,  4.21316999676726 },
      { 1972,  1, 10.0              },
      { 1972,  7, 11.0              },
      { 1973,  1, 12.0              },
      { 1974,  1, 13.0              },
      { 1975,  1, 14.0              },
      { 1976,  1, 15.0              },
      { 1977,  1, 16.0              },
      { 1978,  1, 17.0              },
      { 1979,  1, 18.0              },
      { 1980,  1, 19.0              },
      { 1981,  7, 20.0              },
      { 1982,  7, 21.0              },
      { 1983,  7, 22.0              },
      { 1985,  7, 23.0              },
      { 1988,  1, 24.0              },
      { 1990,  1, 25.0              },
      { 1991,  1, 26.0              },
      { 1992,  7, 27.0              },
      { 1993,  7, 28.0              },
      { 1994,  7, 29.0              },
      { 1996,  1, 30.0              },
      { 1997,  7, 31.0              },
      { 1999,  1, 32.0              },
      { 2006,  1, 33.0              },
      { 2009,  1, 34.0              },
      { 2012,  7, 35.0              },
      { 2015,  7, 36.0              },
      { 2017,  1, 37.0              }
   };

[mhvk]

@mrhso - thanks for reporting! And even more for suggesting fixes. Since this is a problem in SOFA, I think we will have to give specific examples using the ERFA/SOFA routines (just to show it is not astropy wrapping it incorrectly). It looks like you already have been changing the code, so it may be easy for you to provide such examples to forward to SOFA (if not, I can do it).

[mrhso]

@mrhso - thanks for reporting! And even more for suggesting fixes. Since this is a problem in SOFA, I think we will have to give specific examples using the ERFA/SOFA routines (just to show it is not astropy wrapping it incorrectly). It looks like you already have been changing the code, so it may be easy for you to provide such examples to forward to SOFA (if not, I can do it).

#include <sofa.h>
#include <stdio.h>

int main() {
    int j, iy, im, id, ihmsf[4];
    double tai1, tai2;

    j = iauD2dtf("UTC", 9, 2441317.0, 0.49999999999998834, &iy, &im, &id, ihmsf);
    printf("%d %d %d %d %d %d %d\n", iy, im, id, ihmsf[0], ihmsf[1], ihmsf[2], ihmsf[3]);
    j = iauUtctai(2441317.0, 0.49999999999998834, &tai1, &tai2);
    j = iauD2dtf("TAI", 9, tai1, tai2, &iy, &im, &id, ihmsf);
    printf("%d %d %d %d %d %d %d\n", iy, im, id, ihmsf[0], ihmsf[1], ihmsf[2], ihmsf[3]);

    return 0;
}

It also shows #91.

(SOFA 20210512) Result:

1971 12 31 23 59 59 999999999
1972 1 1 0 0 10 2

(SOFA 20210512 with d2dtf.c modified)

      leap = (fabs(dleap) > 0.00390625);

Result:

1971 12 31 23 59 60 107757999
1972 1 1 0 0 10 2

[mrhso]

(SOFA 20210512 with modified) d2dtf.c

   double a1, b1, fd, dat0, dat12, w, dat24, dlod, dleap;
…
   /* Any sudden change in TAI-UTC (seconds). */
      dlod = 2.0 * (dat12 - dat0);
      dleap = (dat24 - (dat0 + dlod)) * (DAYSEC / (DAYSEC + dlod));

   /* If leap second day, scale the fraction of a day into SI. */
      leap = (fabs(dleap) > 0.00390625);
      if (leap) fd += fd * dleap/DAYSEC;

dat.c

/* Dates and Delta(AT)s */
   static const struct {
      int iyear, month;
      double delat;
   } changes[] = {
      { 1960,  1,  1.41781799014226 },
      { 1961,  1,  1.42281799021726 },
      { 1961,  8,  1.37281798946726 },
      { 1962,  1,  1.84585798946726 },
      { 1963, 11,  1.94585799076726 },
      { 1964,  1,  3.24012999076726 },
      { 1964,  4,  3.34012999226726 },
      { 1964,  9,  3.44012999376726 },
      { 1965,  1,  3.54012999526726 },
      { 1965,  3,  3.64012999676726 },
      { 1965,  7,  3.74012999826726 },
      { 1965,  9,  3.84012999976726 },
      { 1966,  1,  4.31316999976726 },
      { 1968,  2,  4.21316999676726 },
      { 1972,  1, 10.0              },
      { 1972,  7, 11.0              },
      { 1973,  1, 12.0              },
      { 1974,  1, 13.0              },
      { 1975,  1, 14.0              },
      { 1976,  1, 15.0              },
      { 1977,  1, 16.0              },
      { 1978,  1, 17.0              },
      { 1979,  1, 18.0              },
      { 1980,  1, 19.0              },
      { 1981,  7, 20.0              },
      { 1982,  7, 21.0              },
      { 1983,  7, 22.0              },
      { 1985,  7, 23.0              },
      { 1988,  1, 24.0              },
      { 1990,  1, 25.0              },
      { 1991,  1, 26.0              },
      { 1992,  7, 27.0              },
      { 1993,  7, 28.0              },
      { 1994,  7, 29.0              },
      { 1996,  1, 30.0              },
      { 1997,  7, 31.0              },
      { 1999,  1, 32.0              },
      { 2006,  1, 33.0              },
      { 2009,  1, 34.0              },
      { 2012,  7, 35.0              },
      { 2015,  7, 36.0              },
      { 2017,  1, 37.0              }
   };

utctai.c

/* Separate TAI-UTC change into per-day (DLOD) and any jump (DLEAP). */
   dlod = 2.0 * (dat12 - dat0);
   dleap = (dat24 - (dat0 + dlod)) * (DAYSEC / (DAYSEC + dlod));

Result:

1971 12 31 23 59 60 107757999
1972 1 1 0 0 9 999999999

[mhvk]

@PTW19 - while you happen to have been looking at our ERFA repository, do you mind if I draw your attention to this one, about leap seconds? Here it looks to me like @mrhso has uncovered a genuine bug; some test code and suggested fixes are provided in the thread.

p.s. @mrhso - apologies for not reacting to this earlier; your examples are very useful!!

[PTW19]

The pre-1972 changes are approximated to 1e-7, leading to strange results.

This is how they were defined: see Table 3 in the Explanatory Supplement (3rd edition). Consequently things don't quite tie together, up to and including the final transition to the leap second era as 1972 Jan 1 begins.

Another factor in the strange-looking results is that 9 digits for ndp in the iauD2dtf call is pushing it. If you drop back to 8 you do get the expected result:

1972 1 1 0 0 0 0 1972 1 1 0 0 10 0

I'm aware the D2DTF comments say 9 should be safe, but evidently rounding errors have started to creep in at that point on the platforms we are both using.

As a first step, please see if these comments explain the reported effects.

[mrhso]

This is how they were defined: see Table 3 in the Explanatory Supplement (3rd edition). Consequently things don't quite tie together, up to and including the final transition to the leap second era as 1972 Jan 1 begins.

TimeSteps.history, the table of the frequency offsets and steps which should be the original definition, is available from the IERS EOC. Then combined with the comment of dat.c:

**  2) The 5ms timestep at 1961 January 1 is taken from 2.58.1 (p87) of
**     the 1992 Explanatory Supplement.

Obtaining a complete table is possible.	Date	Frequency offsets	Steps
1960-01-01	$-150\times10^{-10}$	N/A
1961-01-01	$-150\times10^{-10}$	$-0.005\kern{3 mu}\text{s}$
1961-08-01	$-150\times10^{-10}$	$+0.050\kern{3 mu}\text{s}$
1962-01-01	$-130\times10^{-10}$	$0$
1963-11-01	$-130\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1964-01-01	$-150\times10^{-10}$	$0$
1964-04-01	$-150\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1964-09-01	$-150\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1965-01-01	$-150\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1965-03-01	$-150\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1965-07-01	$-150\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1965-09-01	$-150\times10^{-10}$	$-0.100\kern{3 mu}\text{s}$
1966-01-01	$-300\times10^{-10}$	$0$
1968-02-01	$-300\times10^{-10}$	$+0.100\kern{3 mu}\text{s}$
1972-01-01	$0$	$-0.1077580\kern{3 mu}\text{s}$

Let $o$ be the offset, $s$ be the step, $\text{UTC}_-$ / $\text{UTC}_+$ be the UTC scale before / after jumping. For frequency offsets: $\text{TAI}-\text{UTC}=(\text{TAI}_0-\text{UTC}_0)-86400o\kern{3 mu}\text{s}\kern{3 mu}(\text{MJD}_\text{UTC literal}-\text{MJD}_{0\text{ UTC literal}})$ For steps: $\text{UTC}_+=\text{UTC}_-+s$ $\text{MJD}_{+\text{ UTC literal}}=\text{MJD}_{-\text{ UTC literal}}+\dfrac{s}{86400\kern{3 mu}\text{s}}$ (Here 'literal' is opposite to 'quasi'; $\text{TAI}_0$ and $\text{UTC}_0$ are obtained from the table) The table I provided earlier can be generated by using these.

[mrhso]

I'm aware the D2DTF comments say 9 should be safe, but evidently rounding errors have started to creep in at that point on the platforms we are both using.

Compare the tables of 1e-7 one and exact one. The former will be $4.2131700\kern{3 mu}\text{s}-4.21316999676726\kern{3 mu}\text{s}=3.23274\times10^{-9}\kern{3 mu}\text{s}$ later than the latter. $9.999999999+3.23274\times10^{-9}={\color{red}10.000000002}23274$ The analysis results match the example, indicates that the rounding result of D2DTF is correct and the table is not accurate.

[PTW19]

Obtaining a complete table is possible.

SOFA has tended to rely on the USNO Time Service, especially as the Board included members from that department. Here's the table they publish:

[]https://web.archive.org/web/20191019051734/http://maia.usno.navy.mil/ser7/tai-utc.dat

I wasn't conscious that it is a derived product and includes redundancy, but it doesn't surprise me. It's probably too late for SOFA to go back to the "urtext" of the WWV time and frequency steps, but there's no reason for ERFA not to do so. Indeed, the special licensing conditions for the SOFA DAT routines mean you could even leave it with the "iau" prefix if you wanted to.

[PTW19]

BTW thanks for pointing out that using 9 as the D2DTF ndp argument wasn't the cause of the stray 2. Reassuring.

[mhvk]

@mrhso - thanks for that investigation! I think we should make the adjustments you suggested. Would you be able to open a PR?

@PTW19 - are you sure SOFA should not follow? It feels to me @mrhso essentially uncovered a transcription bug in SOFA (even if the bug is due to the supplement, not SOFA). From the ERFA side, I really prefer there to be no differences between SOFA and ERFA (the dat.c exception really only allows us to not force our users to update the whole library when a new leap second is announced).

[PTW19]

I take it by "transcription error" you mean that the USNO Time Service description of UTC versus TAI is not a rigorous equivalent to the IERS description. Although I'm sure the SOFA Board will be uncomfortable with changing the behavior of the DAT routine, it is worth taking it as far as producing a definitive revised algorithm. So if mrhso can supply actual C code for a rigorous iauDat I'll be happy to take it from there.

[mrhso]

the USNO Time Service description of UTC versus TAI is not a rigorous equivalent to the IERS description.

That's not the case. IERS EOP also provides UTC-TAI.history, which has no difference with USNO's table except for format. In fact, both forms were produced by BIH (Bureau International de l'Heure). BIH Annual Report for 1972 clearly shows this. From this perspective, the 1e-7 accuracy is entirely due to historical limitations.

[PTW19]

From this perspective, the 1e-7 accuracy is entirely due to historical limitations.

Good. So SOFA has the choice of improving on what the BIH published in 1972 or taking it as canonical. I vote for the latter.