Are there any sources with high-confidence state vectors for GRAIL-A?

Question

I have been experimenting with propagation of Moon-orbiting satellites. In particular, of GRAIL-A.

In order to obtain state vectors to use as starting points for my propagation, and to compare the results of it, I used the ephemeris available from JPL Horizons. I have been evaluating my propagation results by calculating the distance at every output time to the ephemerides provided by JPL Horizons at the same times.

However, it has come to my attention that the "accuracy" of my propagation (to be more specific, the distance to the positions returned by JPL Horizons) seems to be very sensitive to the specific epoch from which the state vector is taken, even when comparing results for the same time period.

I focused on the period from 2012-Apr-04 00:29:00 UTC and 2012-Apr-04 10:00:00 UTC, i.e., a period of approximately 10 hours. I did the following experiment:

• Take 30 different state vectors from the ephemerides output by JPL Horizons, at time intervals separated by 1 minute, from 2012-Apr-04 00:00:00 UTC to 2012-Apr-04 00:29:00 UTC.
• Use these 30 different starting points to propagate the trajectory independently in each case until 2012-Apr-04 10:00:00 UTC.
• We therefore obtain 30 different propagated trajectories, some starting earlier than others. To homogenise the compared ranges, I extracted from all of them only the results from 2012-Apr-04 00:29:00 UTC to 2012-Apr-04 10:00:00 UTC. During propagation, I outputted position and velocity every 1 minute.
• For each one, I calculated the distance at each time point to the ephemerides provided by JPL Horizons, and plotted it against time.

The following plot shows the comparison of such distances for the 30 different propagations. The X axis indicates the time in hours from 2012-Apr-04 00:29:00 UTC. The Y axis indicates distance in meters to the corresponding position reported by JPL Horizons. Each line corresponds to a trajectory calculated from a different starting point, whose epochs are listed in the legend.

As it can be seen, the errors differ quite a lot simply by changing the epoch at which the initial state vector for propagation is chosen!

From private discussions, it has been suggested that this could be due to the fact that likely JPL Horizons is storing and providing ephemerides for GRAIL-A and other spaceships in the form of Chebyshev polynomials and doing interpolation for the specific time points requested by users, or some other similar interpolation approach, in a way similar to how celestial bodies ephemerides are provided in JPL DE ephemerides.

Since I am currently interested in exploring the effect of other more subtle forces (such as radiation pressure) in Moon-orbiting spacecraft, this is a bit problematic, since it adds a source of error that overcomes the error due to not modelling such forces.

Therefore, I would like to ask if there is any available source for accurate, interpolation-free state vectors for GRAIL-A? I guess an alternative would be finding out the specific epochs at which the different interpolation intervals begin/end for GRAIL-A in JPL Horizons, since state vectors at those times would be as free as possible of interpolation error, but I have not been able to figure out a way to find out when these would be.

Answer 1

What you want is the SPICE Kernels (SPK) from https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/spk/ . They are named like sssttaaYYYY_DDD_yyyy_ddd.spk, where

sss     3-character spacecraft identifier
            GRA GRAIL-A
            GRB GRAIL-B
            GRX both
tt      Target ID, e.g., LU = Moon
aa      Activity/Experiment ID, e.g. GF = gravity field
YYYY    start year
DDD     start day of year
yyyy    end year
ddd     end day of year

All of GRAIL's mission data products are located somewhere under https://pds-geosciences.wustl.edu/missions/grail/default.htm , and in fact there are two different sets of kernels. The document that told me which SPK to use, and that this set is as good as it gets, is the Grail Data Product Software Interface Specification (pdf), which says (on page 17), "The LGRS RDR data set contains Level 2 products resulting from analysis of the GRAIL science data. The SPK products in this data set differ from those archived by GRAIL navigation; they are created by the GRAIL SDS and make use of the LGRS to provide a more refined solution than those produced by GRAIL Navigation." LGRS stands for Lunar Gravity and Ranging System, the name of the science instrument on the GRAIL spacecraft. The less-refined kernels are the ones archived by NAIF, so don't click the link with that name. Instead, you want the one I listed above, which falls under "LGRS RDR".

EDR, CDR, and RDR, in the names of the Grail data sets, indicate the results of the main steps of the processing chain, from least processed to most processed. The Engineering Data Records are the rawest data product, or "Level Zero" data. These data are exactly what the spacecraft sent back, before anyone did anything to them; this means they are essentially incomprehensible except by the people who designed and operated the spacecraft. The Calibrated Data Records are the first stage of making sense of the return, so they're known as "Level One" data. The procedures used to compute L1 from L0 are described in the Grail Level One Algorithm Theoretical Basis Document (pdf), which sadly has very little to say about orbit determination (Appendix E, 3 pages). The next stage of processing makes "Level Two" Reduced Data Records, also called derived products, which represent the result of understanding and interpreting the raw data the spacecraft sent. The RDRs are "reduced" in the sense of made smaller, but really the term designates "the transformation of... information derived empirically or experimentally into a corrected, ordered, and simplified form". There ought to be a corresponding L2 ATBD, but I haven't found it. A paper that may help fill the gap is Asmar, et al. (2013) "The Scientific Measurement System of the Gravity Recovery and Interior Laboratory (GRAIL) Mission", Space Science Reviews 178(1):25-55.

Answer 2

Partial answer too long for several comments:

tl;dr: The drift you are seeing is just as likely (if not more) to be a difference in your acceleration models as it is an "inaccuracy" (whatever that means) of the starting vectors.

Beautiful analysis! Yes these end-user ephemerides for spacecraft are generated in a very high precision sausage-making process; the ingredients are delay-doppler radio and other observations, an orbital propagator is run over and over from slightly different starting conditions and gravity models until it converges on a solution that can reproduce those observations. The ODE solver chooses it's own time steps, those are then fit with a specific kind of Chebyshev polynomial and those interpolation coefficients are what's stored and retrieved each time someone queries Horizons.

There are even "stitching errors" for some trajectories:

• Is this what station keeping maneuvers look like, or just glitches in data? (SOHO via Horizons)

(note: I think that can see such a glitch in all your traces just after t=7.5) and different solvers and/or interpolators might be used for different periods of time:

• Why do the LRO's orbital elements appear to constantly oscillate - except recently?

Therefore, I would like to ask if there is any available source for accurate, interpolation-free state vectors for GRAIL-A? I guess an alternative would be finding out the specific epochs at which the different interpolation intervals begin/end for GRAIL-A in JPL Horizons, since state vectors at those times would be as free as possible of interpolation error, but I have not been able to figure out a way to find out when these would be.

There will be something close to what you are looking for somewhere in the scientific literature and websites associated with the mission. Hopefully other answers will be able to cite specific links.

But if you want some state vectors that are the closes to being "right" (whatever that means) within a specific gravitational model, you'll need the output from the ODE solver before interpolation, AND the gravitational model and all other assumptions included (solar photon pressure, gravitational effects of the Earth,. Sun, Venus, Jupiter, etc.) before your solver will track another solver.

The drift you are seeing is just as likely (if not more!) to be a difference in your acceleration models as it is an "inaccuracy" (whatever that means) of the starting vectors.

Answer 3

Even though @RyanC and @uhoh already provided excellent answers with all the details I was looking for, I thought it might be useful for future readers to know the outcome of using SPK kernels.

To summarize, I used the SPK kernels available at https://pds-geosciences.wustl.edu/grail/grail-l-lgrs-5-rdr-v1/grail_1001/spk/ , as suggested by @RyanC.

These kernels are of type 1, meaning they provide Modified Difference Arrays (MDAs). How to evaluate these at any arbitrary target time is another story, but the kernels also include reference positions and velocities at a series of epochs which can be used directly and are as interpolation-free as it gets (at least for the general public, I guess).

So I took all those reference positions and velocities, together with their epochs, to cover the same time interval I had been working with before (10 hours after 2012-Apr-04 00:00:00 UTC). I then propagated trajectory using as starting point the 1st reference position, velocity and corresponding epoch. I outputted propagation results at all the epochs of other available reference positions found in the SPK kernels, and calculated the distance to the reference position at each time point. Propagation was performed exactly as when I was using data obtained from JPL Horizons app. A summary of propagation details:

• Propagation done in Moon-centered ICRF frame
• Using JPL DE 440 for positions, velocities and accelerations of Sun, planets and Moon
• Forces modeled: Moon gravity field with GRGM1200B model, Earth gravity field with GGM05C, gravity by Sun and planets modeled as point masses.

The results were indeed way better than using starting points taken from JPL Horizons, with an error of around 6 meters after 10 hours propagation. This is a plot for the distance to the reference positions found in SPK kernels during the 10 hours propagation (X axis is time in hours, Y axis is distance in meters).

Still not perfect, though, so need to look for other error sources/improve force model.

But in summary, it seems indeed that using these interpolation-free starting points seems key to obtaining more accurate propagation results (well, I guess to be precise, key to reproducing better the trajectory available in SPK kernels).

Another interesting observation is that when using these reference positions, the strange discontinuities observed previously and pointed out by @uhoh are not present anymore. Strange, since I would have thought JPL Horizons uses these SPK kernels internally.