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