Thursday, 28th October 2010
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.
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.
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
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
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
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.