28 Oct 2010 Code on Github


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 can be found by clicking the Github link above.


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 let's 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.


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


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:

$\large{\frac{1}{\text{mass of particle}\ +\ \text{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:

$\text{particle.speed}\ =\ \text{particle.speed}\ -\ \text{particle.speed} \left( \frac{1}{\text{particle.mass}\ +\ \text{mass_of_air}}\right)$

Which can be rewritten:

$\text{particle.speed}\ \text{*=}\ \frac{\text{particle.mass}}{\text{particle.mass}\ +\ \text{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:

$\text{particle.speed}\ \text{*=}\ \left(\frac{\text{particle.mass}}{\text{particle.mass}\ +\ \text{mass_of_air}}\right)^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 need to define a constant called mass_of_air, which represents the mass (or momentum) of 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

To actually implement the next system of drag, change the line in the move() function from:

self.speed *= drag


self.speed *= self.drag

Before running the simulation, I think it's 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:

Notice how the denser (darker balls) are deflected less when they collide with the lighter balls and can cause the lighter balls to move a lot faster.

If you change the mass_of_air value to 2.0, the simulation will appear to be running in syrup; if you change it to 0, there will be no air resistance.

Comments (6)

Romain on 3 Jul 2013, 7:51 a.m.


First of all, thanks for the tutorial, it's very useful.

Now, I'm a bit puzzled about the drag. You gave this equation:
particle.speed = particle.speed - particle.speed * (1 / (particle.mass + mass_of_air))

Which can be simplified to:

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

But you wrote that it was equivalent to:

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

Did I miss something?

Anonymous on 3 Jul 2013, 7:19 p.m.

I get an indent error everytime I made a little change on your code. It's at the line where the overlap is, BUT I cannot understand why do I get it.

Could you help me?

Peter on 7 Jul 2013, 7:52 p.m.

Hi Romain,

You're right, that simplification doesn't work. I can't remember what I was thinking. I think I started with the final equation and was trying to justify it. I'll see if I can make some sense of it. But for now, I think using the final equation works reasonably well.

@ Anonymous - it might be that you're using tabs and when you copy and paste, the code uses spaces. Your text editor should be able to convert between tab and spaces, so make sure everything is using the same whitespace.

Ian on 22 Jul 2014, 4:54 a.m.

I've been using this as a guide to write my own gravity simulator in python v3.3 for a game I'm writing and so far it has gone quite well, but the new for collide function needs to be broken down. I haven't been able to figure out what it doing exactly because it isn't explained and all the math is in the function imput making it hard to read. If each argument could be assigned to variables that alone would really help.

Anonymous on 2 Mar 2016, 10:08 a.m.

I am getting an error with the looping of particles on line 137. Thz error i am getting is, Traceback (most recent call last):

File "C:\Python33\", line 137, in <module>

for particle2 in my_particles[i+1]:

TypeError: 'Particle' object is not iterable

Anonymous123 on 14 Jun 2016, 9:28 a.m.


thank you for the very useful tutorial.

I have a question for a actual project: How can I combine this with an accelerometer sense? Instead of the mouse interactions I would use the accelerometer data.

How does this work? Thanks for your help.