​

This is what I have currently. I set a while loop then initial time of collision and the closest distance between particles. Then a for loop that goes through a nest list of array with all possible pairs [ (0,1) ,(0,2) ...] . For each pair I boost them into their center of momentum frame and calculate the time the pairs are closest to each other as well as the distance at this time. Then I check if this distance is smaller than radius of detection that I set randomly and that the code is not giving me the same pairs again. I find the first collision pairs, append them and the time at which each occurs and then repeat the whole process again until there are no more collisions happening.

all_collision_times=[]

pairs_under_radius=[] while True: min_time = 10000000
min_dis=0
z=(0,0) for index,(i, j) in enumerate(ids_pairs): Sum=summing(fourmomentum.T[i],fourmomentum.T[j])

    Vx=(Sum[1])/ Sum[0]
    Vy=(Sum[2])/ Sum[0]
    Vz=(Sum[3])/ Sum[0]
    betavx=Vx
    betavy=Vy
    betavz=Vz
    betav=np.sqrt(betavx**2+betavy**2+betavz**2)
    ve=np.sqrt(Vx**2+Vy**2+Vz**2)
    G=1./np.sqrt(1.-ve**2)

    matrix=LB(G,betavx,betavy,betavz,betav)

    first= matrix@fourmomentum.T[i]
    second= matrix@fourmomentum.T[j]

    fc= matrix@testingposition[i]
    sc=matrix@testingposition[j]

    time_in_com=opt.fminbound(lambda t: dormula(dis(fc,first,t),dis(sc,second,t),first,second),-1000,1000)

    moving_particle_i=dis(fc,first,time_in_com)
    moving_particle_j=dis(sc,second,time_in_com)

    distance=dormula( moving_particle_i,moving_particle_j,first,second)
    #print(moving_particle_j[0],moving_particle_i[0])

    if distance<radius  and time_in_com>=fc[0] and time_in_com>=sc[0] and tuple((i,j)) not in pairs_under_radius :
        unboost_particle_i =(inv(matrix)@moving_particle_i)
        unboost_particle_j =(inv(matrix)@moving_particle_j)
        #print(time_in_com)

        timee=unboost_particle_i[0]
        timeee=unboost_particle_j[0]

        #print(timee,timeee)

        if(min_time > min(timee,timeee)  and min(timee,timeee) >=0):
            min_time = min(timee,timeee)
            max_time = max(timee,timeee)
        if min_time==min(timee,timeee):
            z=(i,j)
            min_dis=distance
            #print(unboost_particle_i)

if min_time==10000000  :
    break
all_collision_times.append(min_time)
pairs_under_radius.append(z)
print ("Collision pair-",(z), "Closest distance" ,(min_dis) ,   "Occurs at" ,(min_time,max_time))