Оптимизация расчета гравитации для частиц в двумерном пространстве с невесомостью

Я создал небольшую визуализацию частиц в питоне. Я отслеживаю движение частиц в двухмерном пространстве с невесомостью. Поскольку каждая частица притягивает все другие частицы в зависимости от массы частицы и расстояния.

Я сделал визуализацию в pygame, и все работает как план (с расчетом), однако мне нужно максимально оптимизировать расчет. Сегодня система может рассчитывать около 100-150 частиц с нулевой частотой кадров. Я поместил все вычисления в отдельный поток, который дал мне немного больше, но не почти то, что я хочу.

Я смотрел на scipy и numpy, но, поскольку я не ученый или mathguru, я просто запутался. Похоже, я на правильном пути, но я понятия не имею, как это сделать.

Мне нужно вычислить все притяжение всех частиц, которые у меня есть, к петле в петле. И так как мне нужно выяснить, не произошло ли столкновение, я должен сделать то же самое снова.

Мне больно писать такой код ....

У Numpy есть возможность вычислять массив с массивом , однако я не нашел, что бы вычислить все элементы в массиве со всеми элементами из того же / другого массива. Есть один? Если бы это было так, я мог бы создать пару массивов и вычислить намного быстрее, и должна быть функция для получения индекса из двух массивов, где их значения совпадают (Collitiondetect iow)

Вот сегодняшнее вычисление привлекательности / столкновения:

class Particle:
    def __init__(self):
        self.x = random.randint(10,790)
        self.y = random.randint(10,590)
        self.speedx = 0.0
        self.speedy = 0.0
        self.mass = 4

#Attraction    
for p in Particles:
    for p2 in Particles:
        if p != p2:
            xdiff = P.x - P2.x
            ydiff = P.y - P2.y
            dist = math.sqrt((xdiff**2)+(ydiff**2))
            force = 0.125*(p.mass*p2.mass)/(dist**2)
            acceleration = force / p.mass
            xc = xdiff/dist
            yc = ydiff/dist
            P.speedx -= acceleration * xc
            P.speedy -= acceleration * yc
for p in Particles:
    p.x += p.speedx
    p.y += p.speedy

#Collision
for P in Particles:
   for P2 in Particles:
        if p != P2:
            Distance = math.sqrt(  ((p.x-P2.x)**2)  +  ((p.y-P2.y)**2)  )
            if Distance < (p.radius+P2.radius):
                p.speedx = ((p.mass*p.speedx)+(P2.mass*P2.speedx))/(p.mass+P2.mass)
                p.speedy = ((p.mass*p.speedy)+(P2.mass*P2.speedy))/(p.mass+P2.mass)
                p.x = ((p.mass*p.x)+(P2.mass*P2.x))/(p.mass+P2.mass)
                p.y = ((p.mass*p.y)+(P2.mass*P2.y))/(p.mass+P2.mass)
                p.mass += P2.mass
                p.radius = math.sqrt(p.mass)
                Particles.remove(P2)