Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Computing DUT1 #176

Closed
corrado9999 opened this issue May 4, 2018 · 11 comments
Closed

Computing DUT1 #176

corrado9999 opened this issue May 4, 2018 · 11 comments

Comments

@corrado9999
Copy link

According to IERS and Wikipedia, DUT1 is the difference between UT1 and UTC, and it is kept in the range (-1,1) seconds by using the leap seconds. I looked both in the documentation and in the code, but I could not find any reference to it. So I tried to calculate it by using a Time instance.

The only way I found was by using a "private" method called _utc_float, but I am not sure if this is correct,

>>> import skyfield.api
>>> import skyfield.timelib
>>> 
>>> ts = skyfield.api.load.timescale()
>>> t = ts.utc(2018, 3, 15)
>>> dut1 = (t.ut1 - t._utc_float()) * skyfield.timelib.DAY_S
>>> print(dut1)
0.0723391771317

The problem is that the IERS bulletin states that, starting from that very same date, the DUT1 value to be disseminated should be 0.1.

@ghost
Copy link

ghost commented May 5, 2018

I think they're rounding to the nearest tenth of a second. Remember, they say "until further notice", so a more precise value would go out of date faster. Try computing it 6 months from now or whatever, or figure out when it hits 0.15, and see if IERS releases a bulletin.

@corrado9999
Copy link
Author

corrado9999 commented May 7, 2018

I am definitely not an expert here, I was expecting to have values interpolated, but also to have the exact value at the release time. Anyway, I followed @barrycarter suggestion and made some plot.
The previous bullettin has been released on 30/11/2017 with value of 0.2. I plotted the calculated value every day since the previous bulletin until end of July. To highlight the slope, I also plotted the finite differences. Here is the result (code):

figure_1

As you can see:

  1. Slope is instable (roundoff errors?)
  2. Slope is also discontinuous; this is not surprising, I was expecting a linear interpolation, but discontinuities are neither in correspondence with bulletins, nor they are in the midpoints nor at fixed distance from when the advertised value has been reached, they simply look like happening every 3 months, starting from January.
  3. July discontinuity goes to 0, which looked pretty strange,

To better understand what happens in the future (none can predict the future, so I was expecting strange things) I did the same plot every months until mid 2019:

figure_2

So, after July discontinuity, slope changes in October (as expected) to a very low value (as to catchup with the previous line) and then goes more or less stable after January.

In the end, I was expecting a more predictable behaviour and more agreement with bulletins dates and values, but probably there are a lot of other things to take into account here (such as e.g. propagation in the future with what I think could be a model). What I would like to know is if the method to calculate DUT1 is correct or if I am misinterpreting something.

@ghost
Copy link

ghost commented May 7, 2018

Just being obnoxious, but could you correct 30/11/2018 in your comment to 30/11/2017?

@brandon-rhodes
Copy link
Member

That's a good question, @corrado9999. The value DUT1 doesn't wind up being terribly important in an application like Skyfield, because it's the difference between the very important value UT1 and the really-not-important-at-all (for astronomy calculations) UTC. While it's nice to have a way to translate to and from UTC to display clock time to the user, it's a very bad time scale for computation, because it has a 1-second discontinuity each time the standards bodies decide to sync it back up with the Earth's irregularly slowing rotation.

Instead, Skyfield needs and uses the value ΔT (delta_t in the code) which is the difference between the rock-solid, continuous, completely regular clock Terrestrial Time and the value UT1 that says where the Earth is pointed.

The series of bulletins you're reading happen to not be the source of data Skyfield uses for ΔT, which is why you're not seeing any relationship between the dates on which they issue bulletins and the inflection points in Skyfield's ΔT estimates. Instead, Skyfield uses linear extrapolation between the points in the data files:

http://maia.usno.navy.mil/ser7/deltat.data
http://maia.usno.navy.mil/ser7/deltat.preds

As you can see, these are much higher precision than the tenth-of-a-second estimates you're seeing from the DUT series of bulletins, so they're more appropriate for astronomical calculations.

I'm not sure, though, why the discontinuities you're seeing have such a clear sawtooth pattern. The precision of the floating point dates in Skyfield should be around 20μs, much smaller than it looks like the sawtooth jumps are in your graph. If you'd like to share the Python program that generates the graph, I'll be happy to look into the computations to try to figure out where the sawtooth is coming from!

@corrado9999
Copy link
Author

corrado9999 commented May 12, 2018

Probably not very clearly, by I already linked the code just before the plot. Sorry, it's very ugly and uncommented, hope it will be still useful. Although not used anywhere, do you think you could expose such value? I am not very happy to use an internal function to do this, and to compare the results with other libraries it's useful to have access to it.

@brandon-rhodes
Copy link
Member

Yes, I'd missed the link! I'll take a look, thanks.

@brandon-rhodes
Copy link
Member

I've now tried out the script! Alas, it errors out:

Traceback (most recent call last):
  File "tmp28.py", line 9, in <module>
    plt.plot(t.utc_datetime(), dut1)
NameError: name 'dut1' is not defined

And I don't see myself, even in the commented code, exactly where the dut1 value is intended to come from in the code. Could you double-check the gist and see whether a line was inadvertently omitted? Thanks!

@corrado9999
Copy link
Author

Unbelievable, I missed the most important line of code! Sorry about that, I just updated the gist.

@brandon-rhodes
Copy link
Member

I understand now that the jagged red sawtooth line belongs with the right axis, not the left axis, of the graph; the legend hadn't made it clear to me which lines corresponded to which axes. So I think all mysteries are now cleared up? To recap, repeating questions from the thread above:

The problem is that the IERS bulletin states that, starting from that very same date, the DUT1 value to be disseminated should be 0.1.

As explained, DUT1 is low precision, rounds to a tenth of a second,
and is ignored for astronomical purposes.
The value Skyfield gave you of 0.0723391771317 is simply
a more exact version of the disseminated value.

Slope is instable (roundoff errors?)

The jitter in the slope
(now that I’m looking at the correct axis!) is in the 40µs range,
which is maybe twice what I would have expected,
but reasonable for representing a Julian date as a 64-bit float.
In general, if your application requires better time accuracy
than around 0.1ms = 100µs, then Skyfield’s default choice
of a 64-bit float for time will cause jitter.

discontinuities are neither in correspondence with bulletins, nor they are in the midpoints nor at fixed distance from when the advertised value has been reached

They occur at the data points in the high precision files
deltat.data and deltat.preds.

July discontinuity goes to 0, which looked pretty strange

The deltat.preds file drops to tenth-of-a-second prediction accuracy
by mid-2018, leading to periods when linear interpolation is flat:

 2018.00      68.99               0.192         0.02
 2018.25      69.14               0.041         0.02
 2018.50      69.3                              0.2
 2018.75      69.3                              0.3
 2019.00      69.5                              0.4
 2019.25      69.6                              0.5

What I would like to know is if the method to calculate DUT1 is correct

Yes, you are predicting it correctly!

If those are indeed most of the outstanding issues,
then the remaining complaint that isn’t addressed
is that you’re having to call a hidden method
to access Skyfield’s interpolation of leap seconds.
Am I correct that if I were to pop that computation out
as a separate calculation of its own,
that you’d have everything you need for your DUT1 routines?

@brandon-rhodes
Copy link
Member

brandon-rhodes commented May 21, 2018

After looking at the various options, I decided that your question was so important that .dut1 should become an attribute of all Time objects. You should now be able to reproduce the graph at the Wikipedia page:

https://en.wikipedia.org/wiki/DUT1

Here's an example notebook that does so:

https://github.com/skyfielders/astronomy-notebooks/blob/master/Skyfield-Notes/DUT1.ipynb

by plotting t.dut1. Enjoy! :)

@corrado9999
Copy link
Author

Great! Thank you for your support.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants