In the last two notebooks about proper motion$^{1,2}$ I saw that the big dipper was placed in really strange positions in the night sky in the far-away future, and in the ancient past. We resolved that the proper way to compute the positions of stars while mimicking an observer on earth was to vary both year and epoch together to get the proper celestial orientation.
I filed an issue in the PyEphem issue tracker$^3$ just to get confirmation, that varying time and epoch together was the best route. It was, but there is an intrinsic limitation of the IAU precession model that PyEphem inherits from libastro. Skyfield has a similar limitation. As Brandon Rhodes says$^4$
Far-past and far-future epochs are broken
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
This particular notebook wound up being more of personal research on the subject and less code, because it is my first time reading about precession matrices. I originally intended to utilize data from HORIZONS to correct for precession, because it uses a long-term precession model that can be used for many hundreds of thousands of years. I found out this was a misconception of mine, because the precession model is "baked into" this JPL ephemeris. I include this anyway, because I think the notes are helpful to anyone who may be in the same place as me.
The follow-up vondrak.ipynb (nbviewer) does implement a long-term precession matrix from beginning to end.
Well, I'm going to have to find a different model or formula to intercede when past the bounds of J2000 ±1000 years.
Summary of conclusions:
I looked into Barry's tip. Thanks! We're getting into stuff that is way beyond me calculation-wise, but I found a paper by Owen (1990) that is extremely helpful. More on it below, but we can indeed select time spans from 'BC 9998-03-20 to AD 9999-12-31' for the earth ephemeris using the HORIZONS interface.
HORIZONS uses two different precession models. One between 1799-Jan-1 to 2202-Jan-1, and the other for outside those bounds.
For the time-span of 1799-Jan-1 to 2202-Jan-1, the official IAU precession model [16] of Lieske is used. As published, this model is valid for only ~200 years on either side of the J2000.0 epoch. This is due to round-off error in the published coefficients and truncation to a 3rd order polynomial in the expressions for the Euler rotation angles. Therefore, outside this interval, the long-term precession and obliquity model [17] of Owen is used to maintain accuracy in the calculation of apparent ("of-date") quantities.
This model is a rigorous numerical integration of the equations of motion of the celestial pole using Kinoshita's model for the speed of luni-solar precession.
Precession (IAU) from 1799-Jan-1 to 2202-Jan-1:
[16]: Lieske, J., "Precession Matrix Based on IAU (1976) System of Astronomical Constants", Astron. Astrophys. 73, 282-284, 1979. (pdf)
Precession (long-term) before 1799-Jan-1 and after 2202-Jan-1:
[17]: Owen, William M., Jr., (JPL) A Theory of the Earth's Precession Relative to the Invariable Plane of the Solar System, Ph.D. Dissertation, University of Florida, 1990. (pdf)
Chapter 4 by Owen, pg 84, is dedicated to the long-term theory.
Pg 85, on Kinoshita's luni-solar precession used to integrate the motion of the celestial pole:
... sufficed to give numerical accuracies on the order of a nanoarcsecond after 500,000 years.
A Laskar file of eccentricity and inclination variables spanning a million years (500,000 to either side of J2000) was used for orbital elements and numerical integrations.
Pg 242:
Chebyshev polynomials to degree nine were fit over 8000-year intervals to the angles obtained from a million-year numerical integration based on the work of Laskar (1990) for the motion of the ecliptic and Kinoshita (1977) for the motion of the equator.
Fig 4-1, pg 105-110, shows the change in precession angles over these 10000 centuries.
Ch 4.7, Sources of Error in the Long-Term Theory, pg 116:
The dominant uncertainty in the long-term theory is in the initial speed of general precession in longitude. The currently-accepted value of 5029".0966/century may in error by a significant fraction of an arcsecond per century; in other words, by one part in 10⁴. And since this uncertainty is in a rate, the error in the accumulated angle pₐ will grow roughly linearly with time. After 500,000 years, an error of 0".36/century in p₁ would cause an error of half a degree in pₐ. An error of this magnitude easily swamps all other effects.
For this reason it is fruitless at the moment to extend the integration much past the million-year interval covered here.
Pg253:
Figure 5-1: J2000 Equatorial Coordinates of the Ecliptic Pole and of the Pole of the Invariable Plane. The positions of the ecliptic pole at T = -5000, T = 0, and T = +5000 are marked with an open square, open circle, and filled circle, respectively. The position of the pole of the invariable plane is marked with an asterisk.
Where T is in centuries.
The above-mentioned time spans that HORIZONS limits us to for the Earth stem from Ephemeris File DE431MX.
ID Name Available Time Span Ephemeris File
399 Earth B.C. 9998-03-20 to A.D. 9999-12-31 DE431MX
From JPL Planetary and Lunar Ephemerides:
DE431 :
Created April 2013; includes librations and 1980 nutation.
Referred to the International Celestial Reference Frame version 2.0.
Covers JED -0.3100015.5, (-13200 AUG 15) to JED 8000016.5, (17191 MAR 15).DE430 and DE431 are documented in the following document:
http://ipnpr.jpl.nasa.gov/progress_report/42-196/196C.pdf
I haven't read The Planetary and Lunar Ephemerides DE430 and DE431 (2014) as thoroughly as I did A Theory of the Earth's Precession Relative to the Invariable Plane of the Solar System. However, I think we can ignore the upper and lower limits of B.C. 9998-03-20 to A.D. 9999-12-31 without further proof that this ephemeris is tied into the precession theory of Owen's.
Of interest though is pg 10.
F. Orientation of the Earth
Only the long-term change of the Earth's orientation is modeled in the ephemeris integration. The Earth orientation model used for the DE430 and DE431 integration is based on the International Astronomical Union (IAU) 1976 precession model [21,22] with an estimated linear correction and on a modified IAU 1980 nutation model [23] including only terms with a period of 18.6 years.
[21]: J. H. Lieske, T. Lederle, W. Fricke, and B. Morando, “Expression for the Precession Quantities Based on the IAU (1976) System of Astronomical Constants,” Astronomy and Astrophysics, vol. 58,pp. 1–16, 1977.
[22]: J. H. Lieske, “Precession Matrix Based on IAU (1976) System of Astronomical Con- stants,” Astronomy and Astrophysics, vol. 73, pp. 282–284, 1979.
[23]: P. K. Seidelmann, “1980 IAU Theory of Nutation: The Final Report of the IAU Working Group on Nutation,” Celestial Mechanics, vol. 27, pp. 79–106, 1982.
All of these pre-date Owen's theory. At the time of writing, I am unsure whether I can garner information on purely the precession model by using HORIZONS. I corresponded with Jon D. Giorgini, who cleared up my misconceptions on how DE431, the short-term, and long-term models were implemented by HORIZONS.
the long-term precession model is used only to compute a rotation matrix to find apparent positions with respect to the "true-of-date" coordinate system. It does not affect the dynamical model (an exception being any asteroid or comet that comes within 0.01 au of Earth, for which Earth oblateness acceleration is included in the numerical integration, requiring knowledge of pole direction). The long-term precession model is not available separately in the sense of being able to get a table of Euler angles from Horizons as a function of time.
I don't think the bounds of B.C. 9998-03-20 to A.D. 9999-12-31 are limiting. They are due to the "baked-in" aspects of the IAU76 precession model with DE431.
Note that Horizons is not integrating planetary bodies (only asteroids and comets); it interpolates (reads) planetary positions and velocities from DE431, which was developed and integrated by W. Folkner using the standard IAU76 precession model. The IAU76 precession model, in the sense of physics/dynamics, was "baked into" DE431 via the lunar tidal model.
Giorgini worked with Jay Lieske, who developed the IAU76 precession model. They included a higher precision form of it than was published and adopted by the IAU (or used in DE431).
Then, outside the 1799-2200 interval, Owens long-term numerical integration model is used, as indicated in the documentation. Lieske was Owens supervisor for the development of the long-term model.
The Chebyshev polynomials published by Owens were converted and stored in a binary direct-access form Horizons accesses to construct a transformation matrix at any instant that transforms from the planetary ephemeris coordinates to a true-equator and equinox of date.
I'm not totally sure, but I'd like to try to duplicate the above figure 5-1 from Owen's dissertation. If I can duplicate this for the date range of -500,000 to 500,000 about J2000 as Owen did, by using data from HORIZONS or other, then I think we can conclude that we can use such a model.
py-NASA-horizons and HorizonJPL (docs) provide nice Python wrappers to the telnet interface. pip install horizonjpl
.
from horizon.interface import Interface
i = Interface()
i.get_version()
'3.95'
earth = i.get(399)
print(earth['Name'])
Earth
i.query(399)
{'Atm. pressure': '1.0 bar', 'Atmos': '5.1 x 10^18 kg', 'Escape velocity': '11.186 km s^-1', 'Fluid core rad': '3480 km', 'Geometric albedo': '0.367', 'ID': '399', 'Inner core rad': '1215 km', 'Love no., k2': '0.299', 'Magnetic moment': '0.61 gauss Rp^3', 'Mass, 10^24 kg': '5.97219+-0.0006', 'Name': 'Earth', 'Sidereal period': '365.25636 days', 'Vis. mag. V(1,0)': '-3.86'}
After Giorgini cleared up my misconceptions, I realize that I can not use HORIZONS in the original way I wanted to. However, I found a second long-term precession model that seemed easier for an astrophysics layman. This follow-up, vondrak.ipynb (nbviewer), implements a long-term precession matrix from beginning to end.
A helpful read is The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models (2005).