3

We have discussed how to calculate positions from JPL development ephemerides in this question.

Following this answer, I am now using SPICE's algorithm for calculating the derivatives of the Chebyshev sums, which are in fact the derivatives of the positions, i.e., velocities.

However, for the case of Earth and Moon, as previously discussed, JPL DE in fact provide the Chebyshev coefficients for the Earth-Moon barycenter (EMB), not for the Earth alone, as well as for the Moon itself. For calculating positions, this is not much problem, since we can use the Earth-Moon mass ratio together with the calculated positions of Moon and EMB to calculate the position of the Earth.

However, this leaves me wondering, if we use the derivative of Chebyshev sums, we would in fact get the velocity of the EMB with respect to the Solar System Barycenter, as well as the velocity of the Moon with respect to the Earth. How can we use this information to calculate the velocity of Earth with respect to the SSB?

Rafa
  • 351
  • 2
  • 10

3 Answers3

4

It should be quite obvious that you compute the velocity of the Earth the same way you compute the position. It's multiplying a function by a constant, and that constant will be unaffected by the derivative.

Here is the code I use for computing the Earth position and velocity from the Moon and EMB position and velocity. The variables emb and moon are arrays containing the x,y,z postion and dx,dy,dz velocity (in that order).

    def getEarthPositionFromEMB(emb,moon):
        earthMoonRatio=Decimal(0.813005600000000044E+02)
        earth=[0,0,0,0,0,0]
        for i in range(6):
            earth[i]=emb[i]-moon[i]/(Decimal(1)+earthMoonRatio)
        return earth
Greg Miller
  • 5,857
  • 1
  • 14
  • 31
  • 2
    The full code is here: https://github.com/gmiller123456/jpl-development-ephemeris/blob/master/de.py – Greg Miller Mar 02 '22 at 21:13
  • 1
    Where did you get that Earth/Moon mass ratio? It's slightly low compared to the value I get from dividing the GM values given in Horizons body data files, 81.30056908. BTW, it's better to initialize Decimals with a string or int. If you initialize with a float, the value will contain spurious digits from the decimal to binary conversion performed at float precision (53 bits). See https://docs.python.org/3/library/decimal.html#decimal.Decimal – PM 2Ring Mar 03 '22 at 01:38
  • Horizons body data – PM 2Ring Mar 03 '22 at 01:51
  • Duh, that makes sense... I feel so dumb now! @PM2Ring The value I get for EMRAT by the way is closer to what you get, although still not exactly. I'm using 0.813005682214972154E+02, which is directly taken from file header.440 here. Under section GROUP 1040, we see EMRAT is the 11th value. So if we go to section GROUP 1041, we need to get the 11th value also (that would be 4th row, column 2). That is for JPL DE440. It will probably be slightly different for other DE I guess, and I believe Horizons uses by default JPLDE441 – Rafa Mar 03 '22 at 01:53
  • In fact just checked header.441 and it has the same value. I guess the difference might be to rounding errors from division? – Rafa Mar 03 '22 at 01:56
  • Ok, so it seems @GregMiller's code is using DE405 and DE438. The EMRAT value in DE405's header.405 is indeed 0.813005600000000044E+02, so I guess that's where it comes from – Rafa Mar 03 '22 at 02:06
  • 1
    Horizons uses DE441 & DE440. You can get the Horizons news file by giving a body name of news in that body data script. OTOH, the values given in the body data files aren't from the DE, so your values from the DE file headers are possibly better. ;) FWIW, here's the script I used for the mass ratio calculation. https://sagecell.sagemath.org/?z=eJxFzr0KgzAUhuE9V3GggwoSjya1ZuimdHJyFySNrbRJJEbw8utPacfzvPBxemc13JUcdPeGQY_WeSi_ZzdBGcNDeWmNV4sn5AS3GmwPzWxiqDrnnzHU1hoYDLx0y5KpzYieG7hCGQYpyy5pxjkiT6lgghVBtNbqqEwUOSLl7MxZvof6CFxgRgtEzDcm0i-r_78Io43o6JRcPUVCRjcYH27DCawr0Q_qHaroA_9JQIM=&lang=python – PM 2Ring Mar 03 '22 at 02:08
  • 1
    Yea, that value comes from DE405 which is what the code is written to compute, I guess I should make that a parameter. – Greg Miller Mar 03 '22 at 02:37
  • 1
    The trailing 00000000044 at the end of the DE405 Earth-Moon mass ratio is an artifact of IEEE floating point. The value back then was 81.30056, period. Any more digits after that would not have been real. With multiple missions to the Moon since the 1997 release of DE405, including the two GRAIL satellites, the value now is better stated as 81.300568. All of the extra digits in the DE440 header are superfluous. – David Hammen Mar 03 '22 at 08:51
  • @PM2Ring This means that your comment about converting a float to a Decimal is a bit superfluous as well. The 1997 value had just a bit over IEEE single precision accuracy. The current value remains well within the scope of an IEEE double precision number. – David Hammen Mar 03 '22 at 08:54
  • 1
    @DavidHammen Sure, we don't really need Decimal for this computation, since the values are all well within the 53 bit precision of Python floats (equivalent to C double). And of course the precision of the calculation is limited to that of the least precise value. My earlier comment was mostly just general advice about the proper initialisation of Decimal values. There are quite a few questions on SO where the OP is puzzled by those spurious digits. ;) – PM 2Ring Mar 03 '22 at 09:47
  • @PM2Ring It's interesting to me that JPL switched from double precision to quad precision in their integration with DE405. I suspect there was a relative error on the order of $10^{-13}$ prior to DE405 just from the integration alone due to the use of double precision arithmetic. – David Hammen Mar 03 '22 at 10:31
  • @David Right. As we've discussed previously, higher precision is certainly useful when integrating differential equations, especially if the step size is small. Also, it can be necessary when computing Chebyshev coefficients, as mentioned in the mpmath docs. FWIW, that also applies to computing minimax polynomial coefficients with the Remez algorithm. – PM 2Ring Mar 03 '22 at 10:59
2

The JPL Development Ephemerides provides the Earth-Moon mass ratio to several digits. (A handy and easy number to remember: The Moon's mass is about 0.0123 Earth masses. This is low accuracy, but easy to remember.)

The combination of the vector from the solar system barycenter to the Earth-Moon barycenter, the vector from the Earth to the Moon, the Earth-Moon mass ratio, and (for example) a vector from the solar system barycenter to another point in the solar system enables one to compute the vector from the Earth to that other point in the solar system.

David Hammen
  • 33,900
  • 3
  • 74
  • 125
  • I understand how to use the Earth-Moon ratio (nice mnemonics by the way!) to calculate the position of Earth and from there the Earth-centered position of coordinates in the solar system. However, I am unsure about how to calculate the velocity of Earth. This is because, if I understand correctly, by calculating the derivative of the Chebyshev sum for the EMB, we actually get the velocity of the EMB, not of Earth. Or did I miss something? – Rafa Mar 02 '22 at 18:15
  • In fact this document seems to state that the same relationship that applies to positions (R_earth = R_earthmoon - EMRAT1 * R_moon), also applies directly to velocities (VR_earth = VR_earthmoon - EMRAT1 * VR_moon). Interesting! I thought the calculation for velocity would be somehow different than for positions – Rafa Mar 02 '22 at 18:49
  • 1
    Actually, that value is pretty good, since the next few digits are zeroes. Using Horizons data, I get 0.01230003690 – PM 2Ring Mar 03 '22 at 02:12
  • 2
    @Rafa Unless you want to go full-bore general relativistic, it is as simple as that. JPL's velocities are full-bore general relativistic, so you might need to be careful about time scales. JPL uses its own home-brewed relativistic time scale $T_\text{eph}$, which is essentially the same as Barycentric Dynamic Time (TDB). – David Hammen Mar 03 '22 at 09:11
  • 1
    @PM2Ring That's why I like that value (0.0123). It's so easy to remember and it is good to six places of precision. The Moon-Earth mass ratio is known to about eight places, which is why I wrote "low accuracy". – David Hammen Mar 03 '22 at 11:24
1

If you're using SPICE, you can directly have Earth's Chebyshev coefficients from bsp files. The SPICE "brief" utility allows you to display the contents of an SPK bsp file.

Example with de440.bsp file, available bodies are :

BRIEF -- Version 4.1.0, September 17, 2021 -- Toolkit Version N0067
Summary for: de440.bsp

Bodies: MERCURY BARYCENTER (1)  SATURN BARYCENTER (6)   MERCURY (199)
    VENUS BARYCENTER (2)    URANUS BARYCENTER (7)   VENUS (299)
    EARTH BARYCENTER (3)    NEPTUNE BARYCENTER (8)  MOON (301)
    MARS BARYCENTER (4)     PLUTO BARYCENTER (9)    EARTH (399)
    JUPITER BARYCENTER (5)  SUN (10)
    Start of Interval (ET)              End of Interval (ET)
    -----------------------------       -----------------------------
    1549 DEC 31 00:00:00.000            2650 JAN 25 00:00:00.000

You can then extract Earth's Chebyshev coefficients and compute derivative to get velocities. I wrote a CSPICE algorithm to extract Chebyshev coefficients of a particular body if you are interested.

GuillaumeJ
  • 21
  • 2
  • 1
    Thanks for your answer! Do you have by any chance a source confirming that indeed bsp files contain Chebyshev coefficients for Earth directly, and not for the EMB? Welcome to Astronomy SE! – Rafa Mar 16 '22 at 00:55
  • 1
    Yes it's in the link I gave just above SPK bsp file slide 35. – GuillaumeJ Mar 16 '22 at 08:24
  • 1
    I also checked by myself that we found the position of the Earth directly with these coefficients by comparing with Horizons. Note that each bsp file can contain different information, the "brief" command that I gave above allows you to display the content. – GuillaumeJ Mar 16 '22 at 08:32
  • 1
    Interesting! It seems the binary versions of the ephemerides (.bsp files) contain Chebyshev coefficients for Earth directly as well (this is confusingly documented; I have not found a documentation as good as for the ASCII version), but not for Moon libration angles, for example, which are indeed present in the ASCII versions... – Rafa Mar 16 '22 at 09:14
  • 2
    Yes SPK bsp files are not containing sames datas as ASCII versions. I have found an other documentation of datas in bsp files. For the Moon libration angle I think you can use PCK kernel which gives orientation of the Earth and Moon. – GuillaumeJ Mar 16 '22 at 10:21
  • 1
    Wow! That's really interesting. I guess now the only thing left to do is test and see if it produces the same result as computing the Earth's center from the EMB and Moon positions. – Greg Miller Mar 16 '22 at 13:55