Add support for DE440.bsp

[akoumjian]

A pull request to address matthewholman#106 , adding support for the spice kernel version of DE440.

Overview of Changes

• The ephem struct now contains two new optional pointers to the spk_s struct pointed to the planets file and a new spk_global struct that contains the constants and masses loaded from DE440.bsp
• When calculating ephemeris, forces.c will conditionally use the assist_spk_calc_planets functions for the planets when the spk planets file is detected, otherwise using the existing assist_jpl_calc.
• assist_spk_calc_planets differs from assist_spk_calc mainly in that it calculates velocity and acceleration. It also maps the assist planet indices to their spk target codes, rather than directly using an index.
• When using DE440.bsp, there are new functions for joining the mass data to the targets in each file
• There is a new assist_get_constant function that will conditionally get the constant from whichever of the two files the data was initialized from.
• Introducing a truncate_double function, to truncate the significand. This is explained in "Differences in Coefficient Values between ASCII and SPICE", below.

Verification

I have added unit tests for initializing the spk files, loading the constant data, and joining masses to the targets from the loaded constant data.

I have also added a couple of identical unit tests for Apophis and Ceres to make sure tolerances were kept between the two implementations. I manually checked most of the other unit tests as well, but didn't want to duplicate each one until I received opinions on how to structure the unit tests. I could update them all to be parameterized from command line arguments if that seems better.

Earth & Moon Residuals

There is still a disagreement on order 1e-15 AU for the Earth and Moon ephemeris. This difference causes the SPICE and ASCII implementations to perform differently. In particular with the Apophis unit test, the SPICE implementation has a delta from JPL small body of ~47m, whereas the ASCII implementation is ~93m.

The ASCII file contains data for Moon (Earth geocenter frame) and EMB (SSB frame). The SPICE file contains data for Moon (EMB frame), Earth (EMB frame) and EMB (SSB frame). As far as I can tell, this means I cannot use an identical implementation, as the lunar frames are different. For the SPICE implementation, I believe I should be able to simply add the vectors (e.g. Moon (SSB) = EMB (SSB) + Moon (EMB) ).

I have verified that the following values are identical and not the source of the disagreement:

• Earth and Moon masses
• EMB (SSB) ephemeris
• Earth-Moon mass ratio

I'd love suggestions on how to get these number into total agreement. In the case of the ASCII file, Is it possible that using the Earth-Moon mass ratio to calculate the Moon (SSB) and Earth (SSB) vectors is not as accurate?

Differences in Coefficient Values between ASCII and SPICE

This is a bit of a separate topic, and potentially a big one. After a lot of digging, I believe that there is a difference in the double values loaded directly from the files into the memory maps. These differences are quite small and likely below the usable precision of the coefficients. However, the residuals are enough to affect the results.

Here is a table of coefficients for the Sun at J2000. The left column is the value extracted from the plaintext readable ASCII version. The middle column is the residual of the binary to its ASCII version (they are identical, as expected). The right column is the residual of the SPICE kernel to the ASCII file.

| ascp01950.440                  | linux_p1550p2650.440         | de440.bsp                    |
|--------------------------------|------------------------------|------------------------------|
| -1.06808883301021694e+06       | 0.0                          | 2.328306e-10                 |
| 6.43167333734120621e+03        | 0.0                          | 0.0                          |
| 2.01221489458825502e+01        | 0.0                          | 3.552714e-15                 |
| -5.05048275205624633e-02       | 0.0                          | 0.0                          |
| 3.66293150750863095e-03        | 0.0                          | 4.336809e-19                 |
| 6.87050023573884547e-05        | 0.0                          | 0.0                          |
| 1.04739050584889000e-07        | 0.0                          | -3.970467e-23                |
| -9.62554122375419690e-08       | 0.0                          | 0.0                          |
| 1.41137100915376000e-08        | 0.0                          | -4.963084e-24                |
| -1.97668373379263398e-09       | 0.0                          | 0.0                          |
| 5.86673164381741959e-11        | 0.0                          | 0.0                          |
| -3.95512524278347497e+05       | 0.0                          | 0.0                          |
| -8.09283375554062332e+03       | 0.0                          | 0.0                          |
| 1.80149861187318088e+01        | 0.0                          | -3.552714e-15                |
| -8.42604441851760450e-02       | 0.0                          | 0.0                          |
| 1.13093369546329297e-04        | 0.0                          | 2.710505e-20                 |
| 5.03867718918667103e-05        | 0.0                          | 0.0                          |
| 5.93445081607049939e-06        | 0.0                          | 0.0                          |
| -9.87525623150012930e-08       | 0.0                          | 0.0                          |
| 2.33800156279631085e-08        | 0.0                          | 0.0                          |
| 1.59321012060578690e-10        | 0.0                          | 0.0                          |
| 1.56840618138031596e-11        | 0.0                          | 3.231174e-27                 |
| -1.37831006431904796e+05       | 0.0                          | -2.910383e-11                |
| -3.63161643242317587e+03       | 0.0                          | -4.547474e-13                |
| 7.26636726809482969e+00        | 0.0                          | 0.0                          |
| -4.21481400983880414e-02       | 0.0                          | 0.0                          |
| -9.53718391925790533e-05       | 0.0                          | 0.0                          |
| 1.45128782935172393e-05        | 0.0                          | 0.0                          |
| 3.09572793936219392e-06        | 0.0                          | 4.235165e-22                 |
| -4.05084185115523090e-08       | 0.0                          | 0.0                          |
| 1.10153616306458299e-08        | 0.0                          | -1.654361e-24                |
| 2.91662353472744123e-10        | 0.0                          | -5.169879e-26                |
| 1.67388559412173401e-12        | 0.0                          | -4.038968e-28                |

I wanted to make sure that the disagreement didn't have to do with the way that ASSIST was memory mapping the files. I used Brandon Rhode's jplephem to separately load the coefficient values from the SPICE file and it is reading them identically to the ASSIST implementation.

Precision

You may notice that the residuals are proportional the value itself. This could be explained by one or two units of precision difference for the float significand in the two files. I couldn't locate how precise they ought to be (e.g., are the real values using the full precision offered by the double). Perhaps someone more familiar with the specification knows.

Amelioration

I was able to get identical behavior from both implementations by using a truncate_double function. This function masks the bits to a desired precision of the float significand. After manually exploring different levels of truncation, I found that masking down to 40 bits eliminated the difference between the coefficient values. Interestingly, adding this truncation also appears to put the unit tests in closer agreement with the JPL Small Body Code, at least for the Apophis and Ceres unit tests. Truncating further causes them to agree less with JPL Small Body Code, as one might expect.

While the truncation offers a solution for agreement between SPICE and ASCII, I am concerned and curious about why it makes them agree more with Horizons and Small Body. Does this imply that small body, for example, are using truncated or different values than are located in the spice kernels? I don't know, but I want to find out.

[akoumjian]

@moeyensj I'd love to get your eyes on this before I open a pull request against the real assist repo.