Precession, Proper Motion, and Epochs for Ancient or Future Observers

I have a question regarding epochs. I want to show the night sky as a super-ancient observer would see it or someone in the far flung future. Should I set the year and epoch to be the same?

I have an ipython notebook here (proper-motion.ipynb / nbviewer) with the equivalent code.

The following discussion on equatorial telescopes in the precession and epochs tutorial and the below example on Halley's Comet seem sure thing to me to use both year and epoch as the same. However, I found some conflicting results from a few other people.

When planning to observe with an equatorial telescope, you may want to use the current date as your epoch, because the rotation of the sky above your telescope is determined by where the pole points today, not where it pointed in 2000 or some other convenient epoch. Computing positions in the epoch of their date is accomplished by simply providing the same argument for both date and epoch:

>>> j = ephem.Jupiter()
>>> j.compute(epoch=ephem.now())   # so both date and epoch are now
>>> print("%s %s" % (j.a_ra, j.a_dec))  # doctest: +SKIP
8:44:29.49 19:00:10.23
>>> j.compute('2003/3/25', epoch='2003/3/25')
>>> print("%s %s" % (j.a_ra, j.a_dec))
8:43:32.82 19:03:32.5

uma_proper_motion

Ian Ridpath gives a talk, Pictures in the sky: the origin and history of the constellations, that displays an animation of the same transition. See the video at t=2563. Jones got his images from Rick Pogge's video and lecture which was output from his Fortran code.

I read the PyEphem tutorial section on fixed objects, precession, and epochs. This computes the position of Halley's Comet in 1066.

halley.compute(epoch='1066')
This is probably useless: it computes the current position of halley, but returns coordinates relative to the direction the earth’s axis was pointing in the year 1066. Unless you use a Conquest-era star atlas, this is not useful.

halley.compute('1066', epoch='1066')
This is slightly more promising: it computes the position of halley in 1066 and returns coordinates for the orientation of the earth in that year. This might help you visualize how the object was positioned above contemporary observers, who considered it an ill omen in the imminent conflict between King Harold of England and William the Bastard. But to plot this position against a background of stars, you would first have to recompute each star’s position in 1066 coordinates.

halley.compute('1066')
This is what you will probably use most often; you get the position of halley in the year 1066 but expressed in the 2000 coordinates that your star atlas probably uses.

I am a little confused by Ridpath's, Pogge's and Jones' (paper here) results. Although I am happy that we can get the stellar dustpan and butterfly net they talk about, why are they using coordinates relative to our current position? Wouldnt using the ancient epoch and far flung future epoch and year give results for how ancient or future man would see these stars instead? Based on the second Halley's Comet example Brandon Rhodes gives, I am going to go with the following as what I think the depictions of the precession of the Big Dipper are.

uma_transition

%pylab inline
import ephem

UMa = {
'dubhe': 'αUMa', #HIP 54061
'merak': 'βUMa', #HIP 53910
'phecda':'γUMa', #HIP 58001
'megrez':'δUMa', #HIP 59774
'alioth':'εUMa', #HIP 62956
'mizar': 'ζUMa', #HIP 65378
'alcor': '80UMa',#HIP 65477
'alcaid':'ηUMa', #HIP 67301
}

def const(const,year=None,epoch='2000',title=None):
    s = []

    for star in const.keys():
        s.append(ephem.star(star.capitalize()))

    for i in range(len(s)):
        if(year!=None):
            s[i].compute(year,epoch=epoch)
        else:
            s[i].compute(epoch=epoch)

    tau = 2.0 * pi
    degree = tau / 360.0
    hour = tau / 24.0

    ra_list = [star.a_ra / hour for star in s]
    dec_list = [star.a_dec / degree for star in s]

    mag_array = np.array([star.mag for star in s])
    size_array = (5 - mag_array) ** 1.5 * 4

    scatter(ra_list,dec_list,size_array)
    if(title!=None):
        pyplot.title(title)
    gca().xaxis.grid(True)
    gca().yaxis.grid(True)
    gca().invert_xaxis()
    return s

The function can be called with the following to test:

Setting Epoch and Year to Same Thing:

ancient_uma = const(UMa,year='-100000',epoch='-100000',title='Desk Lamp')
present_uma = const(UMa,year='2014',epoch='2014',title='Current Big Dipper (Epoch 2014)')
future_uma = const(UMa,year='100000',epoch='100000',title='Kite')

Only Changing Year and Leaving Epoch as 2000:

ancient_uma = const(UMa,year='-100000',title='Butterfly Net')
present_uma = const(UMa,year='2014',title='Big Dipper (Epoch 2000)')
future_uma = const(UMa,year='100000',title='Dustpan')

Not sure how helpful this is, but:

The epoch only affects the translation from right ascension and declination to altitude and azimuth. It doesn't affect how the stars or constellations actually look in the sky.

For example, if you know a star's RA and DEC in J2000 coordinates, and convert to alt/az using the standard formulas, your values will be a bit off for today, although not enough to make a difference for a visual observer. If you use these alt/az values to point your telescope, however, you may miss the object you seek or at least see if off-center (depending on the field width of your eyepiece).

The changing shape of the Big Dipper you're seeing is due to two effects:

Proper motion: the stars in the Big Dipper (and all other stars) are moving in physical space. Even if the Earth's pole didn't precess or nutate, the Big Dipper would change shape.
Projection: if you're projecting a sphere into 2 dimensions, there is always some distortion. Try projecting so that the Big Dipper is in the center of the projection to minimize distortion.

Let me know if that helped?

On Sat, 15 Nov 2014, Tommy wrote:

Date: Sat, 15 Nov 2014 11:56:51 -0800 From: Tommy notifications@github.com Reply-To: brandon-rhodes/pyephem <reply+0004c41dea837d5cdf3a4ee39ce7e5d4b81fbd6aeac26a7092cf00000001107f720 392a169ce02ea9a7a@reply.github.com> To: brandon-rhodes/pyephem pyephem@noreply.github.com Subject: [pyephem] Precession, Proper Motion, and Epochs for Ancient or Future Observers (#61)

I have a question regarding epochs. I want to show the night sky as a super-ancient observer would see it or someone in the far flung future. Should I set the year and epoch to be the same?

I have an ipython notebook here (proper-motion.ipynb) with the equivalent code.

The following discussion on equatorial telescopes in the precession and epochs tutorial and the below example on Halley's Comet seem sure thing to me to use both year and epoch as the same. However, I found some conflicting results from three other people.

When planning to observe with an equatorial telescope, you may want to use the current date as your epoch, because the rotation of the sky above your telescope is determined by where the pole points today, not where it pointed in 2000 or some other convenient epoch. Computing positions in the epoch of their date is accomplished by simply providing the same argument for both date and epoch:
>>> j = ephem.Jupiter()
>>> j.compute(epoch=ephem.now())   # so both date and epoch are now
>>> print("%s %s" % (j.a_ra, j.a_dec))  # doctest: +SKIP
8:44:29.49 19:00:10.23
>>> j.compute('2003/3/25', epoch='2003/3/25')
>>> print("%s %s" % (j.a_ra, j.a_dec))
8:43:32.82 19:03:32.5
Ian Ridpath gives a talk, Pictures in the sky: the origin and history of the constellations, that displays an animation of the same transition. See the video at t=2563. And I think Jones got his images from Rick Pogge's video which was output from his Fortran code.

I read the PyEphem tutorial section on fixed objects, precession, and epochs. This computes the position of Halley's Comet in 1066.

halley.compute(epoch='1066')
This is probably useless: it computes the current position of halley, but returns coordinates relative to the direction the earth?s axis was pointing in the year 1066. Unless you use a Conquest-era star atlas, this is not useful.

halley.compute('1066', epoch='1066')
This is slightly more promising: it computes the position of halley in 1066 and returns coordinates for the orientation of the earth in that year. This might help you visualize how the object was positioned above contemporary observers, who considered it an ill omen in the imminent conflict between King Harold of England and William the Bastard. But to plot this position against a background of stars, you would first have to recompute each star?s position in 1066 coordinates.

halley.compute('1066')
This is what you will probably use most often; you get the position of halley in the year 1066 but expressed in the 2000 coordinates that your star atlas probably uses.

I am a little confused by Ridpath's, Pogge's and Jones' (paper here) results. Although I am happy that we can get the stellar dustpan and butterfly net they talk about, why are they using coordinates relative to our current position? Wouldnt using the ancient epoch and far flung future epoch and year give results for how ancient or future man would see these stars instead? Based on the second Halley's Comet example Brandon Rhodes gives, I am going to go with the following as the more accurate depictions of the precession of the Big Dipper.
%pylab inline
import ephem

UMa = {
'dubhe': '?UMa', #HIP 54061
'merak': '?UMa', #HIP 53910
'phecda':'?UMa', #HIP 58001
'megrez':'?UMa', #HIP 59774
'alioth':'?UMa', #HIP 62956
'mizar': '?UMa', #HIP 65378
'alcor': '80UMa',#HIP 65477
'alcaid':'?UMa', #HIP 67301
}

def const(const,year=None,epoch='2000',title=None):
   s = []

   for star in const.keys():
       s.append(ephem.star(star.capitalize()))

   for i in range(len(s)):
       if(year!=None):
           s[i].compute(year,epoch=epoch)
       else:
           s[i].compute(epoch=epoch)

   tau = 2.0 * pi
   degree = tau / 360.0
   hour = tau / 24.0

   ra_list = [star.a_ra / hour for star in s]
   dec_list = [star.a_dec / degree for star in s]

   mag_array = np.array([star.mag for star in s])
   size_array = (5 - mag_array) ** 1.5 * 4

   scatter(ra_list,dec_list,size_array)
   if(title!=None):
       pyplot.title(title)
   gca().xaxis.grid(True)
   gca().yaxis.grid(True)
   gca().invert_xaxis()
   return s
The function can be called with the following to test:

Setting Epoch and Year to Same Thing:

ancient_uma = const(UMa,year='-100000',epoch='-100000',title='Desk Lamp')

present_uma = const(UMa,year='2014',epoch='2014',title='Current Big Dipper (Epoch 2014)')

future_uma = const(UMa,year='100000',epoch='100000',title='Kite')

Only Changing Year and Leaving Epoch as 2000:

ancient_uma = const(UMa,year='-100000',title='Butterfly Net')

present_uma = const(UMa,year='2014',title='Big Dipper (Epoch 2000)')

future_uma = const(UMa,year='100000',title='Dustpan')

Reply to this email directly or view it on GitHub: https://github.com/brandon-rhodes/pyephem/issues/61

@barrycarter hey Barry! I definitely like your idea of projection, I'll probably play around with creating a celestial sphere of differing radii or maybe even polygonal celestial spheres from the earth in the future. :) However, changing the epoch definitely changes how the stars look. I am depending on this proper motion / precession. It is what I desire for my depictions.

For example, let's just look at all of the combinations of -100000 years (BCE). I wrote a second IPython notebook here (nbviewer) that will help visualize the differences. I just don't know whether I should use the second, or the third choice.

my_star.compute(epoch='-100000')
my_star.compute(-100000',epoch='-100000')
my_star.compute('-100000')

There's a lot of code in the big-dipper.ipynb example, so I wont post it in the tracker, but I'll make a blog post on it soon. Here is the output of all of the scenarios though.

Present Day: Computes the current position of the Big Dipper stars and uses coordinates relative to the direction the earth's axis is pointing currently.

uma

Present Day, Old Coords: Computes the current position of the Big Dipper stars but uses coordinates relative to the direction the earth's axis was pointing at year -100000. This is rather useless unless you are wondering how an ancient alien on earth was pondering the future maybe?

uma_epoch

Ancient Day, Ancient Coords: Computes positon of the Big Dipper as it was in the year -100000. It returns coordinates for the orientation of the earth that year. At the time of writing, I believe this is the representation I want regarding an ancient ancient observer looking at the night sky.

uma_epoch_and_year

Ancient Day, Present Coordinates: This gets the positon of the Big Dipper in the year -100000 but expressed in 2014 coordinates. To me this is useless, but am I wrong?

uma_year

Keeping Epoch Same, Varying Year

uma_changing_year

Keeping Year Same, Varying Epoch

uma_changing_epoch

Varying Epoch and Year Together

uma_changing_epoch_and_year_together

I believe that the above depiction is the one I want regarding an ancient observer seeing the night sky in his perspective, and a futuristic observer seeing the night sky in his own perspective. This is rather hard to comprehend sometimes, so I would absolutely love some community input on this subject.

Do the coordinates on your graph represent right ascension and declination? Even accounting for precession, I don't think the Big Dipper ever crosses south of the celestial equator, since it's fairly close to the north ecliptic pole.

Yes they do! Personally I think it is crazy. Non-scientifically when it is crazy vs elegant, I'd just automatically say elegant is the right choice.. but I am basing my choice solely off of the Halley Comet descriptions. I made a GIF of the big dipper through time, re-computing for every 10000 years and holy moley it goes all over the place in the "Varying Epoch and Year Together" scenario. In the "Keeping Epoch Same, Varying Year" it stays largely in the same place.

Thanks for reporting this issue — and especially for providing a notebook of examples, letting me see exactly what you are seeing! The extra effort that you took to make your code available in notebook form made it much easier for me to dig in and figure out what is going on.

Here is what is going on:

When the date argument to compute() (the first argument) changes, the stars move, because they have proper motion. This is the only argument that “has the right” to move stars around relative to each other: only the real passage of time can give them the chance to change position.
The epoch argument causes zero relative motion of the stars relative to each other. Which is what we would expect: the shape of the Dipper on, say, January 1, 1970, is going to be the same whatever crazy system of coordinates you use to express it.
In the diagrams beneath your “Keeping Year Same, Varying Epoch” heading, it is very difficult to see that the stars are not moving, because the projection is moving the stars around so much on its own as the celestial pole moves closer to them then farther away.
Finally, the precession model does go crazy at about 20,000 years before or after the present — hence the crazy results like the Dipper in the southern sky!

I have put together a notebook of my own that steps through these results. You will be happy to see that the Basemap library for Matplotlib makes it very easy to switch to a projection that involves very little distortion around the Dipper if you move the Dipper to the diagram’s center:

http://nbviewer.ipython.org/github/brandon-rhodes/pyephem/blob/master/issues/github-issue-61.ipynb

The result of all this, at least for this issue, is that I probably need to put a guard around the precession() call or maybe compute() itself so that PyEphem gives a clear error if the precession calculation is pushed passed its boundary. I will look around for some documentation on what, officially, is considered its working range.

Thanks again!

I have found what looks to be the most recent big paper on the precession theory:

http://www.aanda.org/articles/aa/full/2005/10/aa1908/aa1908.html

The only statement about accuracy that I can find is:

This solution (the accuracy of which is estimated to be about 0.05 mas/cy over a two-millennium interval centered on J2000) is thus proposed as a replacement for the current IAU 1976 model (Lieske et al. 1977) for the precession of the ecliptic.

I see no statements about its accuracy outside of this rough range of AD 1000 through AD 3000, and so I am left without a clear idea of the range of dates to which PyEphem and Skyfield ought to restrict their precession calculations.

HORIZONS will give precession data from 10000 BC to 10000 AD, if that helps any.

On Tue, 18 Nov 2014, Brandon Rhodes wrote:

Date: Tue, 18 Nov 2014 14:05:48 -0800 From: Brandon Rhodes notifications@github.com Reply-To: brandon-rhodes/pyephem <reply+0004c41df0129fc9046b49a7f30d300ff71505dbadf6c46e92cf00000001108384b c92a169ce02ea9a7a@reply.github.com> To: brandon-rhodes/pyephem pyephem@noreply.github.com Cc: barrycarter github@barrycarter.info Subject: Re: [pyephem] Precession, Proper Motion, and Epochs for Ancient or Future Observers (#61)

I have found what looks to be the most recent big paper on the precession theory:

http://www.aanda.org/articles/aa/full/2005/10/aa1908/aa1908.html

The only statement about accuracy that I can find is:

This solution (the accuracy of which is estimated to be about 0.05 mas/cy over a two-millennium interval centered on J2000) is thus proposed as a replacement for the current IAU 1976 model (Lieske et al. 1977) for the precession of the ecliptic.

I see no statements about its accuracy outside of this rough range of AD 1000 through AD 3000, and so I am left without a clear idea of the range of dates to which PyEphem and Skyfield ought to restrict their precession calculations.

Reply to this email directly or view it on GitHub: https://github.com/brandon-rhodes/pyephem/issues/61#issuecomment-63554239

@barrycarter I should probably do some iterative tests to narrow down the range over which HORIZONS will make those predictions, and perhaps adopt the same limits in my own libraries.

Wow! That is really fascinating. Thank you for digging into this.

Summary of conclusions:

The answer to my original question should I vary epoch and time together for an ancient observer on earth or an observer in the far-flung future? is yes.
The IAU precession formulas the code is currently based on gives an accuracy of 0.05 milliarcseconds per cycle from AD 1000 to AD 3000. Past that point we aren't sure (yet). See Improvement of the IAU 2000 precession model (2004).
A quick eyeball shows that anything past J2000 ± 20000 years is completely erratic.
Take stellar depictions involving precession with a grain of salt. Apparently a lot of popular programs don't account for our earth's changing orientation.
Basemap is much easier to use than plotting using an equirectangular projection and provides less distortion.

I looked into Barry's tip. Thanks! We're getting into stuff that is way beyond me calculation-wise, but I did find that you can select time spans from 'BC 9998-03-20 to AD 9999-12-31' for the ephemeris.

HORIZONS uses two different models. One between 1799-Jan-1 to 2202-Jan-1, and the other for outside those bounds.

PRECESSION MODEL:

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)

Thank you! I will take a look at these references and see whether Skyfield could incorporate a similar approach to that of Horizons.

I did a wee bit of research on Owen's long-term precession model and the bounds of BC 9998-03-20 to AD 9999-12-31 for the earth ephemeris using the HORIZONS interface. My background is electrical engineering; astronomy has been a hobby of late, so please feel free to comment on any inaccuracies or misconceptions I have. I didn't know where to place this, but _The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models_ (2005) seems like a good read on these subjects.

On HORIZONS Earth Ephemeris DE431MX

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. I need someone more qualified to weigh in, but I think we can ignore the upper and lower limits of B.C. 9998-03-20 to A.D. 9999-12-31. These are just for the ephemeris and not limits tied to 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.

Thank you Mr. Giorgini for the quick and thorough reply!

On A Theory of the Earth's Precession Relative to the Invariable Plane of the Solar System

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.

Pg 253:

precession

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 long-term precession model by Vondrák was a bit easier for me to wrap my brain around (and less coefficient typing!). I translated his Fortran code into Python in vondrak.ipynb (nbviewer).

If I wanted to supplement the libastro precession function, can I? and still use ephem.compute?

Thanks for doing all of this research! It will take me another week or two to absorb, but I do appreciate it and it is a big help in expanding the range of dates over which these astronomy libraries can operate.

I will have to think about whether I want to be supplementing libastro with any further custom code — it is, so far as possible, supposed to be the libastro from XEphem with only those modifications made that are needed to make it work with Python. If we start changing its behavior in substantive ways, users will not know why PyEphem, given a certain date, starts returning different results than XEphem when they ask it about the same date.

I think it would make more sense to make these more modern precession algorithms available in Skyfield, which has no historical tie to XEphem and can be taken in any direction we want! (And is also written in pure Python, making development easier.)

Okay, that is a really good point about keeping PyEphem close to its libastro roots. I totally agree with that logic.

Time to play with Skyfield today :)

I am going to close this old PyEphem issue as I am convinced that precession solutions should be based on a more modern toolset than PyEphem.

brandon-rhodes / pyephem

Precession, Proper Motion, and Epochs for Ancient or Future Observers #61