utc_strftime("%j") fails when presented with a time near the end of the day

[bismurphy]

Hi,

I'm trying to convert Skyfield dates into the format used for TLE epochs (last two digits of year, followed by a float representing the day of the year). I'm doing this because I want to determine a spacecraft's position throughout the past year, and I have all the TLE's issued for that spacecraft in the past year, so to get the position I first load up whichever TLE was the most current at that time.

I ran into what seems to be a bug - the utc_strftime("%j") is meant to return the day-of-year of a date. 1 for January 1, 31 for January 31, 32 for February 1, etc. If this is passed a Skyfield object which is very close to the end of a day, it will increment to the next day. The following short script demonstrates this:

from skyfield.api import load
ts = load.timescale(builtin=True)
print(ts.utc(2021,1,1,23,59,59.9999798834251798496).utc_strftime("%j"))
print(ts.utc(2021,1,1,23,59,59.9999798834251798497).utc_strftime("%j"))

These two times, both of which are very near to the end of January 1, give different results. The first correctly indicates that it is still January 1, but the second, which is 10^-19 seconds later, thinks the time has shifted to January 2.

I suspect this is a matter of running into the limits of floating point precision. It may be impossible to fix this within Skyfield, but in case it's fixable, I thought I'd post this as a notice of unintended behavior. Thanks!

In my original use case, where I want the day of the year and a fractional date (how far through that day time has progressed), this becomes especially frustrating, because my fractional date value, which is found by (3600 * hour + 60 * minute + second) / 86400 returns 1 - that is, it correctly finds that we are at the end of the day. Therefore when I take the day, which has rolled over to 2, and add the fractional day, I find that it's the start of January 3. If it's not possible to improve the strftime performance, it would be cool if ts.utc(2021,1,1,23,59,59.9999798834251798497).utc would roll over together with it and return 2021,1,2,0,0,0. That way, even though it's jumping the gun to January 2, at least it correctly says it's the very beginning of January 2 rather than the very end.

[bismurphy]

Here's another example of the most troublesome behavior associated with this:

from skyfield.api import load
ts = load.timescale(builtin=True)
for t in ts.utc(2020,7,13,13,range(2050,2150,5)):
    integer_day_of_year = int(t.utc_strftime("%j"))
    hour,minute,second = t.utc[-3:]
    day_frac = (second + 60 * (minute + 60 * (hour))) / 86400
    total_day = integer_day_of_year + day_frac
    print(total_day)

This marches through the end of day 196, and then when crossing midnight, thinks it's on day 198. An excerpt of the output is here:

196.99305555555554
196.99652777777777
198.0
197.00347222222223
197.00694444444446

[brandon-rhodes]

Even though I’ve never mixed them in my own code, I think it’s a reasonable expectation that .utc_strftime('%j') and .utc switch dates at the same time instead of one rounding up before the other one does. I’ll dive into the code and see if there’s a reasonable fix waiting to be discovered!

[brandon-rhodes]

@tj-murphy — Oh, and, you might try combining .toordinal() with divmod(ordinal, 1.0) combined with subtracting out the date of January 1 (or December 31?) to go ahead and be able to move forward on your program.

[bismurphy]

Sounds great, thanks for your support!

[brandon-rhodes]

All right, I tried moving some code around so that the two routines t.utc and t.utc_strftime() are now backed by the same logic and, thus, by the same rounding logic. See whether %j now gives an answer that jumps to the next day at the same moment t.utc does! You can try the development version with:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

[bismurphy]

This looks great!

Thanks so much for the quick fix!