<p>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:</p>
<pre><code>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=&quot;blue&quot;)
plt.show()
</code></pre>
<p>Then I get a plot like this:</p>
<p><img src="https://i.sstatic.net/hwiWG.png" alt="Plot generated by script" /></p>
<p>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.</p>

CC BY-SA 4.0
2022-12-21T11:35:45.590Z
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.
Sten L
CC BY-SA 4.0
2023-01-20T08:26:09.150Z
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.
Imyaf
CC BY-SA 4.0
2022-12-21T00:59:48.610Z
2022-12-21T00:59:48.610Z
Mike G
<p>I'm trying to create a map of the northern sky showing Earth's axial precession, like this wikipedia <a href="https://commons.wikimedia.org/wiki/File:Precession_N.gif" rel="nofollow noreferrer">illustration</a> by Tauʻolunga:</p>
<p><a href="https://i.sstatic.net/BvYZk.gif" rel="nofollow noreferrer"><img src="https://i.sstatic.net/BvYZk.gif" alt="precession of the equinox" /></a></p>
<p>I have figured out how to plot the stars using <a href="https://rhodesmill.org/skyfield/" rel="nofollow noreferrer">Skyfield</a> but I'm stumped when it comes to the axial precession. What I tried:</p>
<pre><code>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=&quot;blue&quot;)
</code></pre>
<p>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:</p>
<p><a href="https://i.sstatic.net/hqXjh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/hqXjh.png" alt="enter image description here" /></a></p>
<p>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)?</p>

CC BY-SA 4.0
2022-12-20T13:02:16.970Z
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.
notovny
CC BY-SA 4.0
2022-12-20T17:30:19.807Z
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)?
Sten L
CC BY-SA 4.0
2023-01-20T06:58:26.173Z
The wikipedia illustration by Tauʻolunga is incorrect. The correct wikipedia illustration is
  https://ru.wikipedia.org/wiki/Файл:Прецессия_северного_полюса_Земли.png
Imyaf
CC BY-SA 4.0
2023-09-22T18:11:06.640Z
In reference to the illustration shared by @Imyaf, here's the English [version](https://commons.wikimedia.org/wiki/File:N_pole_precession.png).
colossatr0n
<precession><skyfield>
How to plot the precession of the equinoxes using Skyfield?
Sten L
# How to plot the precession of the equinoxes using Skyfield?

(Asked by: username_0 on 2022-12-20)

I'm trying to create a map of the northern sky showing Earth's axial precession, like this wikipedia [illustration](https://commons.wikimedia.org/wiki/File:Precession_N.gif) by Tauʻolunga:

[![precession of the equinox](https://i.sstatic.net/BvYZk.gif)](https://i.sstatic.net/BvYZk.gif)

I have figured out how to plot the stars using [Skyfield](https://rhodesmill.org/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](https://i.sstatic.net/hqXjh.png)](https://i.sstatic.net/hqXjh.png)

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

**Tags:** precession, skyfield

> **Comments:**
> **username_1**: 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.  
> **username_0**: 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)?  
> **username_2**: The wikipedia illustration by Tauʻolunga is incorrect. The correct wikipedia illustration is   https://ru.wikipedia.org/wiki/Файл:Прецессия_северного_полюса_Земли.png  
> **username_3**: In reference to the illustration shared by @username_2, here's the English [version](https://commons.wikimedia.org/wiki/File:N_pole_precession.png).  

---

## 1 answers

### Answer by username_4 on 2022-12-21. Score: 1 | Accepted answer

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](https://i.sstatic.net/hwiWG.png)

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.

> **Comments:**
> **username_0**: 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.  
> **username_2**: 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.  

---

