Identifying conjunctions in REBOUND simulations

[nchoksi4456]

Hi all,

I am integrating two-planet systems near first-order resonance:

sim = rebound.Simulation()
sim.integrator = "ias15" 
sim.G = 1
sim.add(m=1)
mu = 3e-5
# Add two planets 1% wide of 3:2 resonance
Delta = 0.01
q=2;
a1 = 1
a2 = (  (Delta+1)*(q+1)/q  )**(2./3)*a1
sim.add(m=mu,a=a1,e=0.0001, inc=0, Omega = 0, omega=0.0, f=0.01 )
sim.add(m=mu,a=a2,e=0.0001, inc=0, Omega = 0, omega=0.0, f=-0.01 ) 
p = sim.particles
sim.move_to_com()

and trying to identify all the times between t = 0 and some given t_maxwhen the longitude of conjunctions = 3 x lambda_prime - 2 x lambda (e.g. equation 3 here) is equal to zero. When this quantity is zero, the longitude of conjunctions points along the reference direction = x-axis.

I've been struggling to come up with a simple and reliable method to identify these events, at some user-specified time precision, for the simulation setup above. I tried modifying the bisecting algorithm used to identify TTVs provided in the REBOUND examples, but I wasn't successful. Any help would be appreciated.

Thanks in advance and apologies if this ends up having a trivial solution!

[Rmelikyan]

I'm not sure if this is the optimal method (i.e. it may be quite slow) but I'd look into using a custom heartbeat function as explained in this example here.

This would allow you to check for conjunction during every time step (you'll probably want to include some reasonable range that is sufficiently close to zero) and then maybe you can save those 'close conjunctions' to a simarchive and finish your more precise analysis from there?

If this becomes prohibitively slow, I'd also encourage you to look into doing this in C since this should be a relatively simple simulation and the heartbeat function is more obviously accessed in the C side of rebound.

I hope that helps! I'm sure someone has a more clever solution than this brute-force check... but that's how I'd start.

Best,
-Robert

[hannorein]

Hi,

I think the Poincare Map example does exactly what you want to do. But let me know if I'm missing something. There might be ways to speed this up, but I'd need to know what accuracy you need.

Hanno

[hannorein]

If you're ok with a time precision of +/- one timestep, then @Rmelikyan's idea should work well. But yes, I'd implement it in C if speed is an issue.

[nchoksi4456]

Thanks Hanno and Robert for your suggestions! The simulations are short, just a few hundred orbits, so I used the heartbeat method In Python with a fixed WHFast timestep of ~1e-3 orbital periods (which gives good enough precision in calculating the time of conjunction).