Universal attraction: gravity

17 Feb 2011 Code on Github

Introduction

Now we have a physics engine/module we can create a number of different physical simulations. One example would be a simple model of a gas cloud collasping to form a star or solar system. In this tutorial we will:

Model gravitational attraction between all particles
Allow two particles to coalesce into a single larger particle

The simulation we're going to build in this tutorial is of a cloud of dust particles. Every particle will attract every other particle due to gravity, which will cause the cloud to collapse and form larger accumulations of particles. These could represent stars and planets and with a bit of luck should have smaller particles orbiting them.

I have since made an online 3D version.

Setting up the simulation

We'll start with a basic Pygame program, similar to the one we made in tutorial 1, or how we started our test program in the last tutorial. Whenever we write a Pygame program we're going to start with something like this.

import random
import pygame
import PyParticles

(width, height) = (400, 400)
screen = pygame.display.set_mode((width, height))
pygame.display.set_caption('Star formation')

running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False

    pygame.display.flip()

Now we create an Environment object before the main loop:

universe = PyParticles.Environment((width, height))
universe.colour = (0,0,0)
universe.addFunctions(['move'])

We have created our universe. It is black, (0,0,0), and only calls the particles' move() function when updating. Particles don't experience drag (there is no air in space) or bounce when the reach the edge of the screen (there are no walls in space).

Now we add some particles to represent dust particles:

def calculateRadius(mass):
    return 0.4 * mass ** (0.5)

for p in range(100):
    particle_mass = random.randint(1,4)
    particle_size = calculateRadius(particle_mass)
    universe.addParticles(mass=particle_mass, size=particle_size, colour=(255,255,255))

This loop adds to our universe 100 white particles with masses between 1 and 4. It also defines the size of the particle, which proportional to the square root of its mass (because this a 2D simulation) as defined by the calculateRadius() function. In this simulation, the mass of the particles is the most important parameter; the size is only really important for the display. I decided to multiply the square root of the mass by 0.4 purely so all the initial particles have a radius less than 1 unit.

Now we can display our particles by adding the following code before we flip the display:

universe.update()
screen.fill(universe.colour)

for p in universe.particles:
    p.move()
    if p.size < 2:
        pygame.draw.rect(screen, p.colour, (int(p.x), int(p.y), 2, 2))
    else:
        pygame.draw.circle(screen, p.colour, (int(p.x), int(p.y)), int(p.size), 0)

This loops through all the particles, drawing a 2 by 2 rectangle, if it's small, or a circle. The reason I did this is that circles with a radius of 1 look odd and I want to be able to see all the particles, even if they are very small. You'll notice that we're also calling the universe's update() function, so all the particles fly off the screen. That is pretty much it for this program, see how easy it was. Now we need to add a couple of extra functions to our PyParticles module.

Gravity

The driving force behind the particles' movement in this simulation will be gravity. Unlike in previous simulations where gravity was a constant downward force (towards the floor - effectively a planet-size, immovable particle), in this simulation gravity will work between particles. For this we need to give the Particle class an attract() function. We can keep adding functions to the PyParticles module as we need them so long as they don't interfere with any of the other functions. Somewhere in the Particle class, add:

def attract(self, other):
    dx = (self.x - other.x)
    dy = (self.y - other.y)
    dist = math.hypot(dx, dy)
        
    theta = math.atan2(dy, dx)
    force = 0.2 * self.mass * other.mass / dist ** 2
    self.accelerate((theta - 0.5 * math.pi, force/self.mass))
    other.accelerate((theta + 0.5 * math.pi, force/other.mass))

This function takes another particle as a parameter and calculates the distance and angle between them. It then calculates the gravitational force between them, as defined by Newton. Both particles then experience and acceleration in the direction of attraction based on that force and their mass. I picked a gravitational constant of 0.2 because I thought it gave a nice result. You can try other values and see how it effects the simulation.

Since the attract() function requires two particles, we add it to the Environment's function dictionary like so:

'collide': (2, lambda p1, p2: collide(p1, p2)),
'attract': (2, lambda p1, p2: p1.attract(p2))}

I added it under the collide() function so you can see the similarity. The only real difference is that collide() takes two particles as parameters whilst attract() is a particle method that takes a second particle as a parameter. It would probably make sense to change collide() into a particle method, but I left it this way to show how the different approaches are both easily described by a lambda function.

We can now return to our simulation and add this functionality:

universe.addFunctions(['move', 'attract'])

Collisions

If you run the simulation now you'll see a random spread of white particles slowly move towards the centre, occasionally spinning around one another. Eventually they move past each other and the cloud starts expanding again. What we want is for particles that get very close to coalesce into a single, larger particle.

We've previously tested for particle collisions, however the current function causes particles to bounce, so we need to write a function to cause particles to combine:

def combine(p1, p2):
  if math.hypot(p1.x - p2.x, p1.y - p2.y) < p1.size+p2.size:
    total_mass = p1.mass + p2.mass
    p1.x = (p1.x*p1.mass + p2.x*p2.mass)/total_mass
    p1.y = (p1.y*p1.mass + p2.y*p2.mass)/total_mass
    (p1.angle, p1.speed) = addVectors((p1.angle, p1.speed * p1.mass / total_mass), (p2.angle, p2.speed * p2.mass / total_mass))
    p1.speed *= (p1.elasticity * p2.elasticity)
    p1.mass += p2.mass
    p1.collide_with = p2

This function tests whether two particles are overlapping as before, but now if they are, it combines them by converting the first particle, p1, into a larger particle. This larger particle has a position and velocity based on a weighted average of the two original particles position and velocity. The average is weighted by the particles' masses, so that more massive particles have more of an effect on the final values. The mass of this particle is the sum of the original particles' masses.

I also decided to reduce the speed of the combined particle using elasticity, which describes energy lost as heat during the collision, but isn't necessary. Finally, we give this large particle a collide_with attribute and set it to equal to p2, which hasn't been changed, but which we need to remove.

Removing particles

Back in the simulation code, we can loop through the particles and test, which have combined. We test particles by seeing whether they have a collide_with attribute. If they do then we add the particle it collided with to a list of particles to remove. We also re-calculate the radius of the first particle to take into account its new mass. Finally, we remove the collide_with attribute since it we have dealt with that collision. Since we have to loop through the particles to display them, we can add this code into that loop.

particles_to_remove = []
for p in universe.particles:
    if 'collide_with' in p.__dict__:
        particles_to_remove.append(p.collide_with)
        p.size = calculateRadius(p.mass)
        del p.__dict__['collide_with']

We can't remove the particles during this loop because it causes problems if you remove an element from a list while iterating through that list. But now we can loop through the list of particles to remove and remove them from the list of particles in the universe. The way we're dealing with colliding particles, isn't perfect as slight inconsistencies will occur if one particle collides with more than one other in a single tick, but this should be relatively rare and the difference should be minor. We do however have to check that the particle is still in universe.particles in case it has been removed already.

for p in particles_to_remove:
    if p in universe.particles:
        universe.particles.remove(p)

Lastly, don't forget that the add combine() to the function list and then add it to our simulation.

'combine': (2, lambda p1, p2: combine(p1, p2)),

universe.addFunctions(['move', 'attract', 'combine'])

Now the simulation is complete. If you run it, you should see the cloud collapse as before, but larger particles will form, and as they do they will attract more and more particles towards them, thus getting ever bigger. Eventually you will probably end with a single large particle and maybe a couple of "planets" orbiting it. When I first made this simulation, I coloured particles over a certain mass yellow to indicate they had become a star.

Comments (8)

Anonymous on 24 Apr 2011, 2:23 a.m.

Sometimes I get this error:

yango@mmaker:~/python$ python estrella.py
Traceback (most recent call last):
File "estrella.py", line 44, in <module>
universe.particles.remove(p)

I haven't time to debug your code but this is a nice tutorial. Thanks.

Peter on 7 May 2011, 11:14 a.m.

Sorry about that. I'd fixed the problem in my program but forgot to update the tutorial. Thanks for pointing it out.

You need to change the code that removes particles to:

for p in particles_to_remove:
if p in universe.particles:
universe.particles.remove(p)

The reason for the error is that sometimes more than two particles will collide at once which will result in a particle being added to particles_to_remove more than once, but you can only remove it once.

I've update the tutorial to include this.

derek on 1 Aug 2012, 5:28 p.m.

The gravitational constant you use in the code and the above code snippet is 0.2, not 0.1 as mentioned in text, correct?

A more important question/observation, in the attract method in your PyParticles.py you have the following condition

if dist < self.size + self.size: return True

Does this mean that the gravitational attraction is not applied if the dist between particles is less that the radius of the particle? This is because once they get this close they will be combined? If I remove this condition from the attract method I don't see much (if any?) difference in behavior of the star formation simulation. Does this avoid or fix a bug?

Peter on 4 Aug 2012, 9:22 p.m.

Hi Derek,

I've updated the text to 0.2 which seems to be the constant I settle on for this simulation.

The line:
if dist < self.size + self.size: return True

Should be:
if dist < self.size + other.size: return True

I think the reason for it, is to stop attraction when particles should collide. I seem to recall that without, under some conditions (particularly if particles bounce off the walls), then the distance between two particles can become so small that the force becomes crazily high, or worse become 0, so simulation crashes. It's quite rare, but can happen if you run the simulation enough times.

Anonymous on 12 Dec 2018, 2:09 a.m.

It seems quite messy to seperately add the particle removing code in main along with the combining code for the particles. Can't the removal be done in combine too?

Peter on 16 Dec 2018, 8:43 a.m.

Yes, it's possible, but it's generally not a good idea to remove items from a list while you're iterating through that list.

Charlo on 12 Apr 2019, 11:03 p.m.

It's seems like PyParticles has no attribute for Environment. What should i change it to?

Peter on 13 Apr 2019, 2:44 p.m.

Hi Charlo,

PyParticles should have a Environment object if you've been following the previous tutorials. You can get the full code here; https://github.com/petercollingridge/code-for-blog/blob/master/pygame%20physics%20simulation/particle_tutorial_12/PyParticles.py

Pygame physics simulation