Calculation of ICRF and J2000 frames difference

[mindaugl]

Hi, firstly just wanted to thanks for all people who contributed to Skyfield development/documentation/testing/etc. - it is a great tool for getting astronomy related answers and for learning astronomy itself.

I was trying to compute the difference in arcsseconds between the ICRF and J2000 frames. In the documentation and in the resources that I have seen online, the difference between the frames is qouted as being less than 0.02 arcseconds (https://rhodesmill.org/skyfield/positions.html):

The ICRS axes are within 0.02 arcseconds of the old J2000 axes, so many scripts simply treat J2000 coordinates as modern ICRS coordinates.

I was thinking to see this difference by just calling radec() on Astrometric position with epoch date set to 2000/01/01, versus using no argument for epoch (which would correspond to using ICRF frame directly). My code is below:

from skyfield.api import load

planets = load('de421.bsp')
earth = planets["earth"]
jupiter = planets["jupiter_barycenter"]

ts = load.timescale()
j2000 = ts.tt(2000, 1, 1, 12, 0, 0)
times = [ts.tt(1900, 6, 21, 12, 0, 0), ts.tt(2000, 6, 21, 12, 0, 0), ts.tt(2050, 6, 21, 12, 0, 0)]

for time in times:
    earth_baryc = earth.at(time)
    jupiter_astro = earth_baryc.observe(jupiter)
    print(time.tt_strftime())
    print("    ICRF: ", jupiter_astro.radec())
    print("    J2000:", jupiter_astro.radec(j2000))
    ra1, dec1, _ = jupiter_astro.radec()
    ra2, dec2, _ = jupiter_astro.radec(j2000)
    print("    ra difference (arcseconds):", ra1.arcseconds() - ra2.arcseconds())
    print("    dec difference (arcseconds):", dec1.arcseconds() - dec2.arcseconds())

The (combined) difference for Jupiter (and other planets) it shows to be around 15 arcseconds though, e.g.:

2050-06-21 12:00:00 TT
ICRF: (<Angle 08h 28m 40.76s>, <Angle +19deg 38' 53.9">, <Distance 6.09988 au>)
J2000: (<Angle 08h 28m 39.72s>, <Angle +19deg 38' 52.6">, <Distance 6.09988 au>)
ra difference (arcseconds): 15.595133518567309
dec difference (arcseconds): 1.2448747462622123

Maybe somebody could comment what is missing in the script above and how to practically see the difference of the two frames? Thank you in advance!

[brandon-rhodes]

Good question! You are seeing a large difference because with .radec(j2000) you are asking Skyfield for RA and dec relative to the Earth's true equator and equinox at the moment of J2000, which includes the crazy effects of nutation. Instead, to compare with ICRS, you want to ask about the Earth's mean equator and equinox as of that moment.

And you could, using the new mean_equator_and_equinox_of_date reference frame that I just added to Skyfield, except for two problems. First, I don't seem to have fully documented how to use it. And, second, I accidentally left out the 0.02-arcsecond rotation between J2000 and the ICRS. Drat. I'll go get that fixed.

In the meantime, to see the amount of the rotation between frames, you can take a look at this matrix:

from .framelib import ICRS_to_J2000

It should account for the 0.02 arcseconds you saw mentioned in the docs.

I'll plan to comment here again once I've gotten the mean_equator_and_equinox_of_date frame fixed.

[mindaugl]

Hi Brandon, thank you for the reply! Yes, checked that ICRS_to_J2000 matrix has about 0.02 arcseconds effect on the calculated RA/DEC, and nutation matrix has the big impact when using true equator and true equinox for the j2000 epoch as specified in the .radec() call. It will be interesting to check when the changes are made to the mean_equator_and_equinox_of_date frame, thanks!