Spikes in semi major axis when increasing output number in a simulation about rebound HOT 6 CLOSED

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.

from rebound.

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

from rebound.

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.

from rebound.

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

from rebound.

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

from rebound.

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!

from rebound.