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 5303->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 (using each coefficients file in turn) 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 now ~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? I can explore this more unless you have an idea.
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 two 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 without applying the solution I describe below.
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 second column is the residual of the binary to its ASCII version (they are identical, as expected). The third column is the residual of the SPICE kernel to the ASCII file. The rightmost column is the relative size of the residual from the SPICE file to the ASCII file.
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 of SPICE to the ASCII are consistently relative to the base value on order of 1e-16. 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. The SPICE file appears to have a higher precision than the ASCII.
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 lower than 40 bits 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.
A pull request to address https://github.com/matthewholman/assist/issues/106 , adding support for the spice kernel version of DE440.
Overview of Changes
assist_spk_calc_planetsfunctions for the planets when the spk planets file is detected, otherwise using the existingassist_jpl_calc.assist_spk_calc_planetsdiffers fromassist_spk_calcmainly 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.assist_get_constantfunction that will conditionally get the constant from whichever of the two files the data was initialized from.truncate_doublefunction, 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 5303->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 (using each coefficients file in turn) 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 now ~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:
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? I can explore this more unless you have an idea.
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 two 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 without applying the solution I describe below.
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 second column is the residual of the binary to its ASCII version (they are identical, as expected). The third column is the residual of the SPICE kernel to the ASCII file. The rightmost column is the relative size of the residual from the SPICE file to the ASCII file.
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
jplephemto 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 of SPICE to the ASCII are consistently relative to the base value on order of 1e-16. 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. The SPICE file appears to have a higher precision than the ASCII.
Amelioration
I was able to get identical behavior from both implementations by using a
truncate_doublefunction. 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 lower than 40 bits 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.