skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.38k stars 208 forks source link

Apply dUT1 to Timescales? #930

Closed wg1562 closed 6 months ago

wg1562 commented 6 months ago

I'm trying to use a UTC time to build a timescale and get GAST. I am loading the IERS finals.daily values, but I don't think the dUT1 value to the time to get the UT1 as expected.

From below code, I print out the t.dut1 value and it is not the expected value from the finals value:

url='https://datacenter.iers.org/data/latestVersion/finals.daily.iau2000.txt' 
with load.open(url) as f:
    finals_data = iers.parse_x_y_dut1_from_finals_all(f)

ts = load.timescale()
iers.install_polar_motion_table(ts, finals_data) 

t = ts.utc(2023, 11, 16, 0, 0, 0)
print('t.dut1: ', t.dut1)
print('UTC date and time:', t.utc_strftime())
hr, min, sec = DDtoDMS(t.gast)
print("GAST (H M S.S): {:.0f}:{:.0f}:{:.3f}".format(hr, min, sec))

Output:

t.dut1:  0.04811789999999405
UTC date and time: 2023-11-16 00:00:00 UTC
GAST (H M S.S): 3:39:14.631

But...the dUT1 value should be 0.0106376. I don't know where the incorrect value came from. The results is a GAST calculation that is off. I have a work around as shown below, and it gives the correct result, but it would be nice if Skyfield applied the dUT1 from the finals file. Or maybe it is, and I'm just doing something wrong?

sec_utc=0
dut1=0.0106376         
sec_ut1 = sec_utc+dut1
t = ts.ut1(2023, 11, 16, 0, 0, sec_ut1)
print('UTC date and time:', t.utc_strftime())
hr, min, sec = DDtoDMS(t.gast)
print("GAST (H M S.S): {:.0f}:{:.0f}:{:.3f}".format(hr, min, sec))

Output:

UTC date and time: 2023-11-16 00:00:00 UTC
GAST (H M S.S): 3:39:14.583

From USNO website (and calculation by me):

#USNO (UT1 2023-Nov-16  00:00:00.0): GAST (H M S.S): 03:39:14.5826
#(UTC + dUT1: GAST 03:39:14.5932)
brandon-rhodes commented 6 months ago

Could you share your DDtoDMS() routine, so that I can try your script out? Thanks!

wg1562 commented 6 months ago
def DDtoDMS(DD):
    d = np.trunc(DD)                  #degrees
    m = abs(np.trunc((DD-np.trunc(DD))*60))  #minutes 
    s = abs(((DD-np.trunc(DD))*60-np.trunc((DD-np.trunc(DD))*60)))*60    #seconds 
    return d, m, s
brandon-rhodes commented 6 months ago

~Thanks! The output that I get is:~

Edit: never mind, I wasn't scrolling far enough back to see the first script output that you shared. When I run the script I get what looks to be the same output, so I'll now study your other code and output, and put together an answer.

brandon-rhodes commented 6 months ago

The problem is that the method install_polar_motion_table(), alas, does only and literally that: reads the xp and yp values from the file you give it. It doesn't read or use any of the other columns.

The file that Skyfield knows how to get dUT1 from is finals2000A.all. Could you try this in your script?

load.download('https://maia.usno.navy.mil/ser7/finals2000A.all')
ts = load.timescale(builtin=False)

If it works, then I will at least update the docs to make clear that the polar-motion method only loads the polar motion data. Maybe I should also add a method beside it for loading a finals2000A.all file, instead of the awkward maneuver of builtin=False? Also, it looks like the URL used by builtin=False might be out of date, which is why my little bit of example code in this comment specifies a manual URL; I should fix the built-in URL.

wg1562 commented 6 months ago

Ha! That worked! I was using the finals.daily vs. finals.all, since the size is much smaller (33K vs. 3.4M).

It also looks like the 'builtin=False' parameter makes a difference in the dUT1 value pulled. Thanks for the tip!