Spikes in semi major axis when increasing output number in a simulation

[VincentPujau]

Environment

• REBOUND Version: 3.28.3
• API interface: Python
• Operating System (including version): Windows 10

Describe the bug Dear all, for my master thesis, I am trying to reproduce close encounters between particles and Earth. I have encounter a problem for which I am not sure what the source is. My personal code is similar to what the example STARMAN is doing therefore I will talk about this code here. I have re ran the example on my computer, firstly the exact same code, and later, I just changed the number of output to get a finer plot. I passed from N=1000 to N=100000. When doing so, I get some kind of huge spikes on the semi major axis either before or after its value change from some constant value to another one (when there is a jump in semi-major axis) like can be seen on the screenshot:

Here is the part of the code I changed: # integrate sim.dt = sim.particles[1].P/60. # small fraction of Mercury's period sim.integrator = "mercurius" N = 10000 times = np.linspace(0.,2.*np.pi*1e5,N) a = np.zeros(N) e = np.zeros(N) for i,t in enumerate(times): sim.integrate(t,exact_finish_time=0) orbit = sim.particles[-1].calculate_orbit(primary=sim.particles[0]) a[i] = orbit.a e[i] = orbit.e

Thanks in advance for your help!

I am not sure why this is happening, and I get a similar problem in my personal code so

[hannorein]

You need to think about what the primary object should be when calculating orbital elements. Right now you use the star (sim.particles[0]). But there are also planets in the system which pull on the star. Calculate the centre of mass of the star and all planets interior to STARMAN. Then use that as the primary when calculating orbital parameters for STARMAN. This procedure is referred to as Jacobi coordinates. It will reduce the spikes but not completely get rid of them because there will still be other interactions with planets that affect the instantaneous orbital elements. I hope this helps.

[hannorein]

I'll close this issue for now but feel free to reopen it if you have further questions.

[VincentPujau]

Hi, thanks for your answer, I have tried to do calculate the center of mass including the star and the planets and use this as reference but I don't get how I can later use orbit = sim.particles[-1].calculate_orbit(primary=com) to have accurate orbital elements if the primary is not the sun. I have tried it and get a = 4000 AU and even negative semi major axis.

[hannorein]

That suggests the centre of mass you calculate is not correct. How do you calculate it?

[VincentPujau]

This is my code

sim = rebound.Simulation()
sim.add(["Sun","Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune"],date="2018-02-10 00:00")
sim.save("ss.bin")

sim.move_to_com()
com = sim.calculate_com()

sim = rebound.Simulation("ss.bin")
sim.add("SpaceX Roadster")

sim.dt = sim.particles[1].P/60. # small fraction of Mercury's period
sim.integrator = "mercurius"  

tmax = 1e3
N = 100000
times = np.linspace(0.,2.*np.pi*tmax,N)

for i,t in enumerate(times):
    sim.integrate(t,exact_finish_time=0)
    orbit = sim.particles[-1].calculate_orbit(primary=com)

I skipped a few lines to make it easier for you

[hannorein]

As I wrote before you need to calculate the centre of mass of all the planets interior to the spacecraft. You can do something like this:

sim = rebound.Simulation()
sim.add(["Sun","Mercury","Venus","Earth"],date="2018-02-10 00:00")
com_inner = sim.calculate_com()
sim.add(["Mars","Jupiter","Saturn","Uranus","Neptune"],date="2018-02-10 00:00")
sim.add("SpaceX Roadster")
sim.save("ss.bin")
a = sim.particles[-1].calculate_orbit(primary=com_inner).a

By the way, these lines don't make sense in your code

sim.save("ss.bin")
sim.move_to_com()
com = sim.calculate_com()
sim = rebound.Simulation("ss.bin")

If you first move the simulation to the centre of mass frame, then the centre of mass will be zero. No need to calculate it! Even more confusingly, you are then loading a simulation which was saved before you moved to the com. I suggest you go through the code line by line and try to understand what it does. Coordinate systems and especially orbital elements can be confusing!