Hi, could you help me please with a numerical propagation problem. I’m trying to simulate Earth-to-Moon transfers and also distant retrograde orbits and I’d like to take into account the nonsphericity of the fields of both Earth and the Moon throughout the trajectory. The SingleBodyAbsoluteAttraction class seemed perfect for this at first but sadly it takes into account only the point-like attraction. I also considered using HolmesFeatherstoneAttractionModel and adding ThirdBodyAttraction as a perturbation, but the latter also takes into account only the point-like attraction. CR3BPForceModel is point-like as well. Is it even possible to do this with the built-in functionality or I need to create my own ForceModel?
I am not sure if it is possible or not!
What you may try is add two different Holmes-Featherstone models to your numerical propagator. Internally, the acceleration method in this class starts by converting from whatever frame the spacecraft state is into the attracting body frame, then compute a vector in this frame, and then convert it back to the inertial frame that is used in the spacecraft state. So nothing would prevent to have several models active at any given time, each one knowing its body frame and performing its conversions.
Let us know if it works, and let us know also if it doesn’t so we can fix it.
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.
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).
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:
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.
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.
Regarding the propagation with the SingleBodyAbsoluteAttraction and the difference with the ThirdBodyAcceleration, I have opened a dedicated post recently (“Propagation without central body”) because I faced some problems with it. You can have a look at it. Basically, what I suspect is that the SingleBodyAbsoluteAttraction only works when integrating in the Solar System Barycenter reference frame. If you don’t propagate in this frame, you should add inertial accelerations which I supsect to have a problem. So my first advice for you would be to try to propagate your orbit in the Solar System Barycenter reference frame (CelestialBodyFactory.getSolarSystemBarycenter().getInertiallyOrientedFrame()). Let me know if this helps…
The ThirdBodyAttraction is not erroneous. As its name suggests, it’s only meant to be used when there is a central body (via NewtonianAttraction), in other words for perturbed two body problem. Basically, once you add other bodies, you need to take into account their pull on the reference frame origin, otherwise it is not inertial anymore. Check out textbooks like Vallado’s or Montenbruck&Gill’s you’ll see this term you’re talking about.
Anyway since you’re doing three body problem stuff, you don’t want to use this class. I’m not very familiar with it, but Orekit has some bricks in this area. Maybe some users have more insight.
Thank you for your advices. They helped a lot. Indeed, for this problem one does not need ThirdBodyAttraction (I checked Vallado’s book and now I can see this class does not have any issues). Also, my old code had a problem of propagating in an Earth-centered inertial (ECI) frame but not taking inertial forces into account. I also missed the fact that in Turner’s paper the initial states are specified in the LVLH frame and didn’t transform from LVLH into ECI properly. In the updated version (attached) I switched to propagating in the ICRF frame, which does not need inertial forces, and fixed the transformation from LVLH to GCRF. Now the initial data from Turner’s table 1 successfully produce distant retrograde orbits.
To additionally check that I’m doing things right I considered the two spacecraft that have so far been injected into DROs, i.e. Orion of the Artemis 1 mission and Chang’e 5. For the latter, unfortunately, it turned out there is no ephemeris available (at least at Horizons the only available data are for the booster and not the orbiter). With Orion I had some success, although the spacecraft spent only a few days in a DRO. I picked an initial DRO state at 2022/11/26 00:03:00 TDB from Horizons and propagated it into a state 4 days later (still a DRO). The discrepancy between the final state produced by Orekit and the actual one from Horizons is of order 10 km for each component of the position vector. The relatively poor agreement can be easily understood by looking at the residuals:
The spacecraft obviously experienced several velocity corrections in the considered time interval (most likely due to reaction wheel desaturations) and these non-gravitational forces are of course not taken into account by Orekit. Just before the first velocity correction at Nov-26 21:03 TDB, 17 hr into the propagation interval, the residuals are 200m for the X and Z position components, and 500m for the Y component. This is still quite large but I believe this could be improved by a better modelling of the SRP force.
