Add support for ecliptic_latlon() at additional epochs beyond J2000

[hjy1210]

As far as I know, spring equinox of year 2016 is '2016-3-20 4:30'.

I check it in HORIZONS Web, get result as follows

Date__(UT)__HR:MN        ObsEcLon    ObsEcLat
**********************************************
$$SOE
 2016-Mar-20 04:29     359.9991687  -0.0001637
 2016-Mar-20 04:30     359.9998585  -0.0001637
 2016-Mar-20 04:31       0.0005482  -0.0001637
$$EOE
*******************************************************************************

But when I execute following scripts

# Skyfield issue
from skyfield.api import load
ts = load.timescale()
t = ts.utc(2016,3,20,4,30,0)
planets = load('de421.bsp')
sun=planets['sun']
earth.at(t).observe(sun).apparent().ecliptic_latlon()[1].degrees

result is:

359.77397918133533

In summary, using SKyfield, I got longitude of sun at '2016-3-20 4:30' is 359.7740, not close enough to 360.

[brandon-rhodes]

The method ecliptic_latlon() produces coordinates relative to the J2000 ecliptic — which is why you are not finding that they line up with the equinox 16 years later. In 257d5ed I have modified the documentation to make this clear — thanks for pointing out this potential for confusion that is created by the ecliptic methods!

[ghost]

Thought: de421.bsp is actually fairly old and may generate some inaccuracy as well.

Also, IS there a method to get ecliptic_latlon of an object based on Jdate, where date is arbitrary?

[brandon-rhodes]

(a) Happily, all of the DE ephemerides are of extremely high precision for bodies like the Earth and Moon, for which we have had laser-ranging data since the Apollo landings. The difference between a 2009 ephemeris like DE421 and the more recent DE430, for example, is less than a milliarcsecond:

from skyfield.api import load
ts = load.timescale()
t = ts.utc(2016,3,20,4,30,0)

planets = load('de421.bsp')
sun=planets['sun']
earth=planets['earth']
lat, lon, distance = earth.at(t).observe(sun).apparent().ecliptic_latlon()
print(lon.dstr(4))

planets = load('de430t.bsp')
sun=planets['sun']
earth=planets['earth']
lat, lon, distance = earth.at(t).observe(sun).apparent().ecliptic_latlon()
print(lon.dstr(4))

359deg 46' 26.3246"
359deg 46' 26.3249"

(b) The USNO Circular 179 from 2005, on which much of Skyfield is conceptually based —

http://aa.usno.navy.mil/publications/docs/Circular_179.pdf

— indicates in its discussion of the new high-precision reference frames that, “However, there is much unfinished business. The apparently familiar concept of the ecliptic plane has not yet been defined in the context of relativity resolutions.” And since their own NOVAS library correspondingly lacks any code to convert to the ecliptic plane, I have not been sure about how to proceed. I will review my copy of the Explanatory Supplement to the Astronomical Almanac to see if I can find at least a low-precision algorithm for converting, and see whether it might be suitable for Skyfield.

And I will re-open this issue to track my progress!

[ghost]

The following may or may not be helpful: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/aareadme.txt (and parent directory for actual files).

If that's not helpful, here are some other things are also probably not helpful:

http://naif.jpl.nasa.gov/pub/naif/LADEE/kernels/fk/ladee_frames_2014159_v01.tf http://naif.jpl.nasa.gov/pub/naif/VEX/kernels/fk/RSSD0002.TF

The discussion regarding the circular you're referring to is (I think) related to:

https://naif.jpl.nasa.gov/pipermail/spice_discussion/2005-February/000157.html