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 for s in xy], [s 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)?