3
$\begingroup$

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

precession of the equinox

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:

enter image description here

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)?

$\endgroup$
4
  • $\begingroup$ You're looking at the North Celestial Pole relative to the position of Stockholm, if i'm not mistaken. That's why you're not seeing any movement outside of floating point precision, as both positions are stapled to the current position and orientation of Earth at the time of observation. $\endgroup$
    – notovny
    Dec 20, 2022 at 13:02
  • $\begingroup$ Thanks. I was hoping that stockholm would be the current position of Stockholm, since I provided a date when I made it. This is probably what I don't understand. How can I make the 'projection' correspond Stockholm fixed at today's date, and then project the position of the celestial axis at a date in the past (or future)? $\endgroup$
    – Sten L
    Dec 20, 2022 at 17:30
  • 1
    $\begingroup$ The wikipedia illustration by Tauʻolunga is incorrect. The correct wikipedia illustration is ru.wikipedia.org/wiki/Файл:Прецессия_северного_полюса_Земли.png $\endgroup$
    – Imyaf
    Jan 20, 2023 at 6:58
  • 1
    $\begingroup$ In reference to the illustration shared by @Imyaf, here's the English version. $\endgroup$ Sep 22, 2023 at 18:11

1 Answer 1

1
$\begingroup$

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:

Plot generated by script

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.

$\endgroup$
2
  • $\begingroup$ Perfect, thank you! I still want it shifted to the observer's location (so that it matches the stars I plotted), but I can figure that out probably. $\endgroup$
    – Sten L
    Dec 21, 2022 at 11:35
  • 1
    $\begingroup$ Discontinuity is must be. It is not error. The pole does not depict the closed circle. At long term the inclination of the Earth's axis is changed about 8 degrees relative to the stars. Because the Earth’s orbit is also precessing about 4 degrees. Also the rotation velocity is changed in 1.5 times. $\endgroup$
    – Imyaf
    Jan 20, 2023 at 8:26

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .