You are right, there are no way to use custom Love numbers, and this would be a great improvement!
Some more findings,
I tried to add some additional force models (relativity, Moon albedo, Moon solid tides k3, modelling GRAIL-A with solar array + box model, adding spherical harmonics for Earth up to 4x4), and most of these either had a very modest effect compared to the still remaining position error, or improved error modestly at some times during the propagation run but made it worse during others (which to me seemed to somehow suggest it was not addressing really the underlying cause for error).
I doubled propagation time to ~3.5 hours and then noticed the total position error seemed to have 2 components, one always growing at a constant rate and another oscillatory one. So for the constant growing component I wondered if it still would be due to some small error in either Moon gravity model or small frame misalignment. I tried to use different GM value for Moon (so far I was always using the value provided with GRGM1200B), trying also the DE430 and DE421 values. Interestingly, both of these made a modest improvement in propagation error, seemingly reducing the constantly growing error.
So I decided to try a sweep of slightly modified GM values (up to around 50 ppb smaller than GRGM1200B’s Moon GM), and the results showed that it’s possible to reduce the position error greatly in this way, see the below plot:
Then I wondered if reducing GM value just was an overfit to this specific part of the trajectory or really made things improve, so I tried propagation for ~10 hours flight time with the same setup. And for this run I’ve also plotted radial, in-track and cross-track components of position error:
This was quite interesting:
- It seems that indeed using the lower GM value was not just an overfit to that part of the trajectory, since it maintains a similar total position error when going from 3.4 hours to 10 hours propagation
- The continuously increasing error component seems to be in-track, and still seems like there is a bit left (blue line), which maybe could be completely removed by more finely “tuning” the GM value
- The periodic error is cross-track, and the period is just under 2 hours, which seems to match closely with GRAIL-A orbit period. So probably related to SRP?
Of course, I think fine-tuning GM value is not good 
But doing so must somehow be compensating for some error during the setup likely related to Moon’s gravity field. Maybe some frame misalignment, small error in Moon’s position relative to Earth (since this transformation is used to construct the Moon-centered ICRF used during propagation) or some small timescale conversion error?
I wonder if that could not be related to relativity. In Orekit, we work internally in TAI, and we often use UTC for input/output. Both of these time scales are based on the flow of time at Earth surface, which depends on the local gravity. As we go further from Earth, time flow changes (see Relativistic Clock Correction in navipedia). This is the reason the base frequency of the on-board clock in MEO GNSS satellites is very slightly below 10.23MHz so the signal on Earth is based on 10.23MHz. I think this is also the reason we use GM=398600.4415km³/s² for Earth near the body but the value for interplanetary computation in DE files is rather GM=398600.4329km³/s².
I know that for future GNSS constellations around the Moon, there is one specific problem which is to set up a proper Moon time scale. You can see for example GNSS-based Lunar Orbit and Clock Estimation With Stochastic Cloning UD Filter section C.1. We should perhaps add TCL (Lunar Coordinate Time) into Orekit? If so, I guess we could still work in TAI but add a relativistic correction for Moon orbits.
Thanks a lot @luc for pointing this out. This indeed sounds very likely to be the reason!
I’m unfortunately quite unfamiliar with relativity corrections, so it’s taking me a while to read into it and understand the details. However, I found this article, which seems quite relevant. From it, it seems to be that in order to define a proper Moon-centered frame for propagation with high precision, it is not enough to just take the GCRF frame and apply a translation to center it around the Moon, and then a rotation to make it be aligned with Moon’s principal axes, which is exactly what I’ve been doing so far using SPK kernels/DE430. If I understand correctly, GCRF frame is an Earth-centered implementation of BCRS system, but relativistic corrections will be different in the proximity of the Moon vs in the proximity of Earth.
From the paper, it seems equation 63 gives an accurate perturbation of acceleration due to relativistic corrections. So I’ve tried to add the first two terms in it as additional perturbing accelerations. On the other hand, equation A24 gives a position correction from BCRS into a proper LCRS. Equations A25-A27 give transforms for position/velocity/acceleration, but these require a scaling factor, and I couldn’t see a scaling transformation incorporated into frame transforms. So instead, I tried to use equation A24 to convert at each propagation step the position into a proper LCRS position, then calculate Moon’s attraction there as a point mass there and also at the uncorrected position, and add this difference to Moon’s attraction previously calculated through full spherical harmonics. My idea in doing this correction through these 2 paths was that if it’s a step in the good direction, they should both lead to similar effects.
And also, since it would be complicated to separate other effects such as inaccurate modelling of SRP etc, I searched in the SPK kernels for GRAIL-A for arcs where GRAIL-A would be in full eclipse of Sun, so SRP should not matter. I took 10 of these full eclipse arcs, and then I’ve done 3 propagations for each arc:
- With the best setup so far but without fine-tuning GM at all (i.e., using GRGM1200B Moon’s GM as-is) (red lines)
- Adding the correction as an additional acceleration perturbation from the first 2 terms of eq 63 of the linked article (blue lines)
- Adding the correction through eq A24 as described (green lines)
See the results here for the 10 arcs:
It seems like interesting results, and looks like a step in the good direction probably to me!
However I’m still not convinced because the amount of improvement varies quite a bit between arcs, although I also checked the residual acceleration by components (radial, along track and cross track) for the 10 arcs and the radial component was similarly corrected in all cases, with the arcs where there is less improvement having, from the beginning, a larger acceleration error in cross-track components.
I’ll continue digging into this but wanted to share these preliminary findings in case someone with proper expertise can figure it out properly!
This looks great!
Indeed, scaling factors are not present in transforms as we consider them to be isometries.
This already implied we departed from official definitions in Earth frame, for the Helmert transforms between one ITRF realization and another. These transforms do include a tiny scaling factor and we ignore it.
I asked one of the people working on these effects (Xavier Collilieux, from IGN) explaining him I had concerns with the scaling factor when applied to satellites and not to points near Earth. He confirmed to me these scaling factors were indeed computed to reduce error on surface reference points and not for points at high altitudes. He therefore confirmed to me that the approach to ignore the scaling was not too bad, but of course this was not really studied, it was just a rough guess.
Some more plots, I’ve extracted all eclipse arcs of GRAIL-A, for a total of 1780. I’ve sorted them by decreasing mean altitude and want to propagate all of them to see how the error evolves, but for a quick test I took 60 sampled homogeneously and propagated them with the correction of first 2 terms of the EIH equation:
I notice the following:
- For very low altitude arcs (with minimum altitude below ~15 km), there is mostly no effect, and the position error is to begin with larger. Maybe because higher SPH degree is needed, I’ve used 600x600 so far, but will have to test higher I guess.
- Leaving those aside, the correction seems to work quite generally, in many cases bringing the error to below 5 cm across different altitudes!
- For cases where the correction has less of an impact in the total final error, they are indeed almost always dominated by cross-track position error.
I’ll try to see if I can get some more insights from the full set of eclipse arcs
This is impressive.
I guess you could write a paper on this!
Hi all,
I have continued the analysis done by @Rafa on Orekit’s propagation comparing it against the LRO reconstructed SPK. I am not using OD residuals here and I am not estimating range/range rate biases. The goal was only to understand how much of the SPK propagation mismatch can be reduced by improving the dynamical model and the Moon-centered relativistic/time-scale handling.
The validation set used for the main numbers below is seven LRO ES arcs, propagated for about two orbital periods each and sampled every 300 s. I also kept a 12 h, nine-arc check so that I do not tune the model only on one short arc.
1. Reference Propagation Model Without Lunar Solid Tides
The first reference model is Moon-centered and propagation-only, with lunar
solid tides deliberately disabled. It contains:
- GRGM900C lunar gravity at 360/360;
- Moon GM read from the GRGM900C field:
4.90279996708864e12 m^3/s^2; - lunar harmonics mounted on the Moon PA frame from the DE440 binary PCK;
- Earth, Sun and Jupiter point-mass third bodies;
- Orekit central-body post-Newtonian term with
new Relativity(MoonMu); - simple LRO box plus solar-array SRP model, including Moon and Earth
occultation; - simple LOF attitude for the plate model;
- no lunar solid tides;
- numerical propagation in the Moon-centered frame and comparison to the
reconstructed SPK states transformed to the same frame.
The lunar PA frame is possible here thanks to Rafa’s SPICE reader work. The Moon PA frame is built by parsing the DE440 binary PCK directly and using it as the transform provider for the rotating lunar frame:
File pckFile = resolveOrekitDataFile("MoonPole/moon_pa_de440_200625.bpc");
DataSource moonPADataSource =
new DataSource("moon_pa", () -> pckFile.toURI().toURL().openStream());
PCK moonPApck = new PCKParser().parse(moonPADataSource);
Frame moonCentered = new Frame(FramesFactory.getGCRF(),
new LME2000PositionProvider(FramesFactory.getGCRF()),
"moonCentered", true);
Frame moonPA = new Frame(moonCentered,
new MoonPATransformProvider(moonPApck),
"MoonPA");
The corresponding LRO force-model setup is:
perturbingForces.add(configureLunarGravityField(lunarGravityDegree,
lunarGravityOrder));
perturbingForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getEarth()));
perturbingForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
perturbingForces.add(new ThirdBodyAttraction(CelestialBodyFactory.getJupiter()));
perturbingForces.add(new Relativity(getMoonMu()));
// No lunar solid tides in this first reference.
This first reference gives:
No lunar solid tides, seven two-orbit arcs, 300 s sampling:
overall RMS = 4.338637 m
Per-arc RMS for this no-tides reference is:
| Arc | RMS [m] |
|---|---|
| LRO_ES_19 | 5.651338 |
| LRO_ES_39 | 1.379472 |
| LRO_ES_67 | 4.278672 |
| LRO_ES_72 | 2.627926 |
| LRO_ES_78 | 2.988822 |
| LRO_ES_115 | 5.734987 |
| LRO_ES_120 | 5.533215 |
2. Solid Lunar Tides
Rafa’s issue-1925 branch adds a SolidTides constructor accepting custom LoveNumbers, which avoids using IERS/Earth Love numbers for a non-Earth central body:
Just to do a cross check I have implemented locally a wrapper in the same package as Orekit’s
SolidTidesField. But the force model being built is equivalent to Rafa’s new constructor for the terms used here:
final SolidTidesField raw =
new SolidTidesField(loveNumbers, zeroFrequencyDependence,
0.0, null,
centralBodyFrame, ae, mu, centralTideSystem, bodies);
final NormalizedSphericalHarmonicsProvider provider =
nbPoints < 2 ? raw :
new CachedNormalizedSphericalHarmonicsProvider(raw, step, nbPoints,
OrekitConfiguration.getCacheSlotsNumber(),
7 * Constants.JULIAN_DAY,
0.5 * Constants.JULIAN_DAY);
ForceModel tides =
new HolmesFeatherstoneAttractionModel(centralBodyFrame, provider);
I also checked this numerically. I copied Rafa’s issue-1925 constructor/build path into a test-only helper and compared it against the local wrapper using the same Moon PA frame, same lunar Love numbers, same Earth/Sun tide-generating bodies. The test sampled 18 actual LRO SPK states from LRO_ES_19, LRO_ES_72, and LRO_ES_115.
solid-tide acceleration max norm: 4.370044278474152e-07 m/s^2
local wrapper vs Rafa-style max acceleration: 0.000000000000000e+00 m/s^2
2 h tide-only propagation position difference: 0.000000000000000e+00 m
2 h tide-only propagation velocity difference: 0.000000000000000e+00 m/s
So for the frequency-independent custom Love number model used here, the local
wrapper and Rafa’s proposed constructor give identical numerical results.
The lunar Love numbers currently used by default are order-specific degree-2
values, plus a degree-3 term:
array[2][0] = 0.02408; // k20
array[2][1] = 0.02414; // k21
array[2][2] = 0.02394; // k22
fillDegreeFromProperties(real, 3, 0.0163); // k3
fillDegreeFromProperties(real, 4, 0.0); // k4 disabled by default
The tide-generating bodies are Earth and Sun by default.
The direct no-tides to tides comparison is large. This table uses the same propagation setup as the first reference and toggles the lunar solid tide model:
| Model step | Overall RMS [m] |
|---|---|
| Reference: no lunar solid tides | 4.338637 |
Same reference + custom lunar k20/k21/k22, k3, Earth+Sun tide bodies |
0.808748 |
I also compared this tides implementation against one that uses only the scalar value k2. That reference contains GRGM900C 270/270, Earth/Sun/Jupiter point masses, central lunar relativity, simple SRP/LOF attitude, andGRAIL/TCL state/time mapping:
| Tide model | Overall RMS [m] |
|---|---|
Scalar lunar k2 = 0.024059 |
0.842168 |
Order-specific k20/k21/k22, k3, Earth+Sun tide bodies |
0.831574 |
I also swept the degree-4 lunar term available through this Orekit path:
| k4 | Overall RMS [m] |
|---|---|
| -0.0040 | 0.808746 |
| 0.0000 | 0.808748 |
| 0.0020 | 0.808749 |
| 0.0050 | 0.808751 |
| 0.0100 | 0.808754 |
| 0.0163 | 0.808757 |
So in this validation the degree-4 term is neutral; I just dropped it since it didn’t affect the final results by measurable margins.
3. Relativity and Lunar Time-Scale / Coordinate-Scale Handling
The central lunar relativistic force is still the standard Orekit force model:
perturbingForces.add(new Relativity(getMoonMu()));
I did not try to replace this term with a custom acceleration. Instead, based on Rafa’s comments and the GRAIL/TCL discussion, I treated the comparison against SPK kernels as a second issue: the local Moon-centered coordinates should not be only “GCRF translated to the Moon plus a rotation”.
The first-order mapping currently tested is:
r_BCRS ~= r_TCL * (1 - W_ext(Moon) / c^2)
- 0.5 * v_Moon * (v_Moon . r_TCL) / c^2
dt_TCL/dt ~= 1 - (0.5 * |v_Moon|^2 + W_ext(Moon)) / c^2
where W_ext(Moon) includes the external potential at the Moon from Earth, Sun and Jupiter.
The implementation in the validation runner maps the initial SPK state before propagation and maps propagated states back for comparison:
final double externalPotential = externalPotentialAtMoon(date, moonPV);
final double c2 = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT;
final Vector3D positionGcrf = pvGcrf.getPosition()
.scalarMultiply(1.0 - externalPotential / c2)
.subtract(moonVelocity.scalarMultiply(
0.5 * moonVelocity.dotProduct(pvGcrf.getPosition()) / c2));
final Vector3D velocityGcrf = pvGcrf.getVelocity()
.scalarMultiply(1.0 - externalPotential / c2)
.subtract(moonVelocity.scalarMultiply(
0.5 * moonVelocity.dotProduct(pvGcrf.getVelocity()) / c2))
.scalarMultiply(localCoordinateTimeRate(date));
To keep this separate from the force-model changes, I used a second reference for this section: the same propagation force model with and without the GRAIL/TCL state and elapsed-time mapping. This reference contains GRGM900C 270/270, scalar lunar k2, Earth/Sun/Jupiter point masses, central lunar relativity, simple SRP/LOF attitude, and the same SPK comparison procedure. I also tested two extra variants: adding the small acceleration-dependent position terms from the practical transform and adding the endpoint v_Moon.r/c^2 time term. In this SPK comparison those additional terms were neutral/slightly worse:
time-scale reference force model 1.363802 m
GRAIL/TCL state and elapsed-time mapping 0.842168 m
same + acceleration-dependent position term 0.842173 m
same + endpoint v_Moon.r/c^2 time term 0.842174 m
The two-orbit relativity comparison is:
Over longer 12 h arcs the same correction is not uniformly better on every arc,
but the aggregate improves. I also computed the reconstructed-SPK altitude over
the Moon PA ellipsoid for each of these same segments, using the same sample
epochs as the RMS computation:
| Arc | Mean alt [km] | Alt range [km] | Reference RMS [m] | GRAIL/TCL RMS [m] |
|---|---|---|---|---|
| LRO_ES_10 | 110.6 | 48.9-170.5 | 5.322400 | 9.116917 |
| LRO_ES_19 | 109.4 | 55.4-166.2 | 0.554166 | 4.360724 |
| LRO_ES_39 | 98.5 | 25.7-165.8 | 6.233978 | 2.431299 |
| LRO_ES_67 | 93.4 | 67.1-123.0 | 3.406179 | 0.977517 |
| LRO_ES_72 | 95.2 | 37.5-149.1 | 5.914093 | 2.120016 |
| LRO_ES_78 | 93.7 | 63.3-127.4 | 4.367737 | 0.747982 |
| LRO_ES_95 | 95.5 | 66.1-125.7 | 5.738411 | 2.171670 |
| LRO_ES_115 | 93.1 | 78.3-107.1 | 4.889443 | 2.600571 |
| LRO_ES_120 | 93.5 | 84.6-104.8 | 1.680683 | 2.735076 |
On these nine 12 h segments, the complete-model RMS is positively correlated with mean altitude (r ~= +0.85), and the improvement from reference to complete is negatively correlated with mean altitude (r ~= -0.80). I would not over-interpret this yet because the sample is small and altitude is probably also acting as a proxy for mission epoch, geometry, and reconstructed small forces. Still, it is interesting that the two highest-mean-altitude segments (LRO_ES_10 and LRO_ES_19) are also the cases where the GRAIL/TCL correction worsens the 12 h RMS most strongly. Conversely, LRO_ES_39 reaches the lowest sampled altitude, about 26 km, but improves significantly. So this does not look like a simple “lower altitude means worse propagation” effect.
This mixed behavior is why I am trying not to tune on one arc. ES_19 in particular looks suspiciously along-track dominated and may have something going on (I have no clue about that behaviour).
4. Lunar Gravity and Earth Noncentral Harmonics
After the tide and time-scale work, I retested the lunar gravity truncation:
| Model | Overall RMS [m] |
|---|---|
| Custom lunar tides + GRGM900C 270/270 | 0.844806 |
| Custom lunar tides + GRGM900C 360/360 | 0.828982 |
| Custom lunar tides + GRGM900C 450/450 | 0.828246 |
| Custom lunar tides + GRGM900C 600/600 | 0.828943 |
The selected propagation setup now uses 360/360 because it gave the best accuracy/runtime trade-off in the multi-arc validation.
I also rechecked Earth 2x2 harmonics after fixing an earlier implementation mistake. In a Moon-centered propagation the Earth monopole is already modeled by ThirdBodyAttraction(Earth), so the extra Earth harmonics model should add only the differential noncentral part:
final SpacecraftState originState = new SpacecraftState(
new AbsolutePVCoordinates(state.getFrame(), state.getDate(),
PVCoordinates.ZERO),
state.getMass());
return gravityModel.acceleration(state, gravityParameters)
.subtract(gravityModel.acceleration(originState, gravityParameters));
After this fix, Earth 2x2 is very small in these LRO arcs:
best tides, 600 s screen:
0.844806 m -> 0.844766 m
complete model, 300 s screen:
0.808748 m -> 0.808720 m
So I left out the Earth 2x2 harmonics effect.
5. SRP, Attitude, and Eclipse Checks
I also made an eclipse-focused plot for ES_72, showing the propagated error together with Moon eclipse, Earth eclipse, and the overlap interval. I made this analysis with the objective of understanding if the model, as implemented, is close enough with the LRO SPK kernels deliberatively avoiding one large source of error as the unmodelled SRP as @Rafa has done in its posts .
To me, at least, it seems like there is no large or evident effect of the mis-modelling of the SRP force (most probably it has a visible effects over longer arcs).
6. Current Best Propagation-Only Result
The main current two-orbit comparison against the no-solid-tides reference is:
| Step | Overall RMS [m] |
|---|---|
| Reference: no lunar solid tides | 4.338637 |
Same reference + custom lunar k20/k21/k22, k3, Earth+Sun tide bodies |
0.808748 |
| Same reference + optional Earth 2x2 differential harmonics | 0.808720 |
The mean improvement is real, but the arc-to-arc behavior is not perfectly
uniform. That is why I currently suspect the remaining mismatch is a mixture of:
- incomplete SRP/attitude modeling;
- reconstructed small-force or maneuver details in some SPK segments;
- the formal way the local lunar relativistic coordinate/time transform should
be represented in Orekit.
Any comments are highly appreciated 
.
All the Best,
Davide
References
- Rafa’s custom-Love-numbers branch:
Files · issue-1925 · Rafael Ayala / Orekit · GitLab · Orekit - Orekit forum discussion on Earth/lunar spherical harmonics in Moon-centered
propagation:
Usage of spherical harmonics for both lunar and terrestrial gravity in a Moon-centered propagation - Relativistic cislunar time-scale paper discussed in the thread:
https://iopscience.iop.org/article/10.3847/1538-4357/ae25ff
Hi,
Thanks for this in-depth analysis (which I’ve not looked in details).
About Rafa’s issue on Love numbers, there’s a world where it could be considered a bug and cherry picked for a patch release. What do you think?
Cheers,
Romain
Hi @Serrof,
Yeah It would be correct to classify it as a bug, which has a fairly simple fix 
.
All the best,
D
Thank you @DDega for this very impressive analysis !
I also hadn’t had time to look into it in details yet, but I will.
Maybe we should talk about this in a separate discussion, but there is still @Rafa merge request on the forge on SPK kernels handling that is waiting for a review.
It’s quite large and I do not have time to take care of it right now. Is there someone from @orekit_dev_team who could do it before 14.0 release ?
I mentionned this here because it’s part of a larger work by @Rafa about SPICE kernels handling (see issues and MR with the “SPICE” tag on the forge), and @DDega mentionned PCK reading in his analysis.
Cheers,
Maxime
Hi!
About @Rafa merge request (which represent an important amount of work and thank you very much again for it), I would appreciate a status on it. First, because I see two MR. Are they independent? Which one should be merged first? Secondly, the two MR are on draft status, so technically not mergeable. Is this status still on date?
Thank you
Bryan
Hi @bcazabonne so the overall goal of the SPICE work is to provide support for reading/writing/evaluating all kernel types supported by SPICE! At the beginning I had a huge MR but we decided that to make it feasible to review to split it into parts. So we split the current work into 3 MRs:
- An MR for DAF files (generic file types that provide the basis for other files)
- An MR for SPK files (these are ephemerides files)
- An MR for PCK files (which provide mostly angles for frame rotations)
The first one for DAF files was kindly reviewed by @MaximeJ and @Vincent , and it has been merged. But the two remaining ones are still quite large. We decided to start with the SPK merge request, because it covers a larger range of subtypes, so I thought it would be better to start with that since it would be better to change decision choices etc in it, and then propagate those to PCK rather than the other way around. We decided to leave them in draft status to avoid accidental merge until someone has the chance to properly review!
Actually I just checked and I haven’t opened an MR for PCK files yet (although I have the code in one of my branches), since it seemed better to wait until the SPK MR is done.
I see the confusion, it seems there is currently MR 756 and MR 1099 open! 756 is the original massive one that we decided to split, so that one should be ignored (and I guess it could just be rejected? since it will be now split over the 3 smaller ones)