How to plot the precession of the equinoxes using Skyfield?

Question

I'm trying to create a map of the northern sky showing Earth's axial precession, like this wikipedia illustration by Tauʻolunga:

I have figured out how to plot the stars using Skyfield but I'm stumped when it comes to the axial precession. What I tried:

import numpy as np
from matplotlib import pyplot as plt
from skyfield.projections import build_stereographic_projection
from skyfield.api import Star, load, wgs84, N, S, W, E

ts = load.timescale()
stockholm = wgs84.latlon(59.349304*N, 18.025918*E, elevation_m=28).at(ts.tt(2022, 12, 1, 12, 0))
position = stockholm.from_altaz(alt_degrees=90, az_degrees=0)
projection = build_stereographic_projection(position)

eph = load('de441.bsp')
earth = eph['earth']

celestial_pole = Star(ra_hours=(0, 0, 0), dec_degrees=(90, 0, 0))
xy = []
for year in range(1969, 17190, 100):
    pt = ts.tt(year, 9, 1, 12, 0)
    pos = earth.at(pt).observe(celestial_pole)
    xy.append(projection(pos))
plt.scatter([s[0] for s in xy], [s[1] for s in xy], s=5, color="blue")

But this results in a tiny (1e-15) mostly East-West pendulum movement whereas the whole sky ranges from -1 to 1 both N-S and E-W on the scale I used:

I'm not sure I go about the above the right way. Is there a better way to find the direction of Earth's axis at an arbitrary timepoint in the past or future (±26,000 years)?

Answer 1

I think you need the ICRS coordinates of the celestial pole of date. I would leave out the observer's location and compute a new position each time through the loop:

from matplotlib import pyplot as plt
from skyfield.projections import build_stereographic_projection
from skyfield.positionlib import position_of_radec
from skyfield.api import load

ts = load.timescale()
ecliptic_pole = position_of_radec(ra_hours=18.0, dec_degrees=66.56)
projection = build_stereographic_projection(ecliptic_pole)

xy = []
for year in range(-11000, 15000, 200):
    pt = ts.tt(year, 9, 1, 12, 0)
    celestial_pole = position_of_radec(ra_hours=0.0, dec_degrees=90.0, epoch=pt)
    xy.append(projection(celestial_pole))
plt.scatter([s[0] for s in xy], [s[1] for s in xy], s=5, color="blue")
plt.show()

Then I get a plot like this:

I think the discontinuity is due to the precession model which Skyfield implements. IAU 2000A precession is accurate between 1000 and 3000 CE but not necessarily valid outside that interval.