%pylab inline
import ephem
from mpl_toolkits.basemap import Basemap
Populating the interactive namespace from numpy and matplotlib
This notebook demonstrates how the date and epoch arguments
to PyEphem’s compute()
call change the appearance of a constellation
over very long periods of thousands of years.
It was inspired by the question posed by this GitHub issue:
https://github.com/brandon-rhodes/pyephem/issues/61
The Big Dipper is made up of eight principle stars:
names = ['Dubhe', 'Merak', 'Phecda', 'Megrez',
'Alioth', 'Mizar', 'Alcor', 'Alcaid',]
stars = [ephem.star(name) for name in names]
Plotting their position in the sky is simple thanks to the Basemap extension for Matplotlib, which is particularly easy to install if you are using the Anaconda distribution:
$ conda install basemap
Computing where the stars are located at any particular moment and then displaying them can be described in two simple functions, the first dealing only with numbers and the second plugging those numbers into a plot:
def dipper_radec(*args):
for star in stars:
star.compute(*args)
ra = [-degrees(star.a_ra) for star in stars]
dec = [degrees(star.a_dec) for star in stars]
return ra, dec
def plot_dipper(*args):
ra, dec = dipper_radec(*args)
width = 4000000
m = Basemap(width=width, height=width, projection='aeqd',
lat_0=dec[3], lon_0=ra[3])
m.drawparallels(np.arange(-80,81,10))
m.drawmeridians(np.arange(-180,180,10))
x, y = m(ra, dec)
m.plot(x, y, marker='o', linestyle='--', color='b')
plt.title('date={} epoch={}'.format(*args))
The epoch makes no difference to the shape of the constellation.
Indeed, the epoch also has no effect on almost any of a star’s
properties or coordinates — both ra
and dec
as well as alt
and az
are completely unaffected!
The only coordinates that pay attention to the epoch
are the astronometric a_ra
and a_dec
because they, and they alone, express coordinates
in a fixed coordinate system tied to a specific date.
That is why the code above is careful to use a_ra
and a_dec
—
because otherwise it would not pay attention to epoch at all.
The reason that the epoch makes no difference to the shape of the Big Dipper is that the epoch merely choses how we name different positions in the sky. It chooses one naming scheme over another. This is clear if we plot the Big Dipper’s current coordinates using the sky coordinate systems of three different dates: 8,000 years ago, today (well, the year 2000), and 8,000 years in the future.
While the celestial north and south pole were at a different location on each of these three dates, giving all of the points in the sky a different right ascension and declination, that merely tips the Big Dipper in different directions as we name its position according to these different celestial poles:
plot_dipper('2000', '-6000')
show()
plot_dipper('2000', '2000')
show()
plot_dipper('2000', '10000')
The date — the first argument to compute()
—
is what can change the shape of a constellation,
because the stars of our galaxy gradually change
their relative orientation to each other
as we all swing slowly around the galactic center
in our different orbits.
When checking out the way a constellation shape changes, there are two approaches.
The first approach is to use the modern J2000 coordinate system for each of the plots. This means that the “right side up” direction of all of the plots stays the same, so all of the motion between one plot to the next is really truly the motion of the stars themselves.
To see much motion, we will need to go much farther back in time and ahead in time than we did before. We can try 50,000 years in each direction:
plot_dipper('-50000', '2000')
show()
plot_dipper('2000', '2000')
show()
plot_dipper('50000', '2000')
But if we want to imagine how the constellation looked in the sky of its time, we probably want to switch to the coordinate system of that distant (or far-future) era. (If, indeed, the constellation would still have been recognized as such — eventually the motion will carry these stars into groupings with stars that today belong to other constellations!)
Switching to the coordinate system of the other era will rotate the Dipper left or right in our diagram, making it more difficult to see which stars have moved and in which direction.
But the other era’s coordinate system will give us an idea of whether the Dipper’s stars were up closer to the pole, or down closer to the equator, in the sky in which they appeared.
If only PyEphem could compute in far-past and far-future coordinate systems correctly!
If we try computing the location of the Big Dipper in the skies of the far past or future, we get nonsense: 50,000 years ago PyEphem thinks the Dipper was in the southern sky (!), and 50,000 years in the future it places it up at the north celestial pole:
plot_dipper('-50000', '-50000')
show()
plot_dipper('2000', '2000')
show()
plot_dipper('50000', '50000')
What could be causing such chaos?
Let us plot the position of Polaris, which is currently quite near to the north celestial pole. Its path over 26,000 years should be the smooth circle caused by the precession of the Earth’s axis:
polaris = ephem.star('Polaris')
width = 9000000
m = Basemap(width=width, height=width, projection='aeqd',
lat_0=70.0, lon_0=280.0)
m.drawparallels(np.arange(-80,81,10))
m.drawmeridians(np.arange(-180,180,10))
for year in range(-13001, 13000, 500):
polaris.compute(str(year), epoch=str(year))
x, y = m(degrees(polaris.a_ra), degrees(polaris.a_dec))
m.plot(x, y, marker='o', linestyle='--', color='b')
So far, so good!
The path of Polaris in the above diagram is exactly what we would expect from precession, and matches other diagrams that show the phenomenon.
But everything goes haywire if we step out to many tens of thousands of years in the past and future:
width = 18000000
m = Basemap(width=width, height=width, projection='aeqd',
lat_0=70.0, lon_0=280.0)
m.drawparallels(np.arange(-80,81,10))
m.drawmeridians(np.arange(-180,180,10))
for year in range(-50001, 50000, 500):
polaris.compute(str(year), epoch=str(year))
x, y = m(degrees(polaris.a_ra), degrees(polaris.a_dec))
m.plot(x, y, marker='o', linestyle='--', color='b')
Even though the calculations inside PyEphem seems to be stable over what look like two complete precessions of the Earth’s axis, outside of that safe range it goes crazy and stops returning stable results.
Well, drat.
Which libastro
routine is at fault?
This was going to be such a simple notebook, involving only Python code. But now it appears we have a bug on our hands, and need to dive down into the lower layers to see whether it is the precession or nutation algorithm that is dying.
Thank goodness for ctypes
!
Here goes:
from ctypes import *
cdll.LoadLibrary('../ephem/_libastro.so')
lib = CDLL('../ephem/_libastro.so', use_errno=True)
lib.precess.argtypes = [
c_double, c_double, POINTER(c_double), POINTER(c_double)
]
ra = []
dec = []
ra_double = c_double()
dec_double = c_double()
epoch1 = ephem.J2000
year = arange(-50001.0, 50000.0, 100)
for y in year:
epoch2 = epoch1 + 365.25 * y
ra_double.value = 0.0
dec_double.value = 0.0
lib.precess(epoch1, epoch2,
byref(ra_double), byref(dec_double))
ra.append(ra_double.value)
dec.append(dec_double.value)
plot(year, ra)
show()
plot(year, dec)
[<matplotlib.lines.Line2D at 0xb05eeccc>]
lib.nutation.argtypes = [
c_double, POINTER(c_double), POINTER(c_double)
]
deps = []
dpsi = []
deps_double = c_double()
dpsi_double = c_double()
epoch1 = ephem.J2000
year = arange(-50001.0, 50000.0, 100)
for y in year:
epoch2 = epoch1 + 365.25 * y
deps_double.value = 0.0
dpsi_double.value = 0.0
lib.nutation(epoch2, byref(deps_double), byref(dpsi_double))
deps.append(deps_double.value)
dpsi.append(dpsi_double.value)
plot(year, deps)
show()
plot(year, dpsi)
[<matplotlib.lines.Line2D at 0xb056154c>]
And that settles it!
Nutation looks quite stable over long periods: it stays within its own defined little range of tiny numbers, without going wild.
But precession, once it gets out past 10,000 years and 20,000 years ahead or behind the present, goes completely wild and cannot be trusted.
I supposed I should add some sort of guard statements to PyEphem so it will not try computing out that far, so that other users do not run into problems.
What about PyEphem’s replacement, Skyfield?
from skyfield.api import JulianDate
from skyfield import precessionlib
years = arange(-50000, 50000)
jd = JulianDate(tai=(years,))
plot(years, precessionlib.compute_precession(jd.tdb)[0][0])
[<matplotlib.lines.Line2D at 0xb07643cc>]
Wow!
Skyfield, I believe, uses very modern precession formulae derived from more recent IAU resolutions than those that power PyEphem. Yet we see the same limitation: once we are out past 20,000 years before or after the present, the formula goes unstable and ceases to predict accurately.
So maybe Skyfield needs a similar sanity check to protect users who might want to look too far out!
I will update the issue with a link to this notebook.