def doCollision(ncoll, sample1, sample2): for _ in range(ncoll): index1 = np.random.randint(len(sample1)) index2 = np.random.randint(len(sample2)) # By cancelling out 0.5m from the kinetic # energy formula, only v^2 is needed rand1 = sample1[index1] ** 2 rand2 = sample2[index2] ** 2 hdiff = 0.5 * abs(rand1 - rand2) if rand1 > rand2: v1 = sqrt(rand1 - hdiff) v2 = sqrt(rand2 + hdiff) else: v1 = sqrt(rand1 + hdiff) v2 = sqrt(rand2 - hdiff) sample1[index1] = v1 sample2[index2] = v2 return np.concatenate([sample1, sample2])

The function doCollision takes 3 arguments, 
ncoll, sample1 and sample2. 
In the for loop, two random indices are generated, 
one for each sample, and the square of the corresponding 
velocities are calculated. 
Next, the difference of the squares is calculated. 
In the if statement, the kinetic energy is conserved.
If the velocity of the first particle is higher than the second,
the first particle loses energy and the second particle gains energy. 
If the velocity of the first particle is lower than the second,
the first particle gains energy and the second particle loses energy.