Collisions Won't Turn Off

[jflb96]

At least, that's what I assume is going on. Whenever I put two objects too close together, the code gets stuck processing the first step until eventually throwing out the following error:

KeyboardInterrupt Traceback (most recent call last) Cell In[19], line 43 41 print(o) 42 for i in np.arange(-5,200): ---> 43 sim_sol_inner.integrate(i6np.pi/200) 44 get_ipython().run_line_magic('matplotlib', 'inline') 45 op = rb.OrbitPlot(sim_sol_inner, unitlabel="[AU]", color=True)

File ~/.local/lib/python3.10/site-packages/rebound/simulation.py:2114, in Simulation.integrate(self, tmax, exact_finish_time) 2112 pass # User caused exit. Do not raise error message 2113 if ret_value == 6: -> 2114 raise KeyboardInterrupt 2115 if ret_value == 7: 2116 raise Collision("Two particles collided (d < r1+r2)")

KeyboardInterrupt:

However, I have the collisions turned off at the start, as far as I can tell:

sim_sol_inner = rb.Simulation() sim_sol_inner.collision = "none"

I'm just not sure what's gone wrong, or how to make it not-gone-wrong. Any help would be excellent.

[hannorein]

If you read the error message carefully, you'll see that it has nothing to do with collisions. It looks like you interrupted the integration manually. I'm not able to help you out without more information. Show me your code, tell me what you want to do, etc.

PS: I'm much more inclined to respond to this kind go question if you introduce yourself.

[jflb96]

Hi! I'm John! I'm a physics student trying to get a cycler orbit set up between Earth and Ceres, but whenever I put something too close to Earth the code freezes until I manually interrupt it. Maybe it is a coincidence that the next line down is talking about collisions and the problem only comes up when two things are right on top of each other, but that's the only thing that I can think of that would be the source of the problem.

import rebound as rb
import math as maths
import glob
from PIL import Image
import numpy as np
import os
rad_e = 6371
rad_s = 696340
mass_e = 5.972e24
mass_s = 1.989e30
au = 1.496e8
sec_un = 5022440.76
ecc_e = 0#.0167
smjax_e = 1
theta = np.pi/9
quad_a = float(input('Give a semi-major axis for the cycler'))
quad_b = np.cos(theta)/(1+ecc_e*np.cos(theta))
quad_c = (smjax_e*(1-ecc_e*ecc_e)-quad_a)/(1+ecc_e*np.cos(theta))
quad_x = (-quad_b + np.sqrt(quad_b*quad_b-(4*quad_a*quad_c)))/(2*quad_a)
sim_sol_inner = rb.Simulation()
sim_sol_inner.collision = "none"
sim_sol_inner.add(m=1.)                # Central object
sim_sol_inner.add(m=0.16601e-6, a=.39, e=0#.2056, #f=-26*np.pi/180
                 ) # Mercury
sim_sol_inner.add(m=2.4478383e-6, a=.72, e=0#.0068, #f=-26*np.pi/180
                 ) # Venus
sim_sol_inner.add(m=3.04043263333e-6, a=1., e=0#.0167, 
                  ,f=np.pi/9, r=4.2635e-5
                 ) # Earth
sim_sol_inner.add(m=0.3227151e-6, a=1.52, e=0#.0934, f=23.5*np.pi/180
                 ) # Mars
sim_sol_inner.add(m=5.21e-10, a=2.77, e=0#.0785, f=-26*np.pi/180
                 ) # Ceres
#sim_sol_inner.add(m=954.79194e-6, a=5.2, e=0.0484, f=-26*np.pi/180) # Jupiter
sim_sol_inner.add(m=9.40123674e-26, a=quad_a, e=quad_x, f=np.pi/9, r=6.6846e-10
                  #x=1., vy=np.cos(np.pi*25.7/180)*(34.9*sec_un/au), vx=np.sin(np.pi*25.7/180)*(34.9*sec_un/au)
                 ) # new cycler
for p in sim_sol_inner.particles:
    print(p.x, p.y, p.z)
for o in sim_sol_inner.calculate_orbits(): 
    print(o)
for i in np.arange(0,300):
    sim_sol_inner.integrate(i*6*np.pi/100)
    %matplotlib inline
    op = rb.OrbitPlot(sim_sol_inner, unitlabel="[AU]", color=True)
    op.ax.set_title('Frame'+str(i+1))
    #op.particles.set_sizes(r_store)
    #op.primary.set_sizes([100*total_strad/strad_count])
    op.primary.set_color("yellow")
    op.primary.set_edgecolor("black")
    op.fig.savefig('ceres_cyc_e0/cyc_solsys_ceres_frame_'+str(i+1)+'.png', bbox_inches="tight")
    print('cyc_solsys_ceres_frame_'+str(i+1)+'.png')
    print(round(100*(i+1)/300,3),'%')
def make_gif(frame_folder):
    frames = [Image.open(image) for image in sorted(sorted(glob.glob(f"{frame_folder}/*.png"), key = lambda 
                                                           x: (os.path.splitext(os.path.basename(x))[0])), key = len)]
    frame_one = frames[0]
    frame_one.save("cyc_ceres_gif.gif", format="GIF", append_images=frames,
               save_all=True, duration=50, loop=0)
    frame_one.close()
    print('All Done! :)')
make_gif('ceres_cyc')

The above is my code. If you change the f values in Earth and the new cycler to not be pi/9, the code runs. If you don't, it freezes.

[hannorein]

Hi John, I suspect your timestep is just getting very small. Try setting a reasonable value for min_dt (see e.g. https://rebound.readthedocs.io/en/latest/integrators/#ias15)

[jflb96]

Would that not be contraindicated by the code working with the same timestep but different starting positions?

[hannorein]

No. How small the timestep gets can be somewhat random (it depends on how close the encounter is, etc).

[jflb96]

I’d assumed that the timestep would be constant unless it was told not to be. I’ll give that a look.

Thanks for your help!

ETA: Yep, that was it! :-)