Unexpected result from d2dtf() very near day boundary in UTC about erfa

[mhvk]

Just to close the loop here - after puzzlement on all sides, the conclusion is that this is indeed a SOFA bug, which will be fixed in the next SOFA release (possibly as a bug-fix release, TBD). It is nice to see how bugs found via astropy continue to help upstream too! Good that you realized things were unexpected, @taldcroft!

[mhvk]

This is fixed in SOFA 17a, released today.

[mhvk]

I think pre-1972 and per leap seconds, the clocks were allowed to drift; see

erfa/src/dat.c

Lines 239 to 240 in 55fbff2

	/* If pre-1972, adjust for drift. */
	if (i < NERA1) da += (djm + fd - drift[i][0]) * drift[i][1];

I think given that drift, it would be unavoidable to sometimes get 60 seconds (but it would not be a full second).

[taldcroft]

I thought that while the UTC scale was allowed to drift that the UTC representation of the time was, by definition, always having seconds < 60. In other words there was not yet this idea of letting there be more than 60 UTC seconds in a minute.

But in truth I'm not an expert here, especially on the historical side. Here is a different demonstration of what seems unexpected to me:

In [22]: t                                                                                                                                         
Out[22]: <Time object: scale='utc' format='iso' value=1960-01-03 00:00:00.000>

In [23]: erfa.d2dtf('UTC', 6,  t.jd1, t.jd2)                                                                                                       
Out[23]: (1960, 1, 3, (0, 0, 0, 0))

In [24]: erfa.d2dtf('UTC', 6,  t.jd1, t.jd2 - 20e-12/86400)                                                                                         
Out[24]: (1960, 1, 2, (23, 59, 60, 0))

In [25]: t.jd1                                                                                                                                     
Out[25]: 2436936.0

In [26]: t.jd2                                                                                                                                     
Out[26]: 0.5

Clearly when you subtract 20 picosec from the midnight, the year/mon/day hr/min/sec representation to the nearest microsec should not change at all. Instead it is rounding 59.999999999 to 60.0 but not recognizing that this should then wrap to 1960-01-03 00:00:00.000000.

[timj]

Doesn't UTC require that a minute can have between 59 and 61 seconds in it? So if there is a leap second added that does mean that the seconds will indicate 60.

See for example https://www.iausofa.org/2021_0125_C/sofa/sofa_ts_c.pdf page 14.

[mhvk]

Hmm, I think you are right that there is a bug: the code dealing with times beyond 24 hours considers any jump a leap second, however small:

erfa/src/d2dtf.c

Lines 121 to 129 in 55fbff2

	/* Has the (rounded) time gone past 24h? */
	if ( ihmsf1[0] > 23 ) {
	/* Yes. We probably need tomorrow's calendar date. */
	js = eraJd2cal(a1+1.5, b1-fd, &iy2, &im2, &id2, &w);
	if ( js ) return -1;
	/* Is today a leap second day? */
	if ( ! leap ) {

Here leap is a boolean for dleap != 0. My sense is that the comment "is this a leap-second day" tells that the intention is to deal with regular (mini) leap seconds, not ones as small as dleap=1.11022e-16, which is what I find by inserting a print statement. Especially since this equals 2**-53, suspiciously close to what you get for a rounding error for subtracting 0.5 from something that is just one bit off. (np.finfo(np.float).eps is 2.220446049250313e-16).

[mhvk]

I've sent an e-mail to Patrick Wallace to ask about this...

[taldcroft]

Great sleuthing @mhvk!

[mhvk]

fixed by the update to erfa 1.7.3