Garrett

test.py

sim = rebound.Simulation()
sim.integrator = "ias15"
sim.add(m=1.)
sim.add(m=1e-3, a=1., r=np.sqrt(1e-3/3.))
sim.add(m=5e-3, a=1.25, r=1.25*np.sqrt(5e-3/3.))
sim.collision="direct"
sim.collision_resolve =  "halt"
E0 = sim.energy()
A0 = sim.angular_momentum()[2]
try:
    sim.integrate(40)
except rebound.Collision:
    E_before = sim.energy()
    A_before = sim.angular_momentum()[2]
    sim.collision_resolve = "merge"
    sim.step()
    E_collision = E_before-sim.energy()
    A_collision = A_before-sim.angular_momentum()[2]
sim.integrate(40)

print("Rel. energy error:  %e"%((E0-sim.energy()-E_collision)/E0))
print("Rel. angular momentum error:  %e"%((A0-sim.angular_momentum()[2]-A_collision)/A0))