Perhaps earthlib's earth_rotation_angle() is wrongly called at siderial_time()

[trufanov-nok]

I was comparing SOFA and skyfield's implementation of GSMT. And found some discrepency.

Let's take a look in SOFA/ERFA's code: https://github.com/liberfa/erfa/blob/master/src/era00.c#L75 theta is basically f + 0.7790572732640 + 0.00273781191135448 * t where f is a fractional part that's calculated as f = fmod(d1, 1.0) + fmod(d2, 1.0) and I would assume that t is a non-fractional part, but in SOFA/ERFA it's a whole date: t = d1 + (d2- ERFA_DJ00);.
Thus if you're passing UT1 = 2457506.7893796771204674 you end up with t = 5961.789379683789 and f = 0.789379683789 no matter how you split UT1 between dj1/dj2 arguments.

On the other side in skyfield python's package the same code is located here: https://github.com/skyfielders/python-skyfield/blob/master/skyfield/earthlib.py#L129

The function is called from gmst->siderial_time method here: https://github.com/skyfielders/python-skyfield/blob/master/skyfield/earthlib.py#L109 as theta = earth_rotation_angle(t.whole, t.ut1_fraction) for t = 2457506.7893796771204674 the t.whole will be 2457506.0 and t.fraction will be 0.789379683789.

So SOFA is using time t=T and f=T_fract, while skyfield is using t=T_whole and f=T_fract.
Firstly I thought it's a SOFA bug. But after digging into several docs I found out some equations. For ex:
https://www.sai.msu.ru/neb/rw/IERS-2010_tn36.pdf

Notice that if the equation doesn't contain a fraction part there is 1.00273781191135448 coefficient instead of 0.00273781191135448. So the fractional part somehow apperas as an attempt to exclude a whole part of this coefficient from this equation providing that later it shall be mod'ed by 2*pi. (I don't really tried to understand the math behind this change). But they don't redefine Tu - it's still a date with a fractional part.
So I think SOFA approach is correct and in skyfield at https://github.com/skyfielders/python-skyfield/blob/master/skyfield/earthlib.py#L112C31-L112C31 instead of
theta = earth_rotation_angle(t.whole, t.ut1_fraction)
should be
theta = earth_rotation_angle(t, t.ut1_fraction)

P.S. Will doublecheck this later

[brandon-rhodes]

Since Skyfield currently agrees with standard astronomy tools to very high precision, I'm going to close this for now, pending your comment:

P.S. Will doublecheck this later

Please feel free to comment further, when you do get a chance to doublecheck. And please re-open if you find an instance where Skyfield returns an incorrect angle. Thanks!

[trufanov-nok]

I doubleckecked it, ensured that it must be t, not t.whole. But didn't reported as it won't change the result too much anyway and now I almost forget the background. I'm pretty sure that it must be t, not t.whole.
You may search for 0.00273781191135448 constant in github - it'll spit out 150 pages. https://github.com/search?q=0.00273781191135448&type=code Most of them are clearly use t with fractional part. For example very first result, Cesium project:

let era = 0.779057273264 + fractionOfDay + 0.00273781191135448 * (daysSinceJ2000 + fractionOfDay); https://github.com/CesiumGS/cesium/blob/main/packages/engine/Source/Core/Transforms.js#L839

[brandon-rhodes]

You may search for 0.00273781191135448 constant in github - it'll spit out 150 pages.

Oh, wow, yes, look at all the results! The GitHub search feature does make possible some interesting archaeology involving who is using constants from which papers (or are copying them right from other libraries). 🙂

Most of them are clearly use t with fractional part. For example very first result, Cesium project: let era = 0.779057273264 + fractionOfDay + 0.00273781191135448 * (daysSinceJ2000 + fractionOfDay);

That's also Skyfield's formula, right?

https://github.com/skyfielders/python-skyfield/blob/1ac697dca28c0feb23b9947360d827d9e191caeb/skyfield/earthlib.py#L136-L137

Where jd_ut1 - T0 is equivalent to daysSinceJ2000, and fractionOfDay is the same as fraction_ut1? And then the following line adds in fractionOfDay to get the same overall sum.

[trufanov-nok]

indeed, I was most probably wrong.

[brandon-rhodes]

I probably changed things around a bit opaquely, trying to squeeze a little more resolution out. I remember I ever had a script at one time that tested whether small differences in time made it through to outputs like ERA; I should probably revive it and keep it working, so that we could see how tweaks to the formulae might let more time bits-of-precision make it through to the resulting angles.