Numerical propagation in the fields of Earth and Moon with higher-order spherical harmonics

Hi Marc and Luc,

I’m also interested in simulating Earth-Moon retrograde orbits with Orekit and I tried Luc’s suggestion for the Moon’s harmonics. Please have a look at my attempt, it is based on Orekit’s Numerical Propagation tutorial. Technically, it seems to work, i.e. the HolmesFeatherstoneAttractionModel is added ok to the propagator and doing so results in a noticable change of the final state. In my example of a 60-day simulation of a wanna-be distant retrograde orbit the change in the final state is of ~ 1 km. Whether it works properly regarding the Moon’s harmonics from the physical point of view still remains to be verified. I haven’t done so yet but, meanwhile, I’d like to discuss several other related problems that I discovered while writing and testing this piece of code.

  1. The code that loads and sets up the Moon’s multipole coefficients from “GL0660B.gfc” is borrowed from this post by yzokras:
    Loading normalized spherical harmonics for the moon - #5 by yzokras
    I believe yzokras’s example erroneously attaches the Moon’s harmonics (moonHarmonicsProvider) to the Moon-centered inertial frame. They should rather be attached to the Moon-fixed frame, as I do in my code. If you agree, I can suggest this correction in that thread (although it is 4 yr old).

  2. As you can see, I tried both SingleBodyAbsoluteAttraction and ThirdBodyAttraction for the monopoles. I expected them to be equivalent but, surprisingly, they give very different results (~10^5 km for the final state). I investigated the source code for both classes and I think ThirdBodyAttraction has a flaw. In particular, the acceleration method returns:

Vector3D(gm / (r2Sat * FastMath.sqrt(r2Sat)), satToBody,
-gm / (r2Central * FastMath.sqrt(r2Central)), centralToBody)

This sum of two vectors seems meaningless. These two vectors express the accelerations due to gravity of the third body at the locations of the satellite and the central body (in the latter case, its direction is reversed). I think there is no sense in adding these accelerations and the second term should just be removed. Then, however, this class would be identical to SingleBodyAbsoluteAttraction and the question arises why having two different classes for the same job. So, may be I’m missing something. Please let me know if I’m wrong. The ThirdBodyAttraction class has already caused some confusion before (search “What is the direction of ThirdBodyAttraction acceleration?” in this forum, I can’t post more than 2 links as a new user) but it seems that discussion hasn’t arrived to the essence of the problem.

  1. Finally, the major problem is that my code fails to reproduce distant retrograde orbits. I take the initial conditions from this work (open access):
    Turner, G. (2016). Results of long-duration simulation of distant retrograde orbits. Aerospace, 3(4), 37.
    In particular, I tried several initial conditions from Table 1 on page 20. The position and velocity specified in this table are as discussed on page 11 (LVLH frame):
    position = distance from Moon, along the Moon-to-Earth position vector,
    velocity = Moon-relative velocity in the Earth–Moon orbital plane, perpendicular to the Moon–Earth position vector

I was not able to obtain any Moon-bound trajectories with any of those initial conditions, although they should produce orbits that are stable on a time scale of ~30 years. As a simple check, the code outputs spacecraft-Moon distance with a time step of 1hr. After 60 days (and even just a few days) the SC-Moon distance always goes to ~10^6 km. Any ideas why this doesn’t work are very much appreciated.

Best regards,
Dmitry
NumericalPropagationDRO.java (9.1 KB)