Skip to content

Commit d618ee9

Browse files
brandon-rhodesbrandon-rhodes
committed
Fix precision when UTC is given as fractional days
When days with attached fractions were given to UTC contructors, like `ts.utc(2023, 1, 76.01808579)`, the day fraction could lose precision, because the contructor would build the full Julian day number before rescuing the fraction. Now, the full fraction is preserved. This makes the `.epoch` more exact when an `EarthSatellite` is loaded from a TLE file, since it happens that TLE files specify the epoch using a day-of-year with a decimal fraction attached. Fixes #1141.
1 parent a1e6e30 commit d618ee9

5 files changed

Lines changed: 87 additions & 21 deletions

File tree

‎diary.txt‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
2+
## 2026 September 20 Sunday
3+
4+
I am nearly done improving `ts.utc()` so that any fraction attached to
5+
its `day` argument is handled with the same caution as fractions
6+
attached to its `second`, which should fix this issue:
7+
8+
https://github.com/skyfielders/python-skyfield/issues/1141
9+
10+
Having nearly worked out the kinks, I was surprised to find that I’d
11+
broken the test `test_velocity_in_ITRF_to_GCRS2()`. With the old logic,
12+
it receives a time with these attributes, and succeeds:
13+
14+
.whole [2459047. 2459047.]
15+
.tai_fraction [0.8691782407407408, 0.8691898148148148]
16+
.tt_fraction [0.8695507407407408, 0.8695623148148148]
17+
18+
But the new code produces these values and makes it fail:
19+
20+
.whole [2459047. 2459047.]
21+
.tai_fraction [0.8691782407407407, 0.8691898148148148]
22+
.tt_fraction [0.8695507407407407, 0.8695623148148148]
23+
24+
But does this really mean that the new value is worse? The test is a
25+
finicky one, testing whether the velocity matrix returned from the
26+
no-longer-used routine `ITRF_to_GCRS2()` really predicts how far a
27+
position will spin around the frame’s axis in 1 second. Does the test’s
28+
failure really mean that the newly rewritten `Timescale._utc()` logic
29+
needs to be tightened up? Or was it mere coincidence that the old logic
30+
produced a final decimal place that made the test slightly tighter than
31+
otherwise?
32+
33+
Let’s instead judge the new decimal on its own merits. The time 8 AM is
34+
actually 20h because Julian dates start at noon; and 37 leap seconds had
35+
accumulated by 2020. So we can ask, what fraction of a day would the
36+
time 20:51:37 AM ideally produce? Let’s use `dc` to compute it without
37+
the limited precision of machine floating point:
38+
39+
99k 20 24 / p 51 1440 / p 37 86400 / p + + p
40+
41+
The high-precision result is:
42+
43+
.86917824074074074074074074074074074074074074074074074074074074074074\
44+
0740740740740740740740740740739
45+
46+
The `407` looks to be a true repeating decimal. So the old logic that
47+
produced an `8` at the end was flawed, because the final digit should
48+
have rounded down, not up; and so it’s just a coincidence that the old
49+
logic happened to make this test happier. I will relax its bounds.

‎skyfield/tests/test_earth_satellites.py‎

Lines changed: 2 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -110,17 +110,12 @@ def test_appendix_c_satellite():
110110

111111
rTEME, vTEME, error = sat._position_and_velocity_TEME_km(t)
112112

113-
# TODO: This used to be accurate to within 1e-8 but lost precision
114-
# with the move to SGP4 2.0. Is the difference an underlying change
115-
# in the algorithm and its results? Or something else?
116-
epsilon = 1e-4
113+
epsilon = 1e-8
117114
assert abs(-9060.47373569 - rTEME[0]) < epsilon
118115
assert abs(4658.70952502 - rTEME[1]) < epsilon
119116
assert abs(813.68673153 - rTEME[2]) < epsilon
120117

121-
# TODO: Similar to the above, this used to be 1e-9. Then the Time
122-
# object started storing UTC as seconds, and it got worse.
123-
epsilon = 5e-8
118+
epsilon = 1e-9
124119
assert abs(-2.232832783 - vTEME[0]) < epsilon
125120
assert abs(-4.110453490 - vTEME[1]) < epsilon
126121
assert abs(-3.157345433 - vTEME[2]) < epsilon

‎skyfield/tests/test_positions.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -207,7 +207,7 @@ def test_velocity_in_ITRF_to_GCRS2():
207207
relative_error = (length_of(actual_motion - predicted_motion)
208208
/ length_of(actual_motion))
209209

210-
acceptable_error = 1e-11
210+
acceptable_error = 1.8e-11
211211
assert relative_error < acceptable_error
212212

213213
def test_light_time_method():

‎skyfield/tests/test_timelib.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -197,6 +197,9 @@ def test_strftime_with_microseconds():
197197
assert t.tt_strftime('%Y %S %f') == ['1980 00 000000'] * 2
198198
assert t.tdb_strftime('%Y %S %f') == ['1980 59 998471'] * 2
199199

200+
t = ts.utc(2023, 8, 10, 6, 21, 45.9)
201+
assert t.utc_strftime('%Y %S %f') == '2023 45 900000'
202+
200203
def test_tai_fraction_loses_no_precision(ts):
201204
t = ts.tai_jd(2459008.0, 0.0123456789)
202205
assert t.whole == 2459008.0

‎skyfield/timelib.py‎

Lines changed: 32 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
from collections import namedtuple
66
from datetime import date, datetime, timedelta
77
from numpy import (
8-
array, concatenate, cos, float64, int64, isnan, isinf, linspace,
8+
array, concatenate, cos, float64, floor, int64, isnan, isinf, linspace,
99
nan, ndarray, nonzero, pi, rollaxis, searchsorted, sin, where, zeros_like,
1010

1111
)
@@ -181,25 +181,44 @@ def utc(self, year, month=1, day=1, hour=0, minute=0, second=0.0):
181181
def _utc(self, tup):
182182
# Build a Time from a UTC tuple, carefully preserving its exact
183183
# second number in the Time's hidden TAI seconds field.
184+
184185
year, month, day, hour, minute, second = tup
185186
cutoff = self.julian_calendar_cutoff
186187

187-
# Figure out exactly the TAI second number.
188-
seconds = (julian_day(year, month, day, cutoff) - 0.5) * DAY_S
189-
seconds, sfr = divmod(seconds, 1.0) # in case there were any fractions
190-
seconds += interp(seconds, self._leap_utc, self._leap_offsets)
191-
more = hour * 3600.0 + minute * 60.0 + second
192-
seconds2, sfr = divmod(sfr + more, 1.0)
193-
seconds += seconds2
188+
# Carefully preserve the day `fraction` before losing precision
189+
# by building the Julian day number.
190+
whole = floor(day)
191+
fraction = day - whole - 0.5
192+
whole = julian_day(year, month, whole, cutoff)
193+
194+
# Use "Julian seconds" to index into the leap second table.
195+
jseconds = whole * DAY_S
196+
leap_seconds = interp(jseconds, self._leap_utc, self._leap_offsets)
197+
198+
# For accurate reconstruction of UTC later, save the time as TAI
199+
# seconds plus the exact seconds fraction, in a secret tuple.
200+
sfloor = floor(second)
201+
sfraction = second - sfloor
202+
seconds = sfloor + leap_seconds + minute * 60.0 + hour * 3600.0
203+
seconds2 = seconds + fraction * DAY_S
204+
205+
sfloor2 = floor(seconds2) # in case minutes or hours had fractions
206+
sfraction2 = seconds2 - sfloor2
207+
tai_seconds = (
208+
sfloor2 + whole * DAY_S,
209+
sfraction + sfraction2,
210+
)
194211

195-
# For the other timescales, use the usual Julian date + fraction.
196-
whole, fraction = divmod(seconds, DAY_S)
197-
fraction += sfr
198-
fraction /= DAY_S
212+
# And save the time more conventionally, counted in days.
213+
fraction += (seconds + sfraction) / DAY_S
214+
whole2 = floor(fraction)
215+
whole += whole2
216+
fraction -= whole2
199217

200218
t = Time(self, whole, fraction + tt_minus_tai)
201219
t.tai_fraction = fraction
202-
t._tai_seconds = seconds, sfr
220+
t._tai_seconds = tai_seconds
221+
203222
return t
204223

205224
def _jd(self, year, month, day, hour, minute, second):

0 commit comments

Comments
 (0)