<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 User Defined Rebound Collision Resolutions
In the CloseEncounter example, we discuss methods for resolving collisions in REBOUND through exceptions and the use of the sim.collision_resolve = "merge" method.
Using the same 3-Body setup, let us explore how to define and implement the same collision resolution function in python and pass it to the sim.collision_resolve function pointer.
Step1:  To reiterate the previous method, let's run the built-in merge collision resolution method
Step2:  We can see above that two particles merged into one with a combined mass of 0.006.
Let's now try to implement this collision function ourselves!
To do this, we need to write a function which we can pass to sim.collision_resolve. In this case let's define my_merge. 
Now, whenever a collision occurs, REBOUND will pass our function two parameters
Step3:  Now we can set our new collision resolution function in the simulation object.
<ASSISTANT_TASK:>
Python Code:
import rebound
import numpy as np
import matplotlib.pyplot as plt
def setupSimulation():
    ''' Setup the 3-Body scenario'''
    sim = rebound.Simulation()
    sim.integrator = "ias15" # IAS15 is the default integrator, so we don't need this line
    sim.add(m=1.)
    sim.add(m=1e-3, a=1., r=np.sqrt(1e-3/3.)) # we now set collision radii!
    sim.add(m=5e-3, a=1.25, r=1.25*np.sqrt(5e-3/3.))
    sim.move_to_com()
    return sim
Explanation: User Defined Rebound Collision Resolutions
In the CloseEncounter example, we discuss methods for resolving collisions in REBOUND through exceptions and the use of the sim.collision_resolve = "merge" method.
Using the same 3-Body setup, let us explore how to define and implement the same collision resolution function in python and pass it to the sim.collision_resolve function pointer.
End of explanation
sim = setupSimulation()
sim.collision = "direct"
sim.collision_resolve = "merge" # Built in function
print("Particles in the simulation at t=%6.1f: %d"%(sim.t,sim.N))
print("System Mass: {}".format([p.m for p in sim.particles]))
sim.integrate(100.)
print("Particles in the simulation at t=%6.1f: %d"%(sim.t,sim.N))
print("System Mass: {}".format([p.m for p in sim.particles]))
Explanation: To reiterate the previous method, let's run the built-in merge collision resolution method
End of explanation
def my_merge(sim_pointer, collided_particles_index):
    sim = sim_pointer.contents # retreive the standard simulation object
    ps = sim.particles # easy access to list of particles
    i = collided_particles_index.p1   # Note that p1 < p2 is not guaranteed.    
    j = collided_particles_index.p2 
    # This part is exciting! We can execute additional code during collisions now!
    fig, ax = rebound.OrbitPlot(sim, xlim = (-1.3, 1.3), ylim = (-1.3, 1.3), color=['blue', 'green'])
    ax.set_title("Merging particle {} into {}".format(j, i))
    ax.text(ps[1].x, ps[1].y, "1"); 
    ax.text(ps[2].x, ps[2].y, "2")
    # So we plot the scenario exactly at the timestep that the collision function is triggered
    # Merging Logic 
    total_mass = ps[i].m + ps[j].m
    merged_planet = (ps[i] * ps[i].m + ps[j] * ps[j].m)/total_mass # conservation of momentum
    # merged radius assuming a uniform density
    merged_radius = (ps[i].r**3 + ps[j].r**3)**(1/3)
    ps[i] = merged_planet   # update p1's state vector (mass and radius will need corrections)
    ps[i].m = total_mass    # update to total mass
    ps[i].r = merged_radius # update to joined radius
    return 2 # remove particle with index j
Explanation: We can see above that two particles merged into one with a combined mass of 0.006.
Let's now try to implement this collision function ourselves!
To do this, we need to write a function which we can pass to sim.collision_resolve. In this case let's define my_merge. 
Now, whenever a collision occurs, REBOUND will pass our function two parameters:
sim_pointer: a pointer to the simulation object which the collision occurred in.
Because it is a ctypes pointer, you will need to use the .contents attribute to access the simulation object
collision: this structure contains the attributes .p1 and .p2 which are the indices of the two particles involved in the collision
Using these inputs, we can define the necessary logic to handle the collision. The return value of our function determines how REBOUND proceeds afterwards:
0: Simulation continues without changes
1: remove p1 from simulation
2: remove p2 from simulation
Let us look at how this information can be used to implement the logic of the merge method for colliding particles in a totally inelastic collision.
End of explanation
sim = setupSimulation()
sim.collision = "direct"
ps = sim.particles
sim.collision_resolve = my_merge # user defined collision resolution function
sim.integrate(100.)
Explanation: Now we can set our new collision resolution function in the simulation object.
End of explanation
<END_TASK>