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.