Particle order in simulation

[zhexingli]

Hello,

I understand that by default the orbital elements of particles are calculated using Jacobi coordinates, which requires us adding particles to the simulation in order of ascending separation from the host object. I'm curious, what would happen if I don't add particles in ascending separation order? Will rebound still recognize the particle parameters and assumes them to be Jacobian, or is it gonna do something else? (I'm not using the 'primary' flag for the central star.)

Thanks.

[hannorein]

The orbital parameters that you get back will still assume Jacobi coordinates. The fact that they are not order will lead to less physically meaningful orbital parameters, but they are still well defined mathematically. Whether the is a big difference or not depends on the mass ratio of the objects. Note that REBOUND can't possible determine with 100% certainty which particles are closer to the host object. That is why you need to add them in the right order if you want to use Jacobi coordinates.

If this is confusing, do a few test where you add particles in the right/wrong order. To see a large effect choose a mass ratio of 1.

[zhexingli]

Thanks Hanno for the quick answer! Will try to test it out. If I use the 'primary' flag, then there's no need to add them in order, right? How's the use of the two systems different for the purpose of the simulation? Computational speed?

[zhexingli]

Something odd just happened I'm not sure why. When I used this line to initialize a particle: sim.add(m=1.0*9.546256e-4,P=545.76,e=0.01,Omega=1.0,pomega=2.0,inc=np.pi/2) and then use os = sim.calculate_orbits() to get the orbital elements. The values for pomega returned by calling os[0].pomega doesn't match the one I provided in the script. Same thing is happening to my os[0].l call. In this case, when I call os[0].pomega, it gives 0 instead of 2.

Here's the test code I'm using

sim = rebound.Simulation()
sim.integrator = 'whfast'
sim.ri_whfast.safe_mode = 0

sim.units = ["day", "au", "Msun"]

sim.t = 2457812.0
sim.dt = 20

sim.add(m=1.0)

#sim.add(m=0.19*9.546256e-4,P=180.34,e=0.07,\
#        pomega=2.0,l=3.4,inc=np.pi/2)
sim.add(m=1.0*9.546256e-4,P=545.76,e=0.01,\
        Omega=1.0,pomega=2.0,inc=np.pi/2)

ob = sim.calculate_orbits()

sim.move_to_com()

[hannorein]

Because your orbit is in the xy plane, some of the orbital parameters such as the ascending node are ill defined. (wrong explanation)

This is because various angles (l, f, theta, pomega) are defined differently for retrograde orbits. Your orbit has an inclination of exactly pi/2 which puts it right on the discontinuity. Change it to 0.99999pi/2 or 1.00001pi/2 and you should get constant results.

[hannorein]

To answer your other question. Some integrators (IAS15, leapfrog, etc) don't care about the ordering of particles. Other (WHFast, SABA) care. For those who care, you will see a loss in accuracy if they are in an unphysical order. There shouldn't be a significant effect on runtime.

[zhexingli]

Got it. Thank you!