Ephemeris trail
Introduction
Generalities
VSOP87
Build VSOP87
Truncating VSOP87
Testing VSOP87
Pluto99
ELP82
Class organization
Reference frames
CoordTransfo
Handling frames
Handling precision
Units
Handling time
Putting all together
Swiss Ephemeris
Building VSOP87
The purpose of this page is to present how VSOP87 were incorporated in JEphem.

Choice of the version

There are six different versions of VSOP87 :
  • The main version (VSOP87) permits to compute the elliptical elements.
  • Versions A and C give heliocentric cartesian (rectangular) coordinates ;
  • Versions B and D give heliocentric spherical coordinates ;
  • Version E gives barycentric catresian coordinates ;

  • For all the versions, time argument is TDB, and reference frame is BRS.

    Our purpose is to compute :
    - heliocentric geometrical and apparent coordinates, and
    - geocentric apparent equatorial and ecliptic coordinates.
    So we'll have to perform coordinate transformations ; it's then simpler to work with cartesian coordinates and transform to spherical just before displaying the results. Versions A and C are then appropriate.
    For version A, the reference frame is defined by the dynamical equinox and ecliptic J2000.
    For version C, the reference frame is defined by the dynamical equinox and ecliptic of the date.
    (see farther for details about frames)
    Precession has been applied to go from the frame of version A to the frame of version C. I first thought that using version C would permit to avoid precession computation. We'll see farther that it's not possible (because matrix multiplication is not commutative).

    So we'll use version A (version D is used in XEphem, and version A in WinAstro).

    VSOP87 files

    In BDL server, we find one file per planet. Here is an extract of VSOP87 file (for Neptune, version A) :
     VSOP87 VERSION A1    NEPTUNE   VARIABLE 1 (XYZ)       *T**0    772 TERMS    HELIOCENTRIC DYNAMICAL ECLIPTIC AND EQUINOX J2000      
     1810    1  0  0  0  0  0  0  0  1  0  0  0  0 -0.00682678287    30.05889926953    30.05890004476 5.31211340029      38.13303563780 
     1810    2  0  0  0  0  0  0  0  0  0  0  0  0  0.00000000000    -0.27080164222     0.27080164222 3.14159265359       0.00000000000 
    etc...
     1810  771  0  0  0  0  6-16  0  2  0  0  0  0  0.00000001037    -0.00000001976     0.00000002232 2.42169508541     158.37366516480 
     1810  772  0  0  0  0  0  0 19-20  0  0  0  0  0.00000000619     0.00000002392     0.00000002471 3.93681150847     658.18966002270 
     VSOP87 VERSION A1    NEPTUNE   VARIABLE 1 (XYZ)       *T**2    102 TERMS    HELIOCENTRIC DYNAMICAL ECLIPTIC AND EQUINOX J2000      
     1812    1  0  0  0  0  0  0  0  0  0  0  0  0  0.00000000000     0.00005371138     0.00005371138 0.00000000000       0.00000000000 
     1812    2  0  0  0  0  0  0  1 -1  0  0  0  0  0.00004488541     0.00000656405     0.00004536283 5.02700751836      36.64856292950 
    etc...
     1812  101  0  0  0  0  1  0  0 -1  0  0  0  0  0.00000001759    -0.00000001439     0.00000002273 5.59760233612     491.55792945680 
     1812  102  0  0  0  0  0  3  0 -4  0  0  0  0 -0.00000002773    -0.00000000305     0.00000002790 1.90430778774     487.36514376280 
    etc..
    
    For each coordinate (here X, Y and Z), terms are given for each power of time (from 0 to 5).
    From VSOP87.doc, we can find out that coordinates can be computed in two different way ; the simplest way uses only the three last numbers of each row (called A, B, C).
    The formula to use is :
    Derivating with respect to time, we get the velocities :
    Where :
    - Ai, Bi, Ci are the terms in the BDL file.
    - n is the number of terms given for a planet and a power of time (for example, n = 772 for Neptune, for variable X, T 0).
    - T is in thousands of julian day from JD2000, T = (jd - JD2000) / 365250 (see the page dealing about time for more details).

    Building the java classes

    To compute the coordinates, we need :
  • a class which computes the summation : done by jephem.astro.solarsystem.vsop87.VSOP87 ;
  • files to store the data.


  • As we only need the terms Ai, Bi, Ci, the first thing to do is to remove the other terms from the original files. This has been done by a class which is not part of JEphem package : BuildVSOP87.java (you can get it from download area - use of BuildVSOP87.java is described in its javadoc page). BuildVSOP87.java permits to generate files containing terms Ai, Bi, Ci in 3 different formats : raw text; formatted java class, and binary format.

    Full precision version

    The ideal would be to format the data into java classes (a double[][][] array), and compile the classes. This would permit to keep the data in binary format, and let the data usable by applets (which can't read files). But a specification of JVMs (Java Virtual Machines) imposes a limit : the size of a compiled class must not exceed 64 Ko.

    So I used binary format to store the values, generating one file per planet. At execution time, the data are retrieved from the files. The class which performs the computation (jephem.astro.solarsystem.vsop87.VSOP87) stores the data in static variables, so file reading is done only at first computation.

    In JEphem hierarchy, the data files are stored in directory Jephem/data/astro/planets/vsop87/VSOP87A (see page 'JEphem directories' of informatic trail for details). The files are named DataVSOP87A_Full_Xxx (where 'Xxx' stands for the name of a planet).

    Truncated version

    Having lower precision routines is interesting when large amounts of computations need to be done.
    So a question arises : how to truncate the series in order to keep only the terms useful to get a given precision? An answer is given in A&A 202 p 314 : they say it's possible to have a precision of , where :
    - h is a number smaller than 2 (XEphem uses 2) ;
    - n is the number of terms retained ;
    - A is the amplitude of the smallest term retained retained.
    I first found this formula in XEphem comments, but I didn't know if this could be applied to cartesian coordinates. I found the A&A article later. So I tried to answer to the question. This is detailed in the next page, Truncating VSOP87.
    As we'll see in Testing VSOP87, there is an error in my home made truncation : I wanted a precision of 1 arc second, and got a precision of 4 arc seconds. For the moment, I kept this version, because I couldn't get coherent results from the formula found in A&A.

    This time, it was possible to generate java classes which don't exceed 64 Ko. BuildVSOP87.java was also used to generate them. The generated classes are jephem.astro.solarsystem.vsop87.DataVSOP87A_JEphem_Xxx (Xxx designate a planet name) ; the suffix "_JEphem" indicates that the class was generated using JEphem truncation method.

    I made tests to compare execution time using these java classes and binary files of the truncated version. It appeared to be equivalent. So I kept the java classes for the truncated version, as they can be used by applets.

    Computing the summation

    The class which performs the summation is jephem.astro.solarsystem.vsop87.VSOP87 ; the method to call is calcCoord().
    I wrote VSOP87.java from the fortran code found in example.f in BDL's FTP site (XEphem code was also useful to undersand the fortran code).

    This did not present major difficulties.
    There is still a point that I did not understand : how precision control works (variables prec, p, q in the code). The routine is called with a parameter, precision. During the summation, a test is done to know if a given term should be taken into account. This is the test I did not understand.
    VSOP87.calcCoord() takes a parameter, precision. For the moment, if precision < 4 arc seconds, full version is used, otherwise truncated version.

    Comparison of full and truncated versions

  • Computation time

  • Here is the time taken to compute one planetary position ; (average time taken on 100000 computations for Mercury).

    FullTruncated
    Average computation time
    (PC 1,6 GHz, 256Mo RAM)
    3.4 ms0.09 ms

    Computation is 38 times faster for the truncated version.

  • Size


  • FullTruncated
    Total size (Binary format)923 Ko70 Ko

  • Number of terms

  • VSOP87A_Full
    (BDL)
    VSOP87A
    (JEphem)
    VSOP87D
    (XEphem)
    Mercury6359135277
    Venus2357378169
    Earth3538234406
    Mars7073793691
    Jupiter44343601056
    Saturn75124491833
    Uranus52894581380
    Neptune2636123563
    Total3919829306375