# Implementing the VSOP-2000 Ephemeris

For a project I'm doing, I need to implement the VSOP-2000 ephemeris for the Moon. However, unlike VSOP-87, there is little to no documentation online on how to use the data. For example, here is the first line of the x-term for the Moon in VSOP-2000:

  1   0  0  1  0   0  0  0  0   0  1  0  0   0  0  0  0   0  0  0  0     -.1053243070325313D-09     -.2560241101399878D-02


This is very different from VSOP-87, which I know how to use. Are there any resources online that could inform me on how to use this?

• Are you sure you need to use the VSOP-2000 ephemeris and you can not use another like one of the JPL development ephemerides, which are widely used and well supported and well documented? You really must use an ephemeris that has "little to no documentation online"? – uhoh Apr 6 '20 at 15:51
• @uhoh as long as I can implement it in Mathematica, and it gets the Moon's position to >10'' for a couple millennia forward and back, I don't really care what it is. I just heard VSOP 2000 could do this. – John Dumancic Apr 6 '20 at 17:55
• Okay I see. So SPICE (mentioned in this comment) uses the JPL development ephemerides (DE's) as does the Python package Skyfield. The DEs are tables of coefficients that are used with a special form of Chebyshev polynomial interpolation. See this answer I am sure Mathematica can do this but it may be quite a lot of work to do what is designed to do in code. – uhoh Apr 6 '20 at 20:45
• If you know any Python at all then Skyfield is by far the easiest way to do implement an eclipse finding algorithm. I don't know if it will be for sure accurate to 10 arcsec over millennia but the documentation for each of the various DEs will go into detail on the expected accuracy. They use all known data including the centimeter-level laser ranging of the Moon that's been going on regularly since the Apollo missions put corner cube reflectors on the Moon in the 1970's. – uhoh Apr 6 '20 at 20:49
• more about Chebyshev interpolation of DE's also in this answer, possibly also helpful are answers to Is there a way to extract the Chebyshev coefficients for a body from a SPICE kernel? – uhoh Apr 6 '20 at 20:51

Probably the most reliable guide to vsop2000-p11.dat is the reference implementation vsop2000.for. The input line in question is the first of 5399 k=0 records contributing to $$x$$. Those are followed by 5842 k=1 records, 610 k=2 records, and so on up to k=8. Then $$y$$ (k=0..8) and $$z$$ (k=0..7) follow the same pattern.

If I read the Fortran code correctly, the fixed-length (not blank-delimited) record format is:

• Columns 1-7 are an index $$i$$
• Columns 9-63 are 17 per-body coefficients $$b_{kij}$$
• Columns 64-72 are three zeroes, ignored
• Columns 73-126 are coefficients $$s_{ki}, c_{ki}$$ in astronomical units

To compute a position for time $$t$$ in Julian years since J2000:

1. For each body $$j$$, compute $$\theta_j(t) = \dot{\theta}_j t + \theta_j(0)$$, where $$\dot{\theta}_j$$ is the mean motion and $$\theta_j(0)$$ is the J2000 mean longitude from vsop2000-exp.txt
2. For each record $$x_{ki}$$, compute $$\phi_{ki} = \sum\limits_{j=1}^{17} b_{kij} \theta_j(t)$$
3. $$x = \sum\limits_{k=0}^8 t^k \sum\limits_{i=1}^{n_k} (s_{ki} \sin \phi_{ki} + c_{ki} \cos \phi_{ki})$$
4. For $$y$$ and $$z$$, repeat steps 2 and 3 using the respective sets of records.

This should give you a geocentric Moon position in astronomical units.

• A note to any hypothetical future readers: this works like a charm. – John Dumancic Apr 7 '20 at 19:22
• @Tesseract I modified vsop2000.for to give Moon positions for 2020-2023 in spherical coordinates accounting for precession. lat/lon differences from JPL HORIZONS (DE430) were on the order of 15". Moisson and Bretagnon 2001 compare their results to DE403. – Mike G Apr 7 '20 at 19:31
• In the sample line given by John Dumancic in his original question, we find e.g. “-.1053243070325313D-09.” What does the D stand for? I’d be expecting an E here, for E-09 as in one billionth, but a D…? – Pierre Paquette Nov 21 '20 at 20:35
• @PierrePaquette I think that's just double precision floating point, no other difference from E. – Mike G Nov 22 '20 at 1:55
• @PierrePaquette $n_k$ is the number of records in the .dat file for a given value of $k$. – Mike G Nov 27 '20 at 5:00