Single Body Attraction vs. Newtonian

I’m baffled by the behavior I’m seeing with single body attraction forces. I’m hoping I’m just missing something but here it goes. I’m exploring this in order to get propagation about the solar system barycenter. Thus I can’t use the standard Newtonian attraction. We have SingleBodyAbsoluteAttraction and SingleBodyRelativeAttraction. I haven’t figured out when to use those or ThirdBodyAttraction in general, but for now none of them work they way I’m expecting to. I’ve created a very simple earth based earth centered example to demonstrate the problem I’m running into in a simplified way. I assumed I could replace the NewtonianAttraction with one of the SingleBody ones but no dice. Code implementing a simple one orbit two body propagation four ways is here:

I’ve run this code two ways, one in Keplerian mode and another in Cartesian mode. The results are equally different in both but it’s easiest to see what’s happening in Keplerian mode. The acceleration computation for the first step of the SingleBody attractor comes up with this:

{-3.4999169942; -1.9056058735; -0.1176217802}

Back of the envelope calculation this is consistent with standard Newtonian gravity GM/r^2. It is then taken into the NumericalPropagator.addNonKeplerianAcceleration function where it is transformed into a yDot of this:
0 = 1.0160761121369433E-12
1 = 2.259799199979748E-16
2 = -2.2445208940156253E-13
3 = 6.314743949604672
4 = -1.3552527156068805E-20
5 = -6.3147439496046704
6 = 0.0

Against an initial state of:
0 = 1.0E7
1 = 1.0E-4
2 = 0.1
3 = 0.3
4 = 0.2
5 = 0.0
6 = 1000.0

For reference, 0 = sma, 1 = ecc, 2 = inc, 3 = pa, 4 = raan, 5 = ta, 6 = mass

What we should see is that it’s only adjusting true anomaly (or mostly adjusting that anyway). This instead is cranking around the argument of periapsis by 344 degrees and true anomaly also about 1.8 degrees. Not only should the change be just in true anomaly but it should be smaller, about 0.04 degrees. When we look at Newtonian Attraction (which goes through a very specific function for adding Newtonian Acceleration) we get this for yDot:
0 = 0.0
1 = 0.0
2 = 0.0
3 = 0.0
4 = 0.0
5 = 6.314743949604675E-4
6 = 0.0

The second step evaluation of yDot is even worse:
0 = 1.1136507479936313
1 = 5.56825368424313E-4
2 = 8.453974655465556E-13
3 = 2.9769301951916534
4 = -6.776263578034403E-21
5 = -2.9769301951935363
6 = 0.0

I’m continuing to explore but does anything jump out as to why this wouldn’t just work? I know why it’s not going through some specific Newtonian force shortcut but either way it seems the net effect on the state should be the same when going through a general perturbation code path, minus some small deviations for round off error and the like. Again the link to the source code is above.

I believe I’ve figured out why and how to fix this…updates pending…

The problem I found was in the NumericalPropagator code. When there is no Newtonian force it never adds the velocity term to the yDot for position. Therefore when we solely have a Single Body Attractor in the model we never get velocity contributions to state change, we only get the acceleration from the point mass. When in Cartesian mode we have enough information to perform this calculation but we don’t for other state spaces. I’ve committed an interim fix for this with this commit.. I’ve illustrated it working well in the two-body cases for both Earth and Sun in this repo. With this change I think we are one more step closer to being able to do propagations around a barycenter by having the primary gravitational force being modeled with the Single Body Attractor.

I’ve run a test comparing a force model where both have all the planets and moon but the main body is modeled as NewtonianAttraction or a SingleBodyAbsoluteAttraction object for a sun-centered heliocentric propagation. Over 50 years they match to within less than 1.5 meters and 3e-7 m/s although it is highly variable over that 50 years. The below graph is the position, but the velocity chart matches it nearly identically in terms of shape:

Hi @hankgrabowski

Can you open an issue for this bug? I think your fix is the correct one. Therefore, as you have writing access to the repository, you can directly add the changes. As Orekit developer, a cleaner way is to open a merge-request and immediatelly click on merge when pipeline succeeds.

Do you have a unit test for this bug?

I don’t know if this difference is significant or not. In my opinion, a maximum of 1.5 metres, over 50 years of orbit propagation, is a very good order of magnitude.

Regards
Bryan

I don’t even direct commit to development branches on my own repositories much less collaborative ones :). I would normally want an independent review on a collaborative project. Could I assign you as the reviewer? I still need to put a unit test in as well. I just wanted to have someone look at my WIP code to make sure it was in the realm of the right thing to do.

Yeah I’m of the same mind on all of those points. It actually oscillates over the 50 years so it’s more than just aggregation of round off error. I’d like to study the nature of it and I may do that as an artifact of continuing to try to get my heliocentric work going but the fix is still better than what is going on right now.

I understand :slightly_smiling_face: That’s why merge request is an interesting tool to able other contributors to look at the code before merging it.

Yes you can. I’d be happy to help you.