Additional forces double-counted in integrator_bs.c

[alekseygenerozov]

Dear Hanno,

I believe I found an issue in the Burlisch-Stohr integrator in rebound. In particular, I think it is double-counting the accelerations from additional forces.

I am sharing some code that illustrates this (https://gist.github.com/alekseygenerozov/fb28ca567a495100d0fce9b24170e8f4). This is a modified example problem from reboundx (in examples/central_force/problem.c). Basically the setup is such that (i) the additional force dominates (to better illustrate the problem) and (ii) the integration time is very short so you see the initial acceleration of the particle (the planet in this example). I am doing this test with the latest versions of rebound and reboundx

The acceleration in this example is almost exactly twice greater than expected. Here is the output I get

t ax ay az: 1.000000e-15 3.125000e-04 4.784160e-25 0.000000e+00 x y z: 0.800000 0.000000 0.000000

The problem goes away when I change

sim->integrator=REB_INTEGRATOR_BS;

To sim->integrator=REB_INTEGRATOR_IAS15;

I believe the issue is line 327 in src/integrator_bs.c. This line seems to be superfluous since additional_forces is also called in reb_update_acceleration. When I comment out this line the issue goes away.

Best,

Aleksey

[hannorein]

Hi Aleksey! Thanks for the detail report! I think you're correct! But let me double check ...

[hannorein]

Thanks again for tracking down the issue. It should now be fixed in the latest version. I also added a unit test for this case.