Fixed circular orbit and position

[AstroEloy]

Hi! I would like to know if any option is already implemented in rebound to define perfectly circular orbits and not have to integrate the equations of motion to propagate in time. I could overwrite in each iteration the value of the position and velocity by the value of the circular motion, but I think it could be very beneficial to have this option native both to save computational cost and to improve simulation control in certain cases.

I would also like to ask if the main body of the system can be fixed to always be at position [0,0,0].

Thank you very much!

Eloy

[hannorein]

There are several ways you could do that:

• Set, set Nactive=1, use WHFast.
• At every step simply remove the particle and add a new particle with the mean longitude changed.

I'd go with the first option because that makes it trivial to add gravitational interactions back in later.

[AstroEloy]

so it is not possible to do this easily with IAS15?

Thanks!

[hannorein]

I'm not entirely sure what you want to do. You can use IAS15 as well, but that probably doesn't make a lot of sense. The two body problem is integrable, so you don't need an integrator.

[AstroEloy]

My plan is to do a simulation with thousands of particles, but fixing the Sun at position [0,0,0] and fixing the orbits of the planets as perfectly circular (and not integrating them, only the particles).

[hannorein]

Ah. That's a very different problem. Maybe look into test particles?

https://rebound.readthedocs.io/en/latest/ipython_examples/Testparticles/

[AstroEloy]

I may be understanding something wrong... I'm trying to explain again:

I want to position the Sun at the center of my system, and for example, a planet with mass rotating in a perfect circular orbit (forever). I don't want these two bodies to interact with each other, but I do want thousands of particles that I am going to add to feel the gravitational pull of the planet and the sun.

[hannorein]

I can only repeat what I said above. I think you might want to look into test particles. ;-)

[AstroEloy]

I can't figure out how to tell the particles not to interact with each other but with the planet and the sun (test_particles=0), and at the same time tell the planet and the Sun not to interact with each other but affect the particles (being the two N_actives).

I think this is not detailed in the link you have shared with me.

[hannorein]

The example does exactly what you've been asking about. If you've modified it and it doesn't work anymore, feel free to post your code here and we can have a look.

[AstroEloy]

The position of the first active body is not maintained at [0,0,0]. I guess because is affected by the second active object. This is what I am trying to avoid I would like both active body remain in their fixed position (one in [0,0,0] and the other one in a circle orbit)

import rebound

sim = rebound.Simulation()
sim.add(m=1.)
sim.add(m=1e-1, a=1, e=0.05)
sim.move_to_com()
sim.integrator = "whfast"
sim.dt = 0.05

import numpy as np
N_testparticle = 100
a_initial = np.linspace(1.1, 3, N_testparticle)
for a in a_initial:
    sim.add(a=a,f=np.random.rand()*2.*np.pi) # mass is set to 0 by default, random true anomaly

sim.N_active = 2

t_max = 200.*2.*np.pi
N_out = 10
xy = np.zeros((N_out, N_testparticle, 2))
times = np.linspace(0, t_max, N_out)
for i, time in enumerate(times):
    sim.integrate(time)
    print(sim.particles[0].x, sim.particles[0].y)

[hannorein]

REBOUND can be used in any inertial frame. In the example above you work in the inertial centre of mass frame. What you want is a non-inertial frame. The simplest way to do this is to simple move to the heliocentric frame every time you want an output: sim.move_to_hel().

[AstroEloy]

Thank you Hanno. The frame change is what was confusing me.

[hannorein]

Got it! Some other integrator packages work in a non-inertial frame, so this can be confusing. But I think an inertial frame is ultimately much more physical and easier to think about...

[dmhernan]

Hi @AstroEloy @hannorein , I think the WHFast option with Nactive=1 works, but isn't that computing the Kepler advancer, which requires several trig functions plus solve an implicit equation? Maybe what @AstroEloy is looking for is just computing positions via x = cos(w t), y = sin(w t). And adding in a rotation for inclinations. Could this be a simple "circular motion" integrator? Not sure in practice if it saves time.

[AstroEloy]

Thank you @dmhernan . This is in fact what I was looking for: an option to define an orbit via circular motion equation (for avoiding integration). I think this could save some time for very large population propagate long time periods.

[danielk333]

I would say: code up the study using currently available options and THEN if it proves to be too slow, look into optimization.

My gut feeling is telling me that thousands of particles will dominate the computation time whatever you do with the massive planets (my coffee-deprived brain might be off here but step for massive planets should be approx N_massive*(N_massive-1)/2 force calculations while the particles have N_test*N_massive, so at N_test=1000 and N_massive=9 its already 9000 versus 36)

[hannorein]

I agree with @danielk333!

[dmhernan]

Oh, sure, although I guess if each test particle simulation is done in parallel, separately (i.e., 1000 simulations), then it makes an impact. Another consideration is roundoff error. Using Nactive=1 will accumulate roundoff error, while the circular motion integrator avoids roundoff error. One would want to redefine time so it's not increasing without bound.

[danielk333]

Yes parallelization does change the optimization options. But one should keep in mind that 1000 processes in parallel would probably have as much communication overhead time as it would have force computation time unless there is a large amount of massive bodies or each thread uses a cache for massive body states to avoid node communication.

If I ever get some extra time, the propagation of test particles using cached integrator steps is one addition I would like to contribute to rebound (very useful for statistical sampling on HPC's).

The roundoff error is more tricky, just overriding the state with the circular motion at each step will not totally do it due to the WHFast intermediate steps. But then again... the important question is probably if it matter for the time scales in question?

[dmhernan]

Sounds good, although wouldn't the 1000 processes be completely independent without any communication? I could be misunderstanding. As for WHFast, wouldn't the goal be to avoid using it entirely? Anyway, doesn't matter, depends on what @AstroEloy wants to do.

[hannorein]

I think this discussion is a bit too abstract at the moment. We'd need to know exactly what the setup looks like. If someone wants to simulate a bunch of test particles under the influence of one massive body on a circular orbit, then WHFast with N_active = 2 is the way to go. This requires no code changes, and is pretty much ideal in terms of runtime.

I'll close this issue for now. Thanks to everyone who chimed in! @AstroEloy, feel free to open a new issue if something else comes up. For example, if you were to post your full setup, and let us know your precise requirements in terms of precision, number of particles, etc., then we could brainstorm on how to optimize it further.