Pygame physics simulation (Tutorial)

This series of tutorials demonstrates how to create a physical simulation using Python and Pygame. The tutorials start with the very basics and build up to a final simulation of a classical physics problem: the trajectory of a cannonball. Something like this:

If you're not interested in the mechanics of the program, you can skip ahead to tutorial 10 and download the PyParticles module and particle_tutorial_10.txt, which demonstrates how to use it from the bottom of the page.

Requirements

It would also be advantageous if you had a grounding in basic mathematics, specifically trigonometry. The tutorial includes some simple mathematical tricks that can be surprisingly powerful. These tricks took me a while to work out and I'll do my best to explain how how they work. If you want, you can blindly copy the code, but if you later want to change the behaviour of the simulation, it would be a good idea to understand how the mathematics works. If you're interested in learning more about the mathematics or the physics then I can recommend the Khan Academy, which has an amazing collection of educational videos.

You won’t however, require any prior knowledge of Pygame; this tutorial will start from the very basics, demonstrating how to use Pygame to create a graphical display, update it, and make it to respond to user inputs.

Uses for the simulation

Although this simulation is very simple, with a few adaptations, I have expanded for several, surprisingly diverse projects. For example:

Below is an example of a network of Chinese characters built using code from the tutorial (See here for more details).

All of these different projects used the same ideas, so should have resulted in a lot of code reuse, but instead, I seem to have spent a lot of time rewriting the code. This was partially because it was initially written in Java and later in Python, and partially because I've come up with better ways of doing things in each iteration. However, each time I rewrote the code, I spent a lot time frustrated, struggling with getting the maths right (witnessing all sort of odd behaviours from the particles in the meantime). So this time, I'm determined will be the last. I've tried to write the code in such a way that I can import the module and create networks easily. I've also decided to write this tutorial explaining how to create this similar sort of simulation in case anyone else wants to use it. I would like to put the program online for people to play with, but I'm not sure I can without redoing everything (again) in Java (or Canvas or SVG).

AttachmentSize
particle_tutorial_13.txt2.96 KB
PyParticles3.txt6.87 KB

Creating a Pygame window

In this first tutorial, we'll cover the basics and write some boilerplate code for Pygame - code you will almost always write whenever you use Pygame. We will:

  • Create a Pygame window
  • Stop the window from immediately disappearing
  • Close the window in response to a quit event
  • Change the title and background colour of the window

The complete code for this tutorial is shown at the bottom of this page.

The Pygame window

Since we are going to use the Pygame module, the first thing we need to do is import it.

import pygame

We can now create a Pygame window object (which I've called 'screen') using pygame.display.set_mode(). This object requires two values that define the width and height of the window. Rather than use constants, we'll define two variables, width and height, which will make the program easier to change later on. Feel free to use any integers that suit you. In order to display the window, we use the flip() function:

(width, height) = (300, 200)
screen = pygame.display.set_mode((width, height))
pygame.display.flip()

If you now run this program, you’ll see a 300 x 200 pixel window appear and then promptly disappear. The problem is that once as the flip() function has been called, the end of the code is reached, so the program ends. To keep the screen visible for as long as we want, we need to make sure the program doesn't end. We could do this by adding an infinite loop.

while True:
  pass

The problem with an infinite loop is that, by definition, it never ends. The program won't quit even if we want it to. If we try to close the window by clicking on the X, nothing happens. You have to use Ctrl + C in the command line to force the program to quit.

Closing the window

We want our window to persist until the user chooses to closes it. To acheive this, we monitor user inputs (known as 'events') using pygame.event.get(). This function returns a list of events which we can loop through and check to see whether any have the type QUIT. If we find such an event, we exit our loop, which is best done by changing a boolean variable (which I've called 'running').

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

The window now persists whilst 'running' is equal to True, which it will be until you close the window (by clicking the X). Note that if you use an IDE for Python programming, then it may interfere with Pygame. This isn’t normally a major problem but it can stop the Pygame window from closing properly. If so, adding pygame.quit() should solve the problem.

Changing the window's properties

Now we have a usable window, we can change its properties. For example, we can change its title using the set_caption() function.

pygame.display.set_caption('Tutorial 1')

We can change the background colour by filling the screen object. Colours are defined using a 3-tuple of integers from 0 to 255, for the red, green and blue values respectively. For example, white is (255,255,255). Changes need to be made before the flip() function is called.

background_colour = (255,255,255)
screen.fill(background_colour)

The final program

The complete program, after a bit of rearrangement, should now look like this:

import pygame

background_colour = (255,255,255)
(width, height) = (300, 200)

screen = pygame.display.set_mode((width, height))
pygame.display.set_caption('Tutorial 1')
screen.fill(background_colour)

pygame.display.flip()

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

Running the program should create a window that looks like this (on Windows XP):

A Pygame Window

It's not very exciting at the moment, but now we have a window that persists until we close it. In the next tutorial we'll draw some shapes in our window and start our simulation by creating a Particle object.

AttachmentSize
particle_tutorial_1.txt372 bytes

Drawing circles

In the previous tutorial we created a blank white screen; now we want to display some graphics on it. In this tutorial, we will:

  • Draw a circle in our window
  • Define a Particle class
  • Create and display a Particle object

This program builds on the code from the previous tutorial. The line numbers show where code should be added assuming the program is written exactly like the program attached in the previous tutorial (including blank lines). In this tutorial and all the following tutorials, the final code is available at the bottom of the page as a txt file.

Drawing with Pygame

In Pygame there are various functions for drawing simple shapes. We will display our particles as circles so use the pygame.draw.circle() function. This function takes several parameters:

  • A surface where the circle will be drawn (here it's our screen)
  • A colour
  • An (x, y) coordinate
  • A radius
  • A thickness (optional)

For example.

pygame.draw.circle(screen, (0,0,255), (150, 50), 15, 1)

This draws a blue ((0,0,255) is blue in RGB notation) circle centred at (150, 50), with radius 15 and thickness of 1 pixel. Note that this function must be called after the screen.fill(background_colour) command and before the flip() command. If you run the program now you should see something like this (I also changed the title of window, but it's not important):

Tutorial 2 screen shot

If you’re not familiar with how computer displays work, you might have expected the circle to be near the bottom of the screen, since the particle's y-coordinate is 50 and the screen is 200 units high. However, on a computer display, the origin (0, 0) is in the top, left of the screen, with the x-axis increasing from left-to-right and the y-axis increasing from top-to-bottom; the circle is therefore centred 50 pixels down from the top of the screen.

Coordinates on a computer display

The Particle object

In our final simulation we will want several particles, each of which will have the same type of attributes. Thus it makes sense to define a Particle class. When we create each Particle object with the __init__ function (note the double underscores which indicates a built-in method), we will give it an x,y coordinate and a size. The colour and line thickness will be given default values. We'll add more attributes as we go along.

class Particle:
  def __init__(self, (x, y), size):
    self.x = x
    self.y = y
    self.size = size
    self.colour = (0, 0, 255)
    self.thickness = 1

Next we add a display() method to our Particle object that will draw it (as a circle) on the screen.

def display(self):
  pygame.draw.circle(screen, self.colour, (self.x, self.y), self.size, self.thickness)

Now we can remove the pygame.draw.circle() call and replace it by code that creates a particle object (defining its x, y coordinates and size) then calling its display() method. If you're not familiar with object-orientated programming, this may seem more long winded, but it will make building on our program a lot easier later on.

my_first_particle = Particle((150, 50), 15)
my_first_particle.display()

If you run the program now, it should look exactly the same as before. However, in the next tutorial I'll show you how the code now can be easily changed to display many more circles. I'll also introduce the random module which is very handy when creating simulations.

AttachmentSize
particle_tutorial_2.txt755 bytes

Randomness

Most simulations will start with a random (or semi-random) initial configuration. Some will introduce random noise to mimic real life uncertainty. In this tutorial, we will:

  • Use Python’s random module
  • Create multiple Particle objects

While this tutorial isn't specific to Pygame, it is import for our simulation. If you're already familiar with Python's random module then you can skip this tutorial. As before, the final code is at the bottom of this page.

The random module

To generate random numbers in Python (not just Pygame), we need to import the random module. This module is part of the main Python program so you don't need to download anything extra.

import random

We now have access to various functions that generate random numbers. The functions I use most often are:

  • random() - returns a random number between 0.0 and 1.0
  • uniform(a, b) - returns a random number between a and b
  • randint(a, b) - returns a random integer between a and b
  • choice(list) - returns a random element from the list
  • shuffle(list) - rearranges the order of the list

For the full range of functions see the Python documentation. For now we just need random.randint().

The code below demonstrates how to create a circle with a radius between 10 and 20 pixels (inclusive) and with a random position on the screen. To ensure that the no part of the circle is off-screen rather than picking a value between 0 and the window's dimensions, we reduce the range of values by the size of circle.

size = random.randint(10, 20)
x = random.randint(size, width - size)
y = random.randint(size, height - size)
my_random_particle = Particle((x, y), size)
my_random_particle.display()

Now, every time the program is run, unless we are very lucky, it will display a circle with a different size and position.

Creating multiple particles

Since we're using variables instead of constants for the particle's position, we can easily add multiple particles to the screen. This will be useful later because most simulations of particles are likely to have more than one. We deal with multiple Particle objects by creating an array and filling it with random Particles using a for loop. Replace the above code with the following:

number_of_particles = 10
my_particles = []

for n in range(number_of_particles):
    size = random.randint(10, 20)
    x = random.randint(size, width-size)
    y = random.randint(size, height-size)
    my_particles.append(Particle((x, y), size))

To display each of the particles we need to change the line:

my_random_particle.display()

To:

for particle in my_particles:
    particle.display()

For the full program (rearranged to be slightly more logical) see the file at the bottom of this post. If you run the program, you will see ten randomly placed and randomly sized circles, some of which may overlap. You might want to avoid having the circles overlap, but I’ll leave you to work out how to do that.

Tutorial 3 screen shot

In the next tutorial, the simulation will become a bit more interesting as things get dynamic and the particles move.

AttachmentSize
particle_tutorial_3.txt1011 bytes

Movement

Pretty much by definition, all simulations have values which change over time. One value that often changes is the position of an object, in which case the object is moving. In this tutorial, we will:

  • Give our particles speed and direction (a velocity vector)
  • Use basic trigonometry to convert vectors into movement
  • Make the particles move randomly

Our simulation is a discrete time simulation, which means that we split time into individual steps. At each step we update the simulation a bit then display the new situation. We keep updating the simulation until the user exits. We do this by adding code into the infinite while loop that we've already created. The first thing we do therefore is move the following block of code from before the while loop to inside it. This will have no effect on how the program runs but allows us to add additional function calls later.

for particle in my_particles:
    particle.display()
pygame.display.flip()

Representing movement

The simplest way to represent movement is to create two attributes: dx and dy. Then during each pass through the loop, add dx to x and dy to y. So, for example, if a particle has dx=2 and dy=1, it will follow a shallow diagonal, from left to right and top to bottom. This s the method I used for a long time - it is simple and fine for many situations. However, when we come to work out interactions between two particles, things become a lot more complex.

My preferred method now is to create attributes to represent speed and direction (i.e. velocity). This requires a bit more work to start with but makes life a lot easier later one. This method is also good for creating objects that have a constant speed, but varying direction. I actually first started using this method when trying to create an ant simulation. The ants have a creepily realistic movement when given a constant speed and randomly changing their a direction every few seconds.

So let's give our particle object a speed and a direction.

self.speed = 0.01
self.angle = 0

The math module

Since we’re going to use some trigonometry, we need to import Python’s math module. Like the random module, this is part of main Python program so there's no need to download anything extra.

import math

We now have access to various mathematical functions, such as sine and cosine, which we'll use shortly.

Movement vectors

We now need to add a move() function to our particle object. This function will change the x, y coordinates of the particle based on its speed and the angle of its movement. The diagram below illustrates how we calculate the change in x and y. I find it simplest to consider an an angle of 0 to be pointing upwards, despite the fact that y-axis actually pointing downwards in computer graphics.

The mathematics of movement

To calculate the change in x and y, we use some secondary school-level trigonometry as shown in the code below. Remember to minus the y to take into account the downward pointing y-axis. Although it doesn't make much difference at the moment you will have to be consistent with your signs later.

def move(self):
    self.x += math.sin(self.angle) * self.speed
    self.y -= math.cos(self.angle) * self.speed

Another point to bear in mind is that the angles are all in radians. If you’re used to working with degrees then this might be a little confusing; just remember that 1 radians = 180°/π. Therefore if you want to make the particle move forwards (left to right) along the screen, then its angle should be π/2 (90°). In Python we do this like so:

self.angle = math.pi/2

Now we can we call the particles' move() function during the loop immediately before calling their display() function.

for particle in my_particles:
    particle.move()
    particle.display()

If we run this program now, we’ll see the circles moving right, leaving a smear across the screen:

Smeared circles

The reason the particles smear is that once Pygame has drawn something on the window it will stay there unless drawn over. The easiest way to clear the circles from the previous time step is to fill the whole screen with the background colour. We can do this by simply moving the screen.fill(background_colour) command into the loop. The effect of movement is therefore achieved by drawing a particle, then clearing it and drawing a short distance away.

You might have also got a deprecation warning telling you "integer argument expected, got float" followed by the pygame.draw.circle() function. The reason, as you may have guessed, is that this pygame can only draw circles at integer x, y coordinates and after we move the circles, their coordinates become floating point numbers. Although the Python deals with the problem perfectly well, we should convert the x, y coordinates to integers first:

pygame.draw.circle(screen, self.colour, (int(self.x), int(self.y)), self.size, self.thickness)

Random movement

If you run this program now, you’ll see all the circles moving rightward at the same speed, so it will look like they are all draw on a single moving surface. We can making things more interesting by giving each particle a random speed and direction. We could do this by defining the speed and angle attributes as random in the Particle class, but I prefer to set these values (the default behaviour) to 0. We can update them once the individual objects have been created, by altering our particle-creating loop:

for n in range(number_of_particles):
    size = random.randint(10, 20)
    x = random.randint(size, width-size)
    y = random.randint(size, height-size)

    particle = Particle((x, y), size)
    particle.speed = random.random()
    particle.angle = random.uniform(0, math.pi*2)

    my_particles.append(particle)

Notice that we use the random.random() function to generate a speed between 0 and 1, and the random.uniform() function to create a random angle between 0 and 2*pi, which covers all angles. If you run the program now, the circles will fly off the screen at random angles and speeds. In the next tutorial I tell you how to keep the particle within the bounds of the window.

AttachmentSize
particle_tutorial_4.txt1.34 KB

Boundaries

When creating a simulation you will often want to make it bounded, to avoid simulating an infinite region. In this simulation, a realistic, physical solution is to add virtual walls that particles bounce off. In this tutorial, we will:

  • Test whether particles have moved out of the window
  • Move such particles to where they would have bounced
  • Change their direction appropriately

Another solution would be to bend a 2D simulation into virtual torus (ring doughnut shape). We do this my making particles that leave one side of the simulation appear on the opposite side.

Exceeding boundaries

The first thing our bounce() function needs to do is test whether a particle has gone past a boundary. The four boundaries are at:

  • x=0 (the left wall)
  • x=width (the right wall)
  • y=0 (the ceiling)
  • y=height (the floor)

Since this simulation has discrete steps of time, we unlikely to catch a particle at the exactly point that it 'hits' a boundary, but rather at a point when it has travelled a bit beyond the boundary. If the speed of the particle is low, then the particle is unlikely to have gone much beyond the boundary (the maximum distance it will have exceeded the boundaries is, in fact, equal to its speed). We could chose to ignore the discrepancy and simply (especially if our simulation is becoming too computationally intense), and simply reflect the particle angle, but we might as well be accurate for now.

Therefore, if a particle exceeds a boundary, we first calculate by how much it has exceeded it (which I'll call d). For example, if the particle's x value is greater than the window's width minus the particle's size (it has gone beyond the right wall), then we calculate d:

d = self.x - (width - self.size)

We then reflect the position in the boundary (i.e. bounce it) by setting the x coordinate as below:

self.x = (width - self.size) - d

This simplifies to:

self.x = 2*(width - self.size) - self.x

So we don't actually need to calculate variable d. Replacing width with 0 calculates the x position when the particle crosses the left wall. The y coordinates can be calculated in a similar way.

Graph of bouncing

The most important feature of bouncing is reflecting its angle in the boundary, which is where our use of vectors starts to become useful (though its more useful when the boundaries aren't straight). A vertical boundary is has an angle of 0, while a horizontal boundary has an angle of pi. We therefore reflect vertically bouncing particles by subtracting their current angle from 0 and from pi in the case of horizontal bouncing.

The final bounce function should look like this:

def bounce(self):
    if self.x > width - self.size:
        self.x = 2*(width - self.size) - self.x
        self.angle = - self.angle

    elif self.x < self.size:
        self.x = 2*self.size - self.x
        self.angle = - self.angle

    if self.y > height - self.size:
        self.y = 2*(height - self.size) - self.y
        self.angle = math.pi - self.angle

    elif self.y < self.size:
        self.y = 2*self.size - self.y
        self.angle = math.pi - self.angle

Finally, don't forget to call the bounce() function after calling the move() function and before calling the display() function.

    for particle in my_particles:
        particle.move()
        particle.bounce()
        particle.display()

Before running the simulation, I would also increase the width and height to 400 to make things easier to see. When you run this simulation, you should now see ten particles bouncing around at different speeds, but completely ignoring one another. In a later tutorial, we’ll allow the particles to interact with one another.

AttachmentSize
particle_tutorial_5.txt1.9 KB

Gravity

By the end of this tutorial we will have a complete simulation of a classical problem in physics: the trajectory of a projectile. In this tutorial, we will:

  • Create a function to add vectors
  • Create a vector to represent gravity
  • Add drag and elasticity

And after that we should have a simulation that looks a bit like this:

As usual the code is available at the bottom of this article.

Gravity works by exerting an constant downward force to each of the particles in the simulation. If we had stored the particles' movement as two positions, dx and dy (as discussed here), then we could simply add a constant to the dy value. However, since we’re using vectors, it’s a bit more complex because we need to add two vectors. Once we have a function to add vectors (which took me a while to work out, but is probably the most useful thing in these tutorials), everything else will be a lot easier.

Adding Vectors

The addVectors() function takes two vectors (each an angle and a length), and returns single, combined a vector. First we move along one vector, then along the other to find the x,y coordinates for where the particle would end up (labelled (x, y) on the diagram below).

def addVectors((angle1, length1), (angle2, length2)):
    x  = math.sin(angle1) * length1 + math.sin(angle2) * length2
    y  = math.cos(angle1) * length1 + math.cos(angle2) * length2

We then calculate the vector that gets there directly. To do this we construct a right-angle triangle as shown in the image below. Then we use good old trigonometry and a couple of handy functions from Python’s math module, which I’ve only recently discovered. The new vector length (speed of the particle) is equal to the hypotenuse of the triangle, which can be calculated using math.hypot(). This takes an x,y coordinate and calculates its distance from the origin (0,0). Note that while the position of our particle on the screen is not (0,0), all the vectors are relative to the particle's position, so can be considered to begin at 0,0.

Vector addition

The angle of the new vector is slightly more complex to calculate. First, we find the angle in the triangle, by calculating the arctangent of y/x. We could do this using the math.atan() function but then we would then need to deal with the case of x=0 and work out the sign of angle. However, Python provides us with a handy function math.atan2(), which takes the x, y coordinates, works out the sign of the angle for us and behaves correctly when x=0. Once we have the angle of the triangle, we subtract it from pi/2 to calculate the angle of the vector.

length = math.hypot(x, y)
angle = 0.5 * math.pi - math.atan2(y, x)
return (angle, length)

Gravity

Now we can create a vector for gravity: the angle is pi, which is downward and I chose a magnitude of 0.002 purely by experimentation. Feel free to change it:

gravity = (math.pi, 0.002)

Then, in the Particle's move() function, we add the gravity vector to the particle’s vector:

(self.angle, self.speed) = addVectors((self.angle, self.speed), gravity)

Now when we run the simulation, we see the particles bouncing about on the ground endlessly. You may wish to set the number_of_particles to 1 to make it easier to follow a particle.

Friction

To complete the trajectory simulation and stop the particles from bouncing forever we need to add two more physical phenomena: drag and elasticity.

Drag represents the loss of speed particles experience as they move through the air - the faster a particle is moving, the more speed is lost. I find it simpler (and computationally quicker) to define a drag variable that represents the inverse of drag. We multiple a particle's speed by this value at each time unit, thus the smaller the value, the more speed is lost. Elasticity represents the loss of speed a particle experiences when it hits a boundary.

You can play with the values to see what looks reasonable, though both should be between 0 and 1. I found that 0.999 and 0.7 respectively work quite well.

drag = 0.999
elasticity = 0.75

To the Particle's move() function add:

self.speed *= drag

Another option would be to multiple by a factor that varied inversely with the particle's size (e.g. self.speed *= (1-self.size/10000.0)), but I don't think that's necessary.

After each of the four boundary conditions, add:

self.speed *= elasticity
self.speed *= elasticity
self.speed *= elasticity
self.speed *= elasticity

And there you have it: a complete simulation of a projectile's trajectory. And we haven't had to explicitly solve any of the equations of motion. Try running the simulation several times to see what happens. You might find it easier to set the y coordinate to 20, so the particle always starts at the top of the simulation. In the next tutorial, we add some user interaction so you can pick up, drop and throw the ball (particle).

AttachmentSize
particle_tutorial_6.txt2.51 KB

Mouse interactions

In the previous tutorial, we completed a simulation of a particle's trajectory. Unfortunately, once we'd added drag and friction, the simulation came to a halt pretty quickly. In order to try out different starting conditions, we have to edit the code and restart the simulation. It would be much simpler if we could reach into our virtual world and interact with it directly. In this tutorial, we will:

  • Add the ability to select particles with a mouse click
  • Add the ability to move particles with the mouse
  • Add the ability to throw particles with the mouse

As usual the complete code can be found by following the link at the bottom of this post. By the end of this tutorial, you too will be able to fling the particles about like this:

Getting the mouse position

Firstly, set the number_of_particles to 3, which will give us a few to chose between without cluttering the screen.

In order to test whether the user has selected a particle, we need to know where they have clicked and we do this by using pygame.mouse.get_pos(). As the name suggests, this function returns the position of the mouse in the Pygame window as an x, y coordinate. For now we only need to know the position of the mouse when the user clicks. We test this by monitoring Pygame events as we did in tutorial 1 using pygame.event.get(). Since we are already calling this function, looking for QUIT events, we might as well search for MOUSEBUTTONDOWN events at the same time.

for event in pygame.event.get():
    if event.type == pygame.QUIT:
        running = False
    if event.type == pygame.MOUSEBUTTONDOWN:
        (mouseX, mouseY) = pygame.mouse.get_pos()
        print mouseX, mouseY

Now, whenever we click in the Pygame window, the position of the mouse is printed to the command line. Note that the MOUSEBUTTONDOWN event is only found once per mouse click; it is not repeatedly found if the mouse is held down.

Selecting particles

Rather than print the position of the mouse, we want to test whether the mouse is within the boundary of a particle. Since the particles are circular, this is very straightforward: we test all the particles to see whether the distance from the mouse to particle's centre is less than that particle's size. We can do this with our old friend, math.hypot().

def findParticle(particles, x, y):
    for p in particles:
        if math.hypot(p.x-x, p.y-y) <= p.size:
            return p
    return None

The return None command is not strictly necessary as Python will return None if it reaches the end of a function without finding a return command. Also note, that if the function finds that the first particle in the array is selected, it will return that particle without bothering to check the other particles.

We can test our findParticle() function by replacing the print function with:

selected_particle = findParticle(my_particles, mouseX, mouseY)
if selected_particle:
    selected_particle.colour = (255,0,0)

If you click on a particle now, its colour should change to red. You can remove this code now if you want.

When we have selected a particle, we want it to stop moving, so change the code that makes particle move to:

for particle in my_particles:
    if particle != selected_particle:
        particle.move()
        particle.bounce()
    particle.display()

And add selected_particle = None before the while loop starts so the variable exists before a mouse click. We also need to cancel the selection once the mouse button is released, so below the code checking for MOUSEBUTTONDOWN events, add:

elif event.type == pygame.MOUSEBUTTONUP:
    selected_particle = None

Dragging particles

We can now select a particle and stop it, but we want to be able to drag it somewhere. We do this by making the x, y coordinates of a selected particle (if there is one) equal the coordinates of the mouse.

if selected_particle:
    (mouseX, mouseY) = pygame.mouse.get_pos()
    selected_particle.x = mouseX
    selected_particle.y = mouseY

We can now pick up and drag the particles. However, when we let go, the particles drop to the ground in a disappointing sort of way. Instinctively, we expect the particle to be released with a speed equal to the speed the mouse was moving at the time. For this we need to measure the difference between the position of the mouse and the particle and creating a vector that joins them. We can then setting the particle's angle and speed to be this vector. Replace, the above code with the below:

if selected_particle:
   (mouseX, mouseY) = pygame.mouse.get_pos()
    dx = mouseX - selected_particle.x
    dy = mouseY - selected_particle.y
    selected_particle.angle = math.atan2(dy, dx) + 0.5*math.pi
    selected_particle.speed = math.hypot(dx, dy) * 0.1

I've found that setting the particle's speed to 0.1 time the actual length of the vector works best. This means it will actually take ten units of time for the particle to catch up with the mouse. Note also that the angle is the arctangent plus half pi. Just trust me it is. I'll try to draw a diagram explaining why later. We now need to ensure the selected particle moves, so change the code back to move all the particle each turn (sorry).

So now we have a simulation that we can interact with. It's quite fun to fling the particles about and I think this simulation could easily be adapted into some sort of game. However, the particles still don't interact with one another. In the next tutorial, I hope to remedy this.

AttachmentSize
particle_tutorial_7.txt3.23 KB

Collisions

In the previous tutorial we allowed the user to interact with the particles, but the particles didn't interact with one another. Now it's time to make them more tangible. In this tutorial, we will:

  • Test whether any two particles have collided
  • Make particles that have collided bounce
  • Prevent colliding particles from sticking to one another

Things start to get complex once we have multiple particles interacting with one another. In fact, it's mathematically impossible to solve the equations describing three or more interacting objects. In our simulation we'll have to make some simplifications, which we make it 'inaccurate', however, it should still represent a reasonable approximation to reality.

Collision testing

Firstly, we want to check whether any two particles overlap with one other, so we need to compare every particle with every other of particle. We can do this with a nested for loop. Since we already have one loop going through the particle array, we can use that.

for i, particle in enumerate(my_particles):
    particle.move()
    particle.bounce()
    for particle2 in my_particles[i+1:]:
        collide(particle, particle2)
    particle.display()

Note that the second loop is not a full loop; if we had two full loops, we'd compare every particle to every other particle twice over (and compare each particle to itself). Instead, we compare each particle with every particle with a higher index than it in the array. We therefore need to know the index of the current particle which we can get using enumerate. We use this value (i) to slice the my_particle array from i+1 to the end to construct the second loop.

Now we need to actually define the collide() function. The first thing the function has to do is to test whether the two particles have collided. This is very simple as the particles are circles. We simply to measure the distance between them (using math.hypot() again) and test whether this value is less than their combined radius..

def collide(p1, p2):
    dx = p1.x - p2.x
    dy = p1.y - p2.y
    
    distance = math.hypot(dx, dy)
    if distance < p1.size + p2.size:
        print 'Bang!'

So far all that happens when two particles collide is that we print 'Bang!'. I'll complete the function as soon as I get a chance.

The angle of collision

When two particles collide, we want them bounce off each other. Theoretically, when two circular particles collide they contact at an infinitely small point. The angle of this point is the tangent of the particle at this point. As the diagram below I hope shows, this angle is perpendicular to the line that joins the centre of the two particles.

The mathematics of colliding particles

We can treat the collision as though the particles were bouncing off a flat surface with the angle of the tangent. We can find the angle of the tangent with the following code:

tangent = math.atan2(dy, dx)

To reflect the angle of the particles on this surface, we subtract their current angle from two times the tangent. I'll have to explain the logic behind this at a later date. It's at this point that knowing the angle that our particle is travelling starts to become useful. 

p1.angle = 2*tangent - p1.angle
p2.angle = 2*tangent - p2.angle

Then we need to exchange the speeds of the two particles as they transfer their energy to on another. We can do this in a single line by constructing a tuple.

(p1.speed, p2.speed) = (p2.speed, p1.speed)

Note that writing the exchange as two lines as below does not work.

p1.speed = p2.speed
p2.speed = p1.speed

The reason this doesn't work, as you may have realised, is that when we come to assign p2.speed on the second line, we are assigning it to the new value of p1.speed, which we have already made equal to p2.speed. Constructing an tuple, which is immutable, allows us to avoid this problem. I would then add some code to reduce the energy of both particles as a result of the collision:

p1.speed *= elasticity
p2.speed *= elasticity

Sticky problem

If you run the simulation now, you'll probably find that the particles stick to one another, and can even get stuck in mid air. The reason for this is that we have a discrete simulation (i.e. the time between updating the particle positions is not infinitely small). This introduces errors, which, for the most part, aren't important (though you should be aware of them).

What's happening is that when a particle collision is detected, the particles have actually overlapped slightly. When we alter their angle and speed then move them, we can't be sure that the particles are no longer overlapping (particularly because of drag, gravity and elasticity also altering their movement). If they still overlap then they will 'bounce' back the way they came, towards each other. They can then get trapped in a loop of 'internal bouncing', which gradually reduces their speed to nothing.

To avoid this problem we add some code into the collide() function to fudge the fact that the particles should have bounced before they overlapped. We calculate the angle between the two particles (which is 90 degrees to the tangent) and move the particles along this angle one pixel. We could work out exactly how far to move the particles along this angle, but I don't think it's worth the extra calculation. The code below moves the particles away from one another by a sufficiently large amount, and seems to work pretty well. 

angle = 0.5 * math.pi + tangent
p1.x += math.sin(angle)
p1.y -= math.cos(angle)
p2.x -= math.sin(angle)
p2.y += math.cos(angle)

Now the particles should bounce off one another cleanly:

AttachmentSize
particle_tutorial_8.txt3.92 KB

Mass

In the previous tutorial we created a reasonable simulation of billiard balls or projectiles bouncing around. In this tutorial we'll add a couple more physical attributes, most importantly, mass. In this tutorial, we will:

  • Give our particles mass
  • Give particles a range of colours based on their density
  • Alter collisions to take into account mass
  • Alter drag to take into account particle size and mass

As usual the code for the completed simulation at the bottom of the this post.

Mass

Currently, when two particles collide, they swap speeds which is reasonable approximation for two particle of equal mass, but what if we want to have particles with different masses?Firstly, we need to add a mass attribute to our Particle class. We want to be able to define the mass attribute when we create our object, but it is useful to have a default value (which I've arbitrarily set to 1). Not only does this mean that we don't have to set the mass if we don't want to, it also makes our Particle class backwards compatible: our old code can still create valid particles using this new class without raising an error.

class Particle():
    def __init__(self, (x, y), size, mass=1):
        self.x = x
        self.y = y
        self.size = size
        self.mass = mass

If you run the program now you should see no difference from before, and you also shouldn't get any errors - the particles are assigned a default mass of 1.

Now lets give our particles a random range of masses. Rather that generate the particle's mass directly, I decided to generate a value to represent the particle's density. We can then use a particle's density to generate its mass. We do this by multiplying the particle's density by its size squared (since size refers to the particle's radius and this is a 2D simulation). The advantage of this method is that we have a large range of masses (from 100 to 8000), but the distribution favours medium masses.

for n in range(number_of_particles):
    size = random.randint(10, 20)
    density = random.randint(1, 20)
    x = random.randint(size, width-size)
    y = random.randint(size, height-size)

    particle = Particle((x, y), size, density*size**2)

To allow us to see differences in density, we can colour the particles with a range of colours. We could create a dictionary that converts a density into a colour (e.g. a density of 1 is red, a density of 2 is yellow), but that approach is a bit long-winded and not very adaptable. Instead, we can directly convert the density value into a colour intensity. The code below will colour particles with a density of 1 a very pale blue (180, 180, 255) and particles with a density of 20, an intense blue (0, 0, 255). Particles with intermediate densities will have intermediate shades of blue. We can get an idea of the mass of a particle by how intense its colour is and its size.

particle.colour = (200-density*10, 200-density*10, 255)

Note that despite generating a range of blue shades, we keep blue at the maximum intensity of 255, and vary of the red and green values.

Force

We can now see differences in particle density by their colour, but not by their behaviour. The behaviour we want is that, when a heavy particle collides with a lighter particle, the lighter one should be more perturbed than the heavier.  I have to admit that is was beyond my physics ability and I had to check Wikipedia. Translating these equations into code, we can change our collide function to:

    if dist < p1.size + p2.size:
        angle = math.atan2(dy, dx) + 0.5 * math.pi
        total_mass = p1.mass + p2.mass

        (p1.angle, p1.speed) = addVectors((p1.angle, p1.speed*(p1.mass-p2.mass)/total_mass), (angle, 2*p2.speed*p2.mass/total_mass))
        (p2.angle, p2.speed) = addVectors((p2.angle, p2.speed*(p2.mass-p1.mass)/total_mass), (angle+math.pi, 2*p1.speed*p1.mass/total_mass))
	p1.speed *= elasticity
	p2.speed *= elasticity

The code may look complex, but it is actually shorter than before. Rather than reflecting the angle the of particles, we assign the particles a new angle and speed in a single step. We are basically taking a particle's vector and adding a second whose angle is perpendicularly to the angle of collision, and whose magnitude is based the momentum (mass x velocity) of the second particle.

We can also now, improve the code that prevents particles from overlapping, though it's not that important:

overlap = 0.5*(p1.size + p2.size - dist+1)
p1.x += math.sin(angle)*overlap
p1.y -= math.cos(angle)*overlap
p2.x -= math.sin(angle)*overlap
p2.y += math.cos(angle)*overlap

Drag

Finally, since our particles now have different masses, we should update how they are affected by drag, or air resistance. Theoretically we could now model air resistance directly by adding lots of tiny air particles, but that would require a lot of computational power. Instead, we can calculate the approximate effect of millions of tiny particles colliding with each particle.

From the collision equation, we can see that the loss of speed due to a collision is approximately proportional to:

1 / (mass of particle + mass of colliding object)

Since the number of collisions with air particles in a unit of time will be proportional to the distance the particle has travelled in a unit of time, every update, we should make:

particle.speed = particle.speed - particle.speed * (1 / (particle.mass + mass_of_air))

Which can be rewritten:

particle.speed *= particle.mass / (particle.mass + mass_of_air)

However, the number of collisions with air particles will also depend on the size of the particle. If we imagine that the above equation represents the loss of speed for a particle of size 1, then a particle of size X will experience that loss of speed X times over, which gives us:

particle.speed *= (particle.mass / (particle.mass + mass_of_air)) ** X

We now have an equation for the drag of a given particle, which we can add to its __init__() function:

self.drag = (self.mass/(self.mass + mass_of_air)) ** self.size

Now we just need to define a constant called mass_of_air, which represents the mass (or momentum) or air particles that hit a particle of size 1, travelling at speed 1 in 1 time unit. I just picked a number that gave a realistic-looking simulation. This variable can replace the variable for drag, which we no longer require.

mass_of_air = 0.2

If you change this value to 2.0, the simulation will appear to be running in syrup; if you change the value to 0, there will be no air resistance. To actually implement the next system of drag, change the line in the move() function from:

self.speed *= drag

To:

self.speed *= self.drag

Before running the simulation, I think its a good idea to 'turn off' gravity by commenting out the code:

(self.angle, self.speed) = addVectors((self.angle, self.speed), gravity)

The simulation is now like looking down on a table of billiard balls, and should look something like this:

Now the basic simulation is complete, we can use it in various ways as described in the introduction.

AttachmentSize
particle_tutorial_9.txt4.19 KB

Making a module

In the previous tutorials we built up a reasonable physics simulation. Now we can take using this simulation in several different directions. Rather than copy large chunks of code every time we make new simulation, it makes sense to create a separate module that can be used by all of them. In this tutorial we will:

  • Convert our simulation into a module
  • Rewrite our particle simulation using this module
  • Introduce keyword-arguments to create flexible code

This tutorial is different from previous ones in that it is not specific to Pygame. In fact, it's not really specific to Python either: it's a general discussion on how to arrange code to maximise its reuse and flexibility. It's more philosophical (or wishy-washy, depending on your point of view), than the previous tutorials, but it's something I've had to think about on multiple occasions. I thought it might be helpful to write about what I've learnt and the thought process I go through when writing a module, though I'm far from being an expert. At the very least, it's forced me to justify my reasons for arranging the code in the way I have.

Splitting code

In this tutorial we are going to create two files: a PyParticles module and a python program that uses this module. I generally find it a good idea when making a module to make a separate program to test it. In fact it's better to make program that tests it more sytematically (e.g. unit tests), but for now we'll just make sure we can recreate our previous simulation. Since we're not actually going to write much new code, just rearrange it into two files, it's probably easiest to download the code from the bottom of this post to see the changes.

There's nothing particularly special about Python modules - they are no different from any other Python program, though generally they only contain classes and functions (they may have additional code which is used when the file is run not as a module). When creating a module, we need to think about how we expect the code to be used. On one hand we want to include as much code as we can in our module to maximise reuse, but on the other hand we want to make the code as flexible as possible, which might mean leaving out some aspects so they can be defined outside of the module.

We'll start by creating a file called PyParticles.py and adding the Particle object and the addVector(), findParticle() and collide() functions. We'll also need to import the math and random modules. But what about Pygame?

Separating particles from their display

At first glance it might seem obvious to include the code for importing, initialising and creating the Pygame screen: surely we'll always require it. But what if we (or someone else) later wants to display the particles using a different package, like PyOpenGL, or even create a command line display? Or what if you just want some code to tell you where a particle with a given position, mass and velocity will be in t units of time? More likely, what if you want to incorporate this code into a game in which you have already initialised and created a Pygame screen and now want to add the simulation to part of the screen.

In general it's a good idea to separate code that determine the behaviour of your objects (the Model) from code that determines how your objects are displayed (the View); it makes changing either one independently of the other much simpler. This is part of the Model-View-Controller architecture, only we're not going to separate the Controller (i.e. the mouse input) for now.

Our test program will therefore include all the Pygame code, so we'll start by creating a skeleton Pygame program in a new file:

import pygame
import PyParticles

pygame.display.set_caption('Tutorial 10')
(width, height) = (400, 400)
screen = pygame.display.set_mode((width, height))

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

    pygame.display.flip()

We now need to remove the display() function from the Particle object. We'll move this code into the Pygame while running loop later. We can leave the Particle's colour attribute even though it is a display property because it makes sense to associate it with the particle object.

Making a container

In order to display the particles, we need to loop through a list of all the particles, but where best to put this list? It makes sense to put it in the PyParticles module, and I think it's best to put the list within a container: an object that stores all the particles and simulation parameters (such as friction and elasticity). The advantage of this approach is that it makes it easier to change parameters and to run multiple simulations at the same time if we so wish.

Below is an Environment class which will allow us to create instances of our simulation. You could even create multiple environments and display them all on the same screen. The Environment class contains the global variables from the old program, including the list of particles (note that background_colour has become simply self.colour). The idea is that when we create a simulation, we import PyParticles, create an Environment object, add particles to it and set it running.

class Environment:
    def __init__(self, (width, height)):
        self.width = width
        self.height = height
        self.particles = []

        self.colour = (255,255,255)
        self.mass_of_air = 0.2
        self.elasticity = 0.75

Now let's go back to our skeleton Pygame program to test our new code.

env = PyParticles.Environment((width, height))

This code creates an Environment object called 'env' with the same dimension as the Pygame screen. If we wanted we could use different dimension, which would allow to add other elements to the screen.

Adding particles to the environment

Now we've set up the environment, we want to add some particles to it. For this we add an addParticles() function to our Environment class. To decide how we want this function to behave, I find it helps to imagine how we want to call it. For example, we probably want be easily add, say, 5 random particles with a call like:

env.addParticles(5)

But what if we want our particles to all be the same size? It would nice if we could get size particles of size 10 by simply writing:

env.addParticles(5, size=10)

And what if we want more control? For example, to follow the trajectory of a particle with specific properties. Ideally we just want to write as little as possible:

env.addParticles(x=5, y=10, speed=1.5, angle=0.8, size=4, mass=1000)

This means using keyword arguments: a list with an arbitrary number and type of parameters. To the Environment class, add:

def addParticles(self, n=1, **kargs):
  for i in range(n):
    size = kargs.get('size', random.randint(10, 20))
    mass = kargs.get('mass', random.randint(100, 10000))
    x = kargs.get('x', random.uniform(size, self.width-size))
    y = kargs.get('y', random.uniform(size, self.height-size))

    p = Particle((x, y), size, mass)
    p.speed = kargs.get('speed', random.random())
    p.angle = kargs.get('angle', random.uniform(0, math.pi*2))
    p.colour = kargs.get('colour', (0, 0, 255))
    p.drag = (p.mass/(p.mass + self.mass_of_air)) ** p.size

    self.particles.append(p)

To create a function that uses keyword arguments, we add **kargs to the end of the parameter list. This will generate a dictionary called kargs of any passed keywords. We can then use get(keyword, default) to search the dictionary for specific keywords, and use a default value if it's not present. Note that you'll now need to import the math and random module into this module to get the default random values. I tried to pick default values that I thought would be most useful.

I also made, n, the number of particles to be added optional, so if you leave it out the function call, the default behaviour is to create one particle.

Displaying particles

Now we can go back to our skeleton Pygame program and to add code to display the particles. It's pretty simple, we just add the following with the loop:

screen.fill(env.colour)
    
for p in env.particles:
    pygame.draw.circle(screen, p.colour, (int(p.x), int(p.y)), p.size, p.thickness)

This fills the screen with the colour of the environment (currently set to white), then loops through the environment's particles, drawing circles to represent them. Now we just need to get the particle to move. Ideally, we'd like to do this with a simple call:

env.update()

Updating the simulation

The environment update method is essentially the remaining code from the main loop in tutorial 9. To the Environment class add:

def update(self):
    for i, particle in enumerate(self.particles):
        particle.move()
        self.bounce(particle)
        for particle2 in self.particles[i+1:]:
            collide(particle, particle2)

You'll notice a couple of difference. First, we're not calling a display method since we've dealt with that elsewhere. Second, we've changed particle.bounce() to self.bounce(particle). The reason is that it makes sense to make bounce a method of the environment, since it depends on where the boundaries of the environment are and width, height and elasticity are no longer global variables, but belong to the container. To make the code work, you have to move the bounce() function from the Particle class to the Environment class and change width, height and elasticity to self.width, self.height and self.elasticity respectively.

I decided against moving the collide function into the Environment class as it didn't make much sense and theoretically we could create particles and test for collision without needing an environment container. However, there are a couple of line that refer to elasticity, so we need to change these:

elasticity = p1.elasticity * p2.elasticity
p1.speed *= elasticity
p2.speed *= elasticity

Then we need to give the Particle class an elasticity attribute, which I set to 0.9. This gives us more flexibility to have particles and walls with different elasticities.

Selecting particles

Finally, we need to add back the ability select and drag particles about. This is part of the controller in the MVC system, so shouldn't be part of the PyParticles module. Instead the following should go into the event handling code in the test program:

elif event.type == pygame.MOUSEBUTTONDOWN:
    selected_particle = env.findParticle(pygame.mouse.get_pos())
elif event.type == pygame.MOUSEBUTTONUP:
    selected_particle = None

Then add the findParticle() function into the Environment class remembering to change particles to self.particles.

Then, below the event handling code, add:

if selected_particle:
    selected_particle.mouseMove(pygame.mouse.get_pos())

To minimise the amount of code we need to write each time, I created a mouseMove() function to contain the code:

def mouseMove(self, (x, y)):
    dx = x - self.x
    dy = y - self.y
    self.angle = 0.5*math.pi + math.atan2(dy, dx)
    self.speed = math.hypot(dx, dy) * 0.1

Documentation

Since we've made a module that other people can use it's a good idea to add some comments to explain the classes and functions. If you've already add comments then that's great - I should have earlier, but was being lazy. Use triple-quoted comments at the beginning of functions to create a documentation string. This will be displayed by some IDEs and can be extracted by, for example:

print Particle.move.__doc__

You can find the final code attached at the bottom of this post and that includes documentation strings.

We've created a fairly flexible module that allows us to recreate our simulation with very little code (just 34 lines, including blanks ones). However, there is still one major inflexiblity in the module in terms of what behaviours the particles display. For example, at the moment, the code for making particles fall under gravity is commented out; if we uncomment it, then any simulation using PyParticles will be effected (and crash as gravity is no longer defined). In the next tutorial, we'll deal with this problem and generalise it for all particle behaviours.

AttachmentSize
particle_tutorial_10.txt1015 bytes
PyParticles.txt5.01 KB

Anonymous functions

In the previous tutorial, we created module, PyParticles, that allowed to us to recreate our particle simulation relatively easily. However, before we can use this module to create a range of different simulations, we need to be able to choose which behaviours the particles exhibit. In this tutorial, we will:

  • Introduce lambda for creating anonymous functions
  • Give the Environment class a dictionary of functions we can choose from

Like the previous tutorial, this one's not specific for Pygame. It introduces Python's anonymous functions, which are quite an advanced programming technique. It took a while to work out how and when to use them, but I think we now have a prefect example of how they can be useful. As usual, the final code is at the bottom of this post.

Particle functions

At present the Particle's move() function is responsible for changes in velocity due to gravity (or not if you commented it out) and drag. In the next tutorial, we'll model a cloud of gas in space, so we don't want to include drag or a gravitational force pulling all particles in the same direction (nor do we want the particle to bounce off the boundary of the screen). We therefore need to make these functions independent. To start with, move the code controlling gravity out of the move() function into a separate accelerate() function within the Particle class:

def accelerate(self, vector):
    (self.angle, self.speed) = addVectors((self.angle, self.speed), vector)

We don't necessarily need to move drag() into a separate function because by setting the mass_of_air variable to 0, drag becomes 1, so there's no reduction in speed. However, it's inefficient to multiply each particle's speed by 1 every tick of the simulation; it would be better if we could ignore drag unless specified otherwise. If we create a separate drag() function, then we could rewrite the Environment update() method like so:

def update(self):
    for particle in self.particles:
        particle.move()
        if self.mass_of_air != 0:
            particle.experienceDrag()
        if self.acceleration:
            particle.accelerate(self.acceleration)
        if self.hasBoundaries:
            self.bounce(particle)

This update() method tests whether to we should bother calculating drag. It also test whether to accelerate particles due to some constant force, such as gravity, which would be set as a vector belonging to the Environment class called self.acceleration. Finally it tests whether the particles should bounce off the walls or continue off into space. This is determined by a boolean (i.e. true or false), self.hasBoundaries.

The problem is that now we testing whether we should call these various functions for each of the particles, every tick of the simulation. If we have 100 particles, then that's 300 if statements every tick! We could improve the efficiency somewhat (to 3 if statements per tick) by putting each of the calls in a separate loop, with the test outside, like so:

for particle in self.particles:
    particle.move()

if self.mass_of_air != 0:
    for particle in self.particles:
        particle.experienceDrag()

if self.acceleration:
    for particle in self.particles:
        particle.accelerate(self.acceleration)

if self.hasBoundaries:
    for particle in self.particles:
        self.bounce(particle)

However, we still have to carry out each of the tests each tick of the simulation and the code is just inelegant. Ideally we should only have to carry out the test once at the start of the simulation.

Variable functions

My solution to this problem (and others probably exist), is to give the Environment object a list of the particle functions we want to use (move, drag, bounce etc.), and call each of them for each particle:

for particle in self.particles:
    for f in self.particle_functions:
        f(particle)

But how do we get our functions into a list? We can't just type:

self.particle_functions = [particle.experienceDrag, self.bounce(particle)]

We want to define the function list at the beginning of the simulation, before we've even defined particle. Even if we had defined the Environment object's list of particles, we need some way to refer each specific particle object. Furthermore, if we add self.bounce(particle)to the list, what we actual add is the result of calling self.bounce(particle), which is None because the function doesn't return anything.

The solution is to use the lamdba function to create anonymous functions. For convenience, I put the functions in a dictionary, so they can be referred to easily. To the Environment class's __init__() function, add:

self.particle_functions = []
self.function_dict = {
'move': lambda p: p.move(),
'drag': lambda p: p.experienceDrag(),
'bounce': lambda p: self.bounce(p),
'accelerate': lambda p: p.accelerate(self.acceleration)}

Lambda allows us to create a single line function without a name. Just as with normal functions they can take a parameter, which is given before the colon. For example, if we now type:

self.function_dict['move'](self.particles[0])

We will call lambda p: p.move() with the first particle in the Environment object's particle as the parameter. This will take the particle object (referred to as p in the lambda function) and call its move() function.

We can now give the Environment class a function to add specific functions to its particle_functions list.

def addFunctions(self, function_list):
    for f in function_list:
        if f in self.function_dict:
            self.particle_functions.append(self.function_dict[f])
        else:
            print "No such function: %s" % f

Now, when we start a new simulation, it's incredibly simple to define which functions to include:

env = PyParticles.Environment((width, height))
env.addFunctions(['move', 'accelerate', 'drag'])
env.acceleration = (math.pi, 0.002)

When the Environment's update() function is called, it will now call the move, accelerate and drag functions for each of its particles.

Two-particle functions

Sadly, the situation is not quite so simple, because we should also like to define which two-particle functions are called. For example, we might not want the particles to collide and in the next tutorial we'll want to add an attract() function which will take two particles as parameters. The following adaptation is required:

self.particle_functions1 = []
self.particle_functions2 = []
self.function_dict = {
'move': (1, lambda p: p.move()),
'drag': (1, lambda p: p.experienceDrag()),
'bounce': (1, lambda p: self.bounce(p)),
'accelerate': (1, lambda p: p.accelerate(self.acceleration)),
'collide': (2, lambda p1, p2: collide(p1, p2))}
        
def addFunctions(self, function_list):
    for func in function_list:
        (n, f) = self.function_dict.get(func, (-1, None))
        if n == 1:
            self.particle_functions1.append(f)
        elif n == 2:
            self.particle_functions2.append(f)
        else:
            print "No such function: %s" % func

Now the function dictionary values are a tuple that indicates whether the function takes 1 or 2 parameters, which we use to determine which of two function lists to add the function to. The Environment update() function should now be changed to:

def update(self):
    for i, particle in enumerate(self.particles):
        for f in self.particle_functions1:
            f(particle)
        for particle2 in self.particles[i+1:]:
            for f in self.particle_functions2:
                f(particle, particle2)

Now the collide() function can be included or ignored in the same way as for the other particle functions. The following tutorial will have an example of how our updated PyParticles module can be used.

Final note on efficiency

Note, our simulation still isn't as efficient as it could be since if we want to move the particles and make them experience drag, then we must make two function calls instead of one as before. Furthermore, we carry out the nested loop in the update() function regardless of whether we include any two-particle functions. This is a limitation of writing flexible code. If processing speed becomes an issue, then it's best to write a custom program. However, the PyParticles module is a good place to start prototyping a simulation and testing parameters; you can start weeding out further inefficiencies once you've got an interesting simulation running. PyParticles should be efficient enough to get you started.

These last two tutorial may have seemed like a lot of work to essentially get back where we started (and if you made it this far then I'm seriously impressed), but hopefully the following tutorials will demonstrate just how flexible our code now is. We can now use our moduel to create a range of apparently very different programs with relatively little effort. We'll start in the next tutorial, by making a simulation of a gas cloud collasping under its own gravity to form a solar system.

AttachmentSize
PyParticles2.txt5.89 KB

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:

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 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.

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 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:

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 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.

AttachmentSize
particle_tutorial_12.txt1.33 KB
PyParticles3.txt6.87 KB

Keyboard interactions

In the previous tutorial, we made a simulation of a gas cloud collapsing under its own gravity. One problem with the simulation is that, unless you are very lucky, the particles eventually move off the screen and you're left staring at a blank screen. In this tutorial, we will:

  • Create functions to scroll and zoom the display
  • Call these functions in response to keystrokes
  • Add the ability to pause our simulation
  • Link keystrokes to functions using a dictionary

These are highly useful for many simulations and other programs, such as games. Below is a video demonstrating scrolling and zooming during the simulation.

It's important to note that we don't want to actually change the position of any of the particles. This may seem odd, but instead of changing the position of the particles, we change how we display them. Scrolling and zooming are related to the display, and as before, we want to separate the display from the mechanics of the simulation. To illustrate why separating the display from the mechanics is important, imagine we zoom out of out simulation by dividing the x and y attributes of each particle by two. The particles would now be twice as close to one another as before, thus giving the appearance of zooming out. However, we would then have to take into account the difference in distances when checking for collisions and when calculating gravity. It's much simpler to keep the 'actual' position of the particles the same and just change how particles are displayed.

The UniverseScreen

Since we want to change how the universe is displayed, we need to create some variables to keep track of these changes. We could just create some global variables, but it's neater to create an object to tie all the variables and functions together. Create the following object:

class UniverseScreen:
    def __init__ (self, width, height):
        self.width = width
        self.height = height
        (self.dx, self.dy) = (0, 0)
        (self.mx, self.my) = (0, 0)
        self.magnification = 1.0
        
    def scroll(self, dx=0, dy=0):
        self.dx += dx * width / (self.magnification*10)
        self.dy += dy * height / (self.magnification*10)
        
    def zoom(self, zoom):
        self.magnification *= zoom
        self.mx = (1-self.magnification) * self.width/2
        self.my = (1-self.magnification) * self.height/2
        
    def reset(self):
        (self.dx, self.dy) = (0, 0)
        (self.mx, self.my) = (0, 0)
        self.magnification = 1.0

The UniverseScreen object contains the variables and functions required for scrolling, zooming, and resetting the display to its original state. Create a UniverseScreen instance with the same width and height as the universe instance and the screen instance:

universe_screen = UniverseScreen(width, height)

Scrolling

Scrolling involves adding a constant value to either the x or y position of the particles. For example, if we scroll left, then we display every object a few pixels to the right, which means adding a small amount to their x attribute. We store the amount that we change the x and y attributes as the UniverseScreen's dx and dy attributes. We can now alter the display code to:

We therefore update the code for displaying particles to:

x = int(p.x + universe_screen.dx)
y = int(p.y + universe_screen.dy)

if size < 2:
    pygame.draw.rect(screen, p.colour, (x, y, 2, 2))
else:
    pygame.draw.circle(screen, p.colour, (x, y), size, 0)

We can now scroll left, for example, by calling:

universe_screen.scroll(dx=1)

This will move all the particles 40 pixels right (the width of the screen divided by 10). Later we'll attach this function call to a keyboard input.

Zooming

The equation for zooming is more complicated because not only do we have to change the relative position of the particles (the more we zoom in, the further apart they should appear), but we also need to alter the x and y attributes such that a particle in the dead centre of the screen remains in the centre. We could do this by altering the dx and dy attribute of the UniverseScreen, but I found it easier to create new variables, mx and my. 

If you look at the zoom() function above, you can see the first thing it does is to change the magnification attribute, which determines how zoomed in or out we are. A magnification of two means that we are twice as "zoomed in" so the distance between particles should appear twices as big (and so particles should also appear to move twice as fast). In addition, we should double the size of the particles and half the effect of scrolling so that we always move a tenth of the screen. Below is the modified display code:

x = int(universe_screen.mx + (universe_screen.dx + p.x) * universe_screen.magnification)
y = int(universe_screen.my + (universe_screen.dy + p.y) * universe_screen.magnification)
size = int(p.size * universe_screen.magnification)
        
if size < 2:
    pygame.draw.rect(screen, p.colour, (x, y, 2, 2))
else:
    pygame.draw.circle(screen, p.colour, (x, y), size, 0)

Notice that as well as multiplying the particles' x and y attributes by the magnification parameter we also multiply the dx and dy parameters. This is because previous scrolling movements also need to be magnified. If you look back at the scroll() function, you can see that we divide the distance we scroll by the magnification value, so as we zoom in, we still scroll a tenth of a screen at a time.

The two parameters we haven't yet discussed are mx and my. These values are set by the zoom() function and are required so that as you zoom in and out, the screen remains centred.

Responding to keyboard events

We now have functions that change the attributes of universe_screen and a display that uses these attributes to correctly draw the particles. The next set is to call these functions in response to keystrokes. Responding to keystroke events is very similar to responding to mouse events, only instead of responding to MOUSEBUTTONDOWN or MOUSEBUTTONUP events, we respond to KEYDOWN events. For example:

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

One difference is that KEYDOWN events have an attribute, event.key. If you run the program now, you will see a number printed to your terminal whenever you press a key. Each key is associated with a different number, so we can determine which key has been pressed. Rather than work out which number is associated with each key (e.g. pressing 'a' gives us 97), Pygame defines constants with fairly easy-to-remember names that represent each of keys (e.g. pygame.K_a = 97). You can find a full list here. For example, we respond to the arrow keys like so:

if event.type == pygame.KEYDOWN:
    if event.key == pygame.K_LEFT:
        universe_screen.scroll(dx=1)
    elif event.key == pygame.K_RIGHT:
        universe_screen.scroll(dx=-1)
    elif event.key == pygame.K_UP:
        universe_screen.scroll(dy=1)
    elif event.key == pygame.K_DOWN:
        universe_screen.scroll(dy=-1)
    elif event.key == pygame.K_EQUALS:
        universe_screen.zoom(2)
    elif event.key == pygame.K_MINUS:
        universe_screen.zoom(0.5)
    elif event.key == pygame.K_r:
        universe_screen.reset()

We can now call scroll events when the arrow keys are pressed, zoom out when the minus key is pressed, zoom in when the equals key is pressed and call a reset() function when 'r' key is pressed. Note that because of the way the scroll() function is written, it can take two parameters (and can scroll diagonally), even though we don't use it this way. The reset() function is a method of UniverseScreen that resets its attributes:

def reset(self):
    (self.dx, self.dy) = (0, 0)
    (self.mx, self.my) = (0, 0)
    self.magnification = 1.0

Play/Pause

In order to pause the simulation we need to add a conditional to the code that updates the universe:

if not paused:
    universe.update()

Set paused to equal False somewhere outside the simulation loop, then add a new event.key test:

elif event.key == pygame.K_SPACE:
    paused = (True, False)[paused]

This toggles the 'paused' value, so if it's True, then it become False; if it's False then it becomes True. The reason this code works is that True is considered equal to 1, while False is considered equal to 0 (more details here), so we can use the 'paused' variable to select from a tuple containing True and False. If 'paused' is False, then the first value is selected; if 'paused' is True, then the second value is selected. Don't forget to set 'paused' to False (or True if you prefer) somewhere before the simulation loop starts.

Bonus section - more anonymous functions

The simulation is now complete. However, the code is a little inelegant and inefficient, what with having six consecutive elif statements. In some programming languages there is a case ... switch pattern for this sort of situation. Python doesn't have such statements, but we can get around this by using a dictionary of functions. If you've had enough of anonymous functions after the previous tutorial, then feel free to skip this section. Maybe it's just because I've been learning Lisp that I have a compulsion to include the lambda function anywhere I can.

To avoid the elif statements, we create a dictionary that links the keyboard input with the relevant function:

key_to_function = {
    pygame.K_LEFT:   (lambda x: x.scroll(dx = 1)),
    pygame.K_RIGHT:  (lambda x: x.scroll(dx = -1)),
    pygame.K_DOWN:   (lambda x: x.scroll(dy = -1)),
    pygame.K_UP:     (lambda x: x.scroll(dy = 1)),
    pygame.K_EQUALS: (lambda x: x.zoom(2)),
    pygame.K_MINUS:  (lambda x: x.zoom(0.5)),
    pygame.K_r:      (lambda x: x.reset())}

You should be able see what's going on here: the keystroke value is linked to a lambda function, which takes a parameter, x, and calls a function of x. We are going to call each function with the parameter universe_screen, so we will call universe_screen's functions. In fact, this is the main reason for creating a UniverseScreen class. We can now use this dictionary to immediately call the right function (once we've checked that the keystroke is in the dictionary):

if event.key in key_to_function:
    key_to_function[event.key](universe_screen)
elif event.key == pygame.K_SPACE:
    paused = (True, False)[paused]

Note that we can't create a lambda function that changes the value of 'paused'. Lambda functions cannot have 'side-effect', i.e. they can't change values defined outside of their scope. We can get around this by making 'paused' an attribute of the UniverseScreen class (because we pass universe_screen to the lambda function as x, so could access it as x.paused). I've left it this way to highlight a limitation of anonymous functions.

AttachmentSize
particle_tutorial_13.txt2.96 KB
PyParticles3.txt6.87 KB

Specific attraction: springs

A couple of tutorial, we created a simulation in which every particle attracted every other particle. Another classic system that physicists have long considered is one in which bodies are are moved towards or away from one another by a spring that connects them. In this tutorial, we will:

  • Create a Spring object
  • Link particles with springs

Once we have created springs, we can use them to connect particles and simulate complex soft bodies. A simple example is shown in the video below. The springs can also represent edges in a graph, with the particles representing vertices.

Setting up the simulation

We start in almost exactly the same way we started the gas cloud simulation, with some bolierplate code to create our screen. The only differences are that this time we will need to import pi from the math module and we can inlcude the pause control.

from math import pi
import random
import pygame
import PyParticles

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

paused = False
running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
        elif event.type == pygame.KEYDOWN:
            if event.key == pygame.K_SPACE:
                paused = (True, False)[paused]

    pygame.display.flip()

As before, we then create an Environment, only this time, the background colour is white and we want our particles, to move, bounce off the wall, collide with one another, experience drag and experience downwards acceleration (gravity). I've also changed the mass_of_air to 0.02, as that seems to give nicer results.

universe = PyParticles.Environment((width, height))
universe.colour = (255,255,255)
universe.addFunctions(['move', 'bounce', 'collide', 'drag', 'accelerate'])
universe.acceleration = (pi, 0.01)
universe.mass_of_air = 0.02

Adding particles

Next we add some particles; the exact parameters aren't too important.

for p in range(4):
    universe.addParticles(mass=100, size=16, speed=2, elasticity=1, colour=(20,40,200))

Then, in the main loop, we update and display the particles.

if not paused:
    universe.update()
        
screen.fill(universe.colour)
    
for p in universe.particles:
    pygame.draw.circle(screen, p.colour, (int(p.x), int(p.y)), p.size, 0)

And we can add some code to allow us to select and move particles. Outside the main loop we add:

selected_particle = None

And inside the main loop, where we deal with user interactions, we add:

    elif event.type == pygame.MOUSEBUTTONDOWN:
        selected_particle = universe.findParticle(pygame.mouse.get_pos())
    elif event.type == pygame.MOUSEBUTTONUP:
        selected_particle = None

if selected_particle:
    selected_particle.mouseMove(pygame.mouse.get_pos())

If you run the code now, you should see some particles bouncing around and coming to halt. You can pick them and a throw them about. I hope you agree that with our module it's relatively simple to set up such a simulation. Now to add something new.

The Spring object

The force exterted by a spring is given by Hooke's law, which tells us that F = -kx, or that the force is equal to a constant (sometimes called the spring constant), k, multiplied by its displacement from its equilibrium position. Our spring object will therefore have a length attribute, representing the length it 'tries' to obtain, and a strength attribute, representing how quickly it tries to reach that length (i.e. the spring constant). It will also record which particles are at either end.

To the PyParticles file add a new class somewhere:

class Spring:
    def __init__(self, p1, p2, length=50, strength=0.5):
        self.p1 = p1
        self.p2 = p2
        self.length = length
        self.strength = strength

Then to the Environment class we add a list to keep track of spring objects and a method to add spring object to it.

self.springs = []
def addSpring(self, p1, p2, length=50, strength=0.5):
    """ Add a spring between particles p1 and p2 """
    self.springs.append(Spring(self.particles[p1], self.particles[p2], length, strength))

Now we can add springs in our simulation code like so:

universe.addSpring(0,1, length=100, strength=0.5)
universe.addSpring(1,2, length=100, strength=0.1)
universe.addSpring(2,0, length=80, strength=0.05)

These lines of code add a spring between particles 0 and 1, particles 1 and 2, and particles 2 and 0, thus making a triangle. We can represent the springs as lines between the particles by adding the following after the code for drawing particles.

for s in universe.springs:
    pygame.draw.aaline(screen, (0,0,0), (int(s.p1.x), int(s.p1.y)), (int(s.p2.x), int(s.p2.y)))

If you run the simulation now, you should see the particles bouncing around as usual, but with three of them joined up in a triangle.

To make our springs actually exert and effect, we need to add a method that accelerates particles at either end in line with the spring:

def update(self):
    dx = self.p1.x - self.p2.x
    dy = self.p1.y - self.p2.y
    dist = math.hypot(dx, dy)
    theta = math.atan2(dy, dx)
    force = (self.length - dist) * self.strength
        
   self.p1.accelerate((theta + 0.5 * math.pi, force/self.p1.mass))
   self.p2.accelerate((theta - 0.5 * math.pi, force/self.p2.mass))

This function is almost identical to the attract() function that we wrote for gravitational attraction, which is unsurprising. Both measure the distance and angle between two particles and calculate the force acting on the particles using this information. The only difference is the equation for force.

Finally we need to call the springs' update function within the Enivironment's update function:

for spring in self.springs:
     spring.update()

Day dreaming

We now have the building blocks for all kinds of simulations and I recommend playing about with different parameters and building structure. By linking particles we can make complex structures. We could prehaps make an Angry Birds style game if we allow springs to break under a certain force (though Angry Birds actually uses rigid body physics). We could allow colliding particles to make and break connections, thus simulating a virtual chemistry. We could also create some sort of creature in which the springs act as muscles, which can be contracted and relax, so the creature can move. At some point I hope to explore some of these ideas.

As an quick example, below shows the approximation to two solid object: a sphere (actually a dodecagon) and a square, defined by their edges and some diagonals. You can see there are a few problems, such as the circle can slip into the square, and the circle can collapse into a tightly-bound mess.

AttachmentSize
particle_tutorial_14.txt1.61 KB
PyParticles4.txt7.8 KB