Add fast algorithms to compute the position of the planets

[ronisbr]

We need a fast algorithm to compute the position of the planets so that we can use in a numerical orbit propagator.

[helgee]

I am working on a pure Julia version of VSOP87: https://github.com/JuliaAstro/AstroBase.jl/blob/master/src/ephemerides/Ephemerides.jl

Would that fit the bill?

It has the advantage of not needing any external data files (e.g. SPICE kernels) but it is actually slower than JPLEphemeris.jl right now.

[ronisbr]

Hi @helgee ,

Indeed, that will be perfect! Thanks!

Maybe, I can add the very, very simple approach on Vallado's book and let the user choose. However, it will be very good to use the VSOP87 due to the precision.

[helgee]

So, my VSOP87 implementation is complete and agrees with the original Fortran version in it's native frame (dynamical ecliptic frame J2000). They give a rotation matrix to FK5 but with that rotation applied there is still a rather large discrepancy to DE200 (to which VSOP87 was originally fitted).

Any ideas what this could be?

EDIT: This graph looks virtually identical for all bodies so this must be a frame alignment issue.

[ronisbr]

Hi @helgee ,

Can you please provide me more information about the inputs and the expected results? I might be able to help.

[helgee]

Thanks! I have been banging my head against this for a few days now...

Here is the information I have collected so far:

• There are several different versions of VSOP87 (A to E). I am using version E since it returns the state of the planets in barycentric rectangular coordinates whereas the others use heliocentric and/or spherical coordinates.

• The reference frame used by VSOP87E is described as the inertial frame defined by the dynamical equinox and ecliptic J2000 (JD2451545.0). I am not sure what dynamical means in this context. The README supplied with the data advises to use the following rotation matrix to transform from the abovementioned frame to equatorial FK5:

  X        +1.000000000000  +0.000000440360  -0.000000190919   X
  Y     =  -0.000000479966  +0.917482137087  -0.397776982902   Y
  Z FK5     0.000000000000  +0.397776982902  +0.917482137087   Z VSOP87A

• This leads to the results in my post above. It shows the difference between JPL DE200 (evaluated via JPLEphemeris.jl and my VSOP87 implementation with the rotation applied.

• Most other open-source implementations of VSOP87 follow the description in "Astronomical Algorithms" by Meeus which uses version C and provides the following correction for spherical coordinates:

  

• I also found this quote in some online forum about a closed-source astrodynamics software (https://www.orbiter-forum.com/showthread.php?p=499906&postcount=23):

Reference Frame: The VSOP87 ephemeris uses a coordinate system that isn't quite inertial. If I remember correctly, it rotates slowly with respect to FK5 and ICRF inertial reference frames. One or other of the latter inertial reference frames is probably used by GMAT - I don't know which - so you also need to make sure that you rotate coordinates from FK5/IRCF to the dynamical reference frame used by VSOP87. Another little nuisance with which to complicate life.

• Finally, there is an online version of VSOP87: MULTI-SAT. It has the following information in its user manual:

Here is the paper they referenced but I have not understood it yet 😜 : Standish 1982.pdf

[ronisbr]

Good @helgee ! I will process all this information to see if I can help you! Probably I will have time this week to do this.

[ronisbr]

Hi @helgee ,

That matrix that rotates the VSOP frame to FK5 seems correct. It is possible to obtain a very, very close approximation using the information available in [1].

Moreover, there is a good test to see if the problem is just a frame rotation. If both vector representations have the same origin, then those norms must be very close. Is it possible to plot the norm difference? Like, norm(v_de200) - norm(v_vsop87). If this is significantly lower than the first result, then it maybe an error in the reference frame.

However, I was thinking about this error and it does not seem to be that big at all. For example, according to [1], VSOP87 has an accuracy of 1 arcseg (I think so, I did not read the paper carefully). If you account the distance between the Sun and Mars, that angular error will be translated into a linear error of about 1,100 km. Hence, if my interpretation is correct, then it seems fine :)

[1] http://adsabs.harvard.edu/full/1988A%26A...202..309B