Optimizing gravitation calculation for particles in a zero gravity 2d space

Question

I've have created a small visualisation of particles in python. I'm caclulation the movement of particels in a 2D space with zero gravity. As each particle attracts all other particles based on the particle mass and distance.

I made a visualsation in pygame and everything works as plan (with the caluclation), however i need to optimize the calculation extreamly. Today the system can calculate about 100-150 particles in a deacent framerate. I put all the calculation in a seperate thread that gave me some more but not nearly what i want.

I've looked at scipy and numpy but since I'm no scientist or mathguru i just get confused. It looks like I'm on the right track but i have no clue howto.

I need to calculate all the attraction on all particles i have to a loop in a loop. And since I need to find if any have collided, i have to do the same all over again.

It breaks my heart to write that kind of code....

Numpy has the ability to calculate array with array, however i haven't found any what to calculate all items in array with all the items from same/another array. Is there one? If so i could create and couple of arrays and calculate much faster and there must be a function to get index from to 2 arrays where their values match (Collitiondetect iow)

Here is todays attraction/collsion calculation:

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)

Answer 1

You can first try to work with complex numbers: the relevant gravitation and dynamics formulas are very simple in this formalism, and can also be quite fast (because NumPy can do the calculation internally, instead of you handling separately x and y coordinates). For instance, the force between two particules at z and z' is simply:

(z-z')/abs(z-z')**3

NumPy can calculate such a quantity very quickly, for all z/z' pairs. For instance, the matrix of all z-z' values is simply obtained from the 1D array Z of coordinates as Z-Z[:, numpy.newaxis] (the diagonal terms [z=z'] do require some special care, when calculating 1/abs(z-z')**3: they should be set to zero).

As for the time evolution, you can certainly use SciPy's fast differential equation routines: they are much more precise than the step by step Euler integration.

In any case, delving into NumPy would be very useful, especially if you plan to do scientific calculations, as NumPy is very fast.