Thursday, 17th February 2011
Universal attraction: gravity
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. It should look something like this:
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.
The driving force behind particl' 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.
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 why 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:
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.
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 recalculate 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. Eventutally you will probably end with a single large particle and maybe a couple of 'planet's orbiting it. When I first made this simulation, I coloured particles over a certain mass yellow to indicate they had become a star.