Throughout test of ephemerides: JPARSEC vs Horizons, IMCCE, and Astronomical Almanac

# Throughout test of ephemerides: JPARSEC vs Horizons, IMCCE, and Astronomical Almanac

This month I have implemented more than 250 test cases of ephemerides in JPARSEC, using data from Astronomical Almanac (AA) and the ephemerides servers Horizons and IMCCE. I tested ephemerides for the last time a few years ago, and never in depth, so I expected to find some bugs in the results and it has been the case, but the errors are not very important. These are the main improvements in the ephemerides of the JPARSEC version uploaded today.

• The calculation of the heliocentric distance of planets were not completely correct (a secondary light-time correction was not implemented), and this had some little effects in the values of the elongations or plase angles. Both values were dependent of the ephemeris type (apparent, astrometric, and geometric).
• The astrometric position of the Moon was always wrong by a few arcseconds due to a bug in the correction for aberration.
• The sign of the subsolar longitude of the Moon was wrong.
• Ephemerides for asteroids/comets never show orientation of axis/sizes. This occurs for those comets/asteroids with known parameters of the orientation of their rotation axis.
• The calculation of LSR velocities in StarEphem were wrong, and B1950 to J2000 transformation was not completely accurate/well documented.

In addition to these bug fixes, new features have been implemented in JPARSEC. The main ones are listed next.

• Ephemerides calculations now supports two new reduction algorithms: the IAU 1976 method (Lieske precession and obliquity, among other algorithms), required to match results from IMCCE server, and IAU 2006 resolutions to fit AA values. Finally, IAU 2006 resolutions are supported, and some methods to work with matrices (like the NPB matrix) are available, although the implementation in the ephemerides is not based in matrices.
• The orientation of the rotation axis of planets and satellites are supported for IAU 2000, 2006, and 2009 resolutions, independently. The selected orientation parameters depends on the selected reduction algorithm. IAU 2006 improves the accuracy of the rotation of some planets, but IAU 2009 resolutions changes too much the longitude of the central meridian of Jupiter and Pluto. Since the reduction itself is the same (RA and DEC), I think the best current method is IAU 2006.
• JPLEphemeris class has been redesigned, adding support for DE422, improving the speed of calculations (Chebyshev polynomials are now stored in memory), and allowing the use of external files. DE102 no longer available.
• The rise, set, and transit times in the EphemElement object are now arrays that can hold different values in the case that a given object has more than 1 rise/set/transit in a given day, or none at all. Calculations in the RiseSetTransit class will automatically add a new rise/set/transit time after the last one calculated.
• JPL ephemerides were not used as reference to calculate ephemerides for natural satellite in the accurate mode in Ephem class. Now this is done automatically if DE4xx files are detected to be available.

A little set of limitations has been identified and could require some kind of fix in the future. These are:

• Moon orientation seems more accurate with Eckardt theory than with IAU resolutions, since Eckardt's theory is much more compatible with the results of JPL ephemerides. I've kept IAU resolutions by default, but the IAU model for the rotation of the Moon seems to be outdated. I expected a difference of 0.03 deg according to documentation, but found 0.7. An accurate reduction of nutation/librations should be implemented in JPLEphemeris and PlanetEphem classes.
• Planetary positions are respect to barycenter (center of mass). For instance, there's an offset of 2000 km in Pluto, which is noticeable in the RA/DEC values (0.07”) despite the distance of this object.
• The getConstellation method in Constellation class should be used with astrometric positions, instead of apparent ones. This means that the constellation where an object lies could be wrong in very extreme cases when the limit to the next constellation is lower than 30” (aberration + nutation).

But the main effort has been to test ephemerides with different sources. The official ephemerides theory used in the astronomical community is the DE405 integration, but it is noticeable the difference between modern JPL integrations (DE422) and DE405, which is clearly becoming outdated. For instance, the difference in the position of Jupiter (which is not the object where this difference is the highest) varies from 0.05” at J2000 to about 0.15” at J1900. It is also noticeable the extremely old implementation of ephemerides reduction algorithms in astronomical ephemerides servers (Horizons and IMCCE), that produces differences in the RA/DEC from 0.05” at J2000 to 0.3-0.4” at J1900 (mainly due to Lieske's precession), compared to IAU 2006 resolutions. These differences explain the dificulties when comparing the results from different sources and when trying to match them with JPARSEC or test them accurately.

## The Horizons ephemerides server

It is available at http://ssd.jpl.nasa.gov/horizons.cgi. After testing it with JPARSEC it seems to have a very good implementation of the geometry of the Solar System, which means that the values of parameters like planetary elongations and phase angles, or longitudes of the central meridians are fine. I've been unable to fit the RA/DEC values to the milliarcsecond, probably because I don't know the reduction algorithms used internally. It seems to use some IAU 1976 algorithms (precession), but not all of them or at least in a non-standard way. Maybe an accurate apparent RA/DEC output is not the main object. I've found also unexpected results for the longitude of the central meridian in Saturn and Pluto (should follow IAU 2006 resolutions as in other bodies, but it doesn't), and also a possible bug: longitudes of central meridian seems to be wrong far from J2000 (even after or before a few years the difference is noticeable), compared to IMCCE and JPARSEC. There's also a strange value for the distance to the Moon. In fact, Horizons gives two different values for this distance depending on the ephemeris type (VECTOR and OBSERVER), being the one for VECTOR type what I can obtain. Since the difference can reach 100 km (much greater than the difference of 2 km between the geometric and barycentric centers of the Moon), maybe the other is the distance that light cross to rearch the observer (corrected for relativistic effects), but it is confusing. It is also noticeable the difference in the distance to the planets when selecting or not the (quite hidden) extra precision flag, that should be checked to get correct results for all those decimal places given.

## The IMCCE ephemerides server

It is available at http://www.imcce.fr/en/ephemerides/generateur_ephemerides.php. In this case the object is clearly to produce accurate RA/DEC output, but the point that output positions are respect to the ICRF takes no sense considering the effects of using these old algorithms instead of those for IAU 2000/2006 resolutions. In addition, according to IAU resolutions ICRF-based positions should be used with astrometric ephemerides, and apparent places of objects should be based on the mean dynamical equinox of J2000. IAU 1976 standard reduction algorithms are used, but close to J2000 the reduction is different and seems to follow more recent techniques (either the APPLY_WILLIAMS or the APPLY_JPLDE40x set of algorithms implemented in JPARSEC). IAU 2000 resolutions are used to calculate the orientation of the planets, but these values are often wrong, sometimes in the sign of the longitude of the central meridian and in other cases in the value itself. It is noticeable the option to generate a chart with the orientation of the planet, that seems to show always the correct value for the longitude of the central meridian, that as I say is sometimes different from that of the table. I've also seen wrong values for the elongations sometimes. The IMCCE server has also an option to access 'General ephemeris of natural satellites position', that produces extended calculations for natural satellites compared to the main option of Solar System bodies. In this case the output far from J2000 for Martian satellites (only test I have done here) is wrong far from J2000, compared to the output positions of the main page. IMCCE seems the right server to test RA/DEC output to the milliarcsecond level, which clearly has been well tested, but the rest is a little buggy.

## Results of the tests

I have implemented different tests to ensure the quality of JPARSEC. A first test of the rectangular coordinates for different planetary and satellite theories to match the results from Horizons and from the documentation of the IMCCE theories (see ftp://ftp.imcce.fr/pub/ephem) gives as a result a very little difference in the positions of the Uranian satellites using GUST86 theory (< 1 km), that I will see in detail to fix it if necessary. Another test using Horizons is implemented to test apparent/astrometric/geometric RA/DEC to 0.1” of accuracy, and other values (elongation, phase angle, subsolar and subearth positions). This test is repeated for all theories in JPARSEC (VSOP87, ELP2000, every JPL integration, Moshier method) to ensure basically that correct results are always obtained for every field of an EphemElement object when the test is not done to the milliarcsecond (mas) level. Next I extended the tests to the milliarcsecond in planets and satellites, trying to 'fit' the results from Horizons, IMCCE, and AA. IAU 1976 and 2006 resolutions were implemented for this task, and I found it easy for AA and IMCCE and planets, and a little tricky but possible with natural satellites. Uranian satellites with GUST86 were already perfect to the mas level. In Saturn and using TASS 1.7 I had to rotate ecliptic positions from TASS to equatorial using the ecliptic of date (even though TASS positions are J2000). With L1 theory in Jupiter the discrepancy increases also from J2000, so some kind of correction seems to be required, although the difference is of a few tenths of mas at most. Martian satellites results departs from JPARSEC results by 2” in 4 centuries, and JPARSEC results seems wrong here since Horizons gives similar output compared to IMCCE. Maybe the public version of the theory for the Martian satellites cannot be used far from J2000.

Below there's a set of comparative results obtained with JPARSEC, IMCCE, and HORIZONS.

 Date and time (TT) Method Body RA, DEC, DISTANCE COMMENT 1600, Jan 1, 0h IMCCE Mars 09h 24m 41.58178s, 19º 15' 18.9756”, 0.746598316 JPARSEC 09h 24m 41.58188s, 19º 15' 18.9746”, 0.746598317 Using IAU 1976 reduction algorithms and ICRF frame 1900, Jan 1, 0h IMCCE Mars 19h 00m 39.90359s, -23º 38' 58.3076”, 2.400963430 JPARSEC 19h 00m 39.90409s, -23º 38' 58.3100”, 2.4009634300 2200, Jan 1, 0h IMCCE Mars 10h 25m 57.25433s, 13º 29' 37.0809”, 0.842946359 JPARSEC 10h 25m 57.25447s, 13º 29' 37.0803”, 0.84294636204 1600, Jan 1, 0h IMCCE Neptune 10h 00m 31.59378s, 12º 51' 13.3640”, 29.466102577 JPARSEC 10h 00m 31.59382s, 12º 51' 13.3634”, 29.4661025774 1900, Jan 1, 0h IMCCE Neptune 05h 39m 21.87494s, 22º 03' 58.9701”, 28.920271361 JPARSEC 05h 39m 21.87494s, 22º 03' 58.9700”, 28.9202713608 2200, Jan 1, 0h IMCCE Neptune 01h 24m 58.76536s, 07º 07' 44.9398”, 29.615237432 JPARSEC 01h 24m 58.76543s, 07º 07' 44.9404”, 29.6152374318 2000, Jan 1, 0h IMCCE Neptune 20h 21m 39.47908s, -19º 13' 01.9774”, 31.021117421 2000, Jan 1, 0h HORIZONS Neptune 20h 21m 39.4756s, -19º 13' 01.987”, 31.0211174207 JPARSEC 20h 21m 39.47555s, -19º 13' 01.9870”, 31.0211174207 (IAU 1976 reduction algorithms) JPARSEC 20h 21m 39.47906s, -19º 13' 01.9774”, 31.0211174207 Using JPLDE40x reduction algorithms 2011, July 18, 0h AA Mars 05h 11m 32.063s, 23º 05' 16.56”, - JPARSEC 05h 11m 32.06253s, 23º 05' 16.5596”, 2.17658764475 Using IAU 2006 reduction algorithms and J2000 frame 2011, July 18, 0h AA Jupiter 02h 21m 54.955s, 12º 49' 28.91”, - JPARSEC 02h 21m 54.95512s, 12º 49' 28.9066”, 5.06967777879 2011, July 18, 0h AA Neptune 22h 10m 46.150s, -11º 49' 46.06”, - JPARSEC 22h 10m 46.15016s, -11º 49' 46.0553”, 29.1731323295 1600, Jan 1, 0h IMCCE Deimos 09h 24m 44.52675s, 19º 15' 19.6628”, 0.746640479 JPARSEC 09h 24m 44.58900s, 19º 15' 18.2332”, 0.74662669602 Using IAU 1976 reduction algorithms and ICRF 1900, Jan 1, 0h IMCCE Deimos 19h 00m 40.56578s, -23º 38' 59.5043”, 2.400848638 JPARSEC 19h 00m 40.56644s, -23º 38' 59.5203”, 2.4008486685 2200, Jan 1, 0h IMCCE Deimos 10h 25m 59.52983s, 13º 29' 25.4275”, 0.842883773 JPARSEC 10h 25m 59.45245s, 13º 29' 25.0445”, 0.84287562191 1600, Jan 1, 0h IMCCE Titan 13h 45m 15.35365s, -08º 17' 57.0768”, 9.99873438 JPARSEC 13h 45m 15.35371s, -08º 17' 57.0768”, 9.99873438030 Using ecliptic of date, not J2000 ecliptic 1900, Jan 1, 0h IMCCE Titan 17h 50m 10.06797s, -22º 24' 27.4608”, 11.031604402 JPARSEC 17h 50m 10.06800s, -22º 24' 27.4608”, 11.0316044020 2200, Jan 1, 0h IMCCE Titan 22h 02m 43.88507s, -13º 27' 15.7030”, 10.450165195 JPARSEC 22h 02m 43.88515s, -13º 27' 15.7028”, 10.4501651962 1600, Jan 1, 0h IMCCE Callisto 09h 33m 53.33724s, 15º 30' 35.5888”, 4.558908178 JPARSEC 09h 33m 53.33510s, 15º 30' 35.6000”, 4.55890833039 1900, Jan 1, 0h IMCCE Callisto 15h 56m 43.78523s, -19º 36' 26.0122”, 6.125545509 JPARSEC 15h 56m 43.78401s, -19º 36' 26.0106”, 6.12554669865 2200, Jan 1, 0h IMCCE Callisto 22h 26m 48.97572s, -10º 48' 53.3763”, 5.494266711 JPARSEC 22h 26m 48.97557s, -10º 48' 53.3779”, 5.49426633534 1600, Jan 1, 0h IMCCE Oberon 01h 41m 49.40287s, 10º 01' 47.6674”, 19.470081383 JPARSEC 01h 41m 49.40292s, 10º 01' 47.6673”, 19.47008138146 1900, Jan 1, 0h IMCCE Oberon 16h 34m 01.71001s, -21º 54' 48.8151”, 19.838041698 JPARSEC 16h 34m 01.71002s, -21º 54' 48.8149”, 19.83804168796 2200, Jan 1, 0h IMCCE Oberon 05h 46m 52.47223s, 23º 32' 39.3707”, 18.121971509 JPARSEC 05h 46m 52.47231s, 23º 32' 39.3711”, 18.12197150938

Tests were implemented for other parts or the library, which resulted in different minor corrections. Everything seems to be fine now, but there's still margin to improve. B1950 to J2000 transformation can maybe be improved according to J. Bennett (see http://ned.ipac.caltech.edu/help/calc_doc.txt). I would like to test eclipses in more detail, and also the MainEvents class. There are some fails in the tests for occultation of stars by planets, mutual events of natural satellites, and the ephemerides for Triton. Other fails mentioned above include the orientation of Saturn and Pluto axes, GUST86 positions, and ephemerides of Martian, Jovian, and maybe Saturnian satellites to the mas level. Perhaps I will have to contact people at IMCCE or Horizons to solve these fails. Other things remain like the ephemerides of dwarf satellites, maps of solar and lunar eclipses, or to test the polar movement correction, among others, and I would like to improve code readability and documentation even more.

## Discussion

 ______   ___    ____   ___    _  __
/_/    /_/ |_|/_/    /_/    /_/|_|