Matlab Particles 2.0

Contents

1 Genesis
1.1 Particle System Creation
1.2 Particle Creation
1.3 Spring Creation
1.4 Attraction Creation
1.5 Simulation

2 Demos
2.1 Demo 1: Free Fall
2.2 Demo 2: Bullet Time
2.3 Demo 3: Bungee Jumping
2.4 Demo 4: Gimme ya Energy!
2.5 Demo 5: Magic Chain Cable
2.6 Demo 6: Heavy Chain Mail
2.7 Demo 7: Don’t Touch me!
2.8 Demo 8: Keep it up
2.9 Demo 9: Catch me if you Can (There’s a Hole in the Bucket) .
2.10 Demo 10: The Hose
2.11 Demo 11: Polyhedrons
2.11.1 Number of Particles: 3
2.11.2 Number of Particles: 4
2.11.3 Number of Particles: 5
2.11.4 Number of Particles: 6
2.11.5 Number of Particles: 8
2.11.6 Number of Particles: 12
2.11.7 Number of Particles: 20
2.12 Demo 12: Three-Body Eight . .

3 Mathematical Background
3.1 May the Force be with you
3.1.1 Gravity
3.1.2 Drag
3.1.3 Inertial Force
3.1.4 Attraction Force
3.1.5 Spring Force
3.1.6 Damping Force
3.2 Differential Equations

4 Particle System Object
4.1 Class Definition particle system
4.2 Properties particle system
4.3 Constructor particle system
4.4 Method get particles positions
4.5 Method get particles velocities
4.6 Methods kill spring, kill attraction
4.7 Method kill particle
4.8 Method advance time

5 Particle System Object (Private Methods)
5.1 Private Method kill old particles
5.2 Private Method get phase space state
5.3 Private Method compute state derivative
5.4 Private Method set phase space state
5.5 Private Method aggregate forces
5.6 Private Method clear particle forces
5.7 Private Method aggregate springs forces
5.8 Private Method aggregate attractions forces
5.9 Private Method aggregate drag forces
5.10 Private Method aggregate gravity forces
5.11 Private Method get particles accelerations
5.12 Private Method advance particles ages
5.13 Private Method update graphics positions

6 Particle Object
6.1 Class Definition and Properties particle
6.2 Constructor particle
6.3 Private Method append
6.4 Methods add force, clear force
6.5 Method delete
6.6 Method set.fixed
6.7 Method set.position
6.8 Method update graphics position

7 Spring Object
7.1 Class Definition and Properties spring
7.2 Constructor spring
7.3 Private Method append
7.4 Method delete
7.5 Method update graphics position

8 Attraction Object
A Vector Projection
B Particle System Dependency
C Particle, Attraction, and Spring Dependency

When Karl Sims did his award-winning computer animation ”Particle Dreams” twenty(!) years ago, he tortured a Connection Machine CM-2 computer with as many as 65,536 processors, using one processor for the simulation of each particle.

illustration not visible in this excerpt

Figure 1: 1988 Computer Animation ”Particle Dreams”9

Today we simulate tens of thousands of particles in real-time on a single cpu (Figure 2)

- even in a browser plugin (Figure 3) - and advanced particle systems have become common practice for the simulation of snow, rain, dust, smoke, fire, and explosions in most computer games. Modern simulation environments like Processing3can be used to produce such astonishingly addicting games as Falling Sand Game6, sodaplay10, BallDroppings7, and Souptoys13.

In 2006, Traer Bernstein2 wrote a pretty impressing particle physics library for Process- ing, which actually was the inspiration for this particle system toolbox in Matlab. As a matter of fact, object oriented programming in Matlab is not really the fastest lane on the particle system highway; we are back at the good old days of some ten or twenty real-time particles. But - the main purpose of this toolbox has never been to develop state-of-the-art computer games; it was rather planned as an educational, interactive learning-by-doing playground, with the aim to understand the mechanical interactions (and maybe the mathematical background) of the particle system components. Have fun!

illustration not visible in this excerpt

Figure 2: Particle System API 5

illustration not visible in this excerpt

Figure 3: Flash 9 Particle System 1

illustration not visible in this excerpt

Figure 4: Souptoys 13

illustration not visible in this excerpt

Figure 5: Falling Sand Game 6

illustration not visible in this excerpt

Figure 6: Sodaplay 10

illustration not visible in this excerpt

Figure 7: BallDroppings 7

1 Genesis

The main purpose of this section is to give you a brief overview on how particles, springs, and attractions are created and implemented into the particle system.

1.1 Particle System Creation

Every particle system simulation begins with the creation of a particle system object¹ that maintains and manages all particles, springs, and attractions

>> Particle_System = particle_system

illustration not visible in this excerpt

The command particle_system opens an empty 3-D axes in a maximized window and creates an empty particle system with default properties

- Acceleration due to gravity:[0 0 0]
- Aerodynamic/fluid resistance (drag): 0
- Axes limits: ±1

After the creation of the particle system object its properties can be defined

>> gravity = [0 0 -9.81];

and set

>> Particle_System.gravity = gravity

illustration not visible in this excerpt

In many cases it might be more convenient to set all three particle system properties (gravity, drag, and limits) directly during the creation of the object

>> Particle_System = particle_system ([0 0 -9.81], 0, 1)

illustration not visible in this excerpt

1.2 Particle Creation

A particle object with default properties is created and incorporated into the parti- cle system by the particle command using the particle system (handle) as its only parameter

>> Particle = particle (Particle_System)

illustration not visible in this excerpt

Again, the particle properties (mass, initial position and velocity, fixed/free, and lifespan) can either be set individually with an existing particle object

>> Particle.mass = 42

illustration not visible in this excerpt

:

or at once during the creation of the object

>> Particle = particle ...

(Particle_System, 42, [0 0 0], [0 0 0], false, inf)

illustration not visible in this excerpt

:

The particle is automatically incorporated into the particle system and a red dot (symbolizing the particle) is drawn in the axes of the particle system. After its creation, the color (or any other property) of the graphical particle representation can be modified by means of its graphics handle set (Particle.graphics_handle, ’color’, ’blue’);

1.3 Spring Creation

A default spring between two already existing particles Particle_1 and Particle_2 is created and incorporated into the particle system via

>> Spring = spring (Particle_System, Particle_1, Particle_2)

particle_1: [1x1 particle]

particle_2: [1x1 particle]

rest: 1

strength: 1

damping: 0

graphics_handle: 195.0077

Again, the long version directly sets the rest length and the strength of the spring and the damping coefficient of an additional damper parallel to the spring

>> Spring = spring ...

(Particle_System, Particle_1, Particle_2, 2, 1, 0)

illustration not visible in this excerpt

The default graphical representation of a spring is a solid blue line between the corresponding particles.

1.4 Attraction Creation

Just like a spring, an attraction (force) connects two particles

>> Attraction = attraction ...

(Particle_System, Particle_1, Particle_2)

illustration not visible in this excerpt

Its properties are the gravitational attraction strength (which can be made negative for a repulsion) and the minimum distance, below which the attraction force will not grow any further. Analogous to the particle and the spring, the following attraction uses are valid syntax

>> Attraction = attraction ...

(Particle_System, Particle_1, Particle_2, 1, 0);

>> Attraction.minimum_distance = 0.1;

>> strength = Attraction.strength;

>> set (Attraction.graphics_handle, ’visible’, ’off’);

An attraction is graphically represented by a dotted blue line.

1.5 Simulation

After at least one particle has been thrown into life, the simulation can be started

step_time = 0.01;

for i = 1 : inf

Particle_System.advance_time (step_time);

end

Since inf has been used as the upper bound of the for-loop, the simulation can only be stopped by Ctrl-C or by closing the corresponding figure window. In the loop, every call to advance_time integrates the underlying nonlinear differential equation system one single time step further, the length of which is defined in the variable step_time. Decreasing step_time makes the simulation smoother and more exact, but also slower; increase step_time to make the particles move faster. If you get too greedy speed-wise, the simulation becomes bumpier and finally unstable. In the simulation loop the current position of the mouse pointer can be determined

current_point = mean (get (gca, ’currentpoint’));

and e. g. used to move² certain particles ”by hand” Particle.position = current_point;

2 Demos

It might be a good idea to read through the demos in chronological order, since the grade of general parameter explanation detailedness is higher in the first demo descriptions.

2.1 Demo 1: Free Fall

The first - pretty boring - demo simulates the free fall of a single mass in a gravity field. After the creation of a default particle system (no drag, unit limits) with ”earthy” gravity

Particle_System = particle_system;

gravity = [0 0 -9.81];

Particle_System.gravity = gravity

the default 3-D view is reduced to two dimensions

view (0, 0);

A single particle with default properties (unit mass, initial position at the origin, no initial velocity, not fixed, infinite lifespan) is created and incorporated into the particle system

Particle = particle (Particle_System);

A step time of 1 millisecond seems to produce a smooth and visually traceable particle motion on a 3 GHz computer³:

step_time = 0.001;

Finally, the most simple form of a simulation loop can be used

for i = 1 : 451

Particle_System.advance_time (step_time);

end

Choosing 451 as the number of steps makes the simulation stop when the particle reaches the lower axes limit⁴.

illustration not visible in this excerpt

Figure 8: Free fall

As you might have expected, the unspectacular simulation (Figure 8) shows a red dot falling down with (linearly) increasing velocity.

2.2 Demo 2: Bullet Time

The second demo simulates the motion of a bullet in zero gravity with aerodynamic resistance⁵. The few lines of source code are very similar to the first demo. Create a particle system with no gravity but a drag coefficient of 10 and set the view to 2D

Particle_System = particle_system ([0, 0, 0], 10, 1);

view (0, 0);

Position a particle on the left edge (-1) of the axes and give it an initial velocity of 20

to the right

Particle = particle ...

(Particle_System, 1, [-1, 0, 0], [20, 0, 0], false, inf);

Simulate infinitely long with a step time of 1 millisecond

step_time = 0.001;

for i = 1 : inf

Particle_System.advance_time (step_time);

end

The simulation depicted in Figure 9 produces the same low level of enthusiasm as the previous demo.

illustration not visible in this excerpt

Figure 9: Bullet time

The ”bullet” starts with initial kinetic energy at the left hand side of the axes (Figure 9 a), becomes slower and slower due to the aerodynamic resistance (Figure 9 b) until it comes to a complete⁶ stop at the right hand side, when all of its energy has been dissipated (Figure 9 c).

2.3 Demo 3: Bungee Jumping

Wikipedia defines bungee jumping as

an activity in which a person jumps off from a high place (generally of several hundred feet/meters) with one end of an elastic cord attached to his/her body or ankles and the other end tied to the jumping-off point.

Well, what do you need for a half-decent bungee jumping simulation? The jumping person is just a free particle with a certain mass and the elastic cord can be modeled by a spring between the jumper and another (hopefully fixed) ”particle”. Add some aerodynamic resistance and a bit of spring damping and the jump will look quite realistic. In order to open up the second dimension, the jumping-off point should not equal the fixed chord end point⁷.

< Genesis mode on >

In the beginning The User created the particle system and the particles. And The User said, Let there be gravity and aerodynamic resistance, and there was gravity and aerodynamic resistance

Particle_System = particle_system ([0, 0, -9.81], 5, 200);

The User saw all that he had made, and it was very good (even in 2D)

view (0, 0);

< Scientific mode on >

The first particle acts as the anchor point the chord is ”tied to”. Its initial position is fixed (true) somewhere in the upper half of the axes ([0, 0, 150]) and it has an infinite life span⁸. Since the particle is fixed and does not move, its mass is completely irrelevant⁹. The default graphical representation of a fixed particle is a red asterisk

Particle_1 = particle ...

(Particle_System, 1, [0, 0, 150], [0, 0, 0], true, inf);

The second particle represents the jumper. It should have a reasonable mass¹⁰ and an initial position at the height of the first particle, with a horizontal offset equalling the

rest length of the spring

Particle_2 = particle ...

(Particle_System, 70, [100, 0, 150], [0, 0, 0], false, inf);

The two particles are now used to define both ends of the spring to be created

Spring = spring ...

(Particle_System, Particle_1, Particle_2, 100, 6, 0.1);

The spring has a rest length¹¹ of 100, a strength¹² of 6, and a damping factor¹³ of

0.1. To make the graphical representation (which is a solid blue line by default) look a bit more like a real spring, you can increase its line width and change its line style to dotted

set (Spring.graphics_handle, ’linewidth’, 10, ’linestyle’, ’:’)

The particle system is now ready to be pushed into life

step_time = 0.1;

for i = 1 : inf

Particle_System.advance_time(step_time);

end

illustration not visible in this excerpt

Figure 10: Bungee jumping

Departing from Figure 10 a with a relaxed spring, the free particle is accelerated down- wards in the gravity field and immediately to the left by the horizontal component of the spring force. At about 12 sec, it reaches its lowest point (Figure 10 b) and is driven back up by the extended spring. After a few more vertical and horizontal oscillations the atmospherical and the spring damping have eaten up all kinetic energy. Figure 10 c shows the ”final” steady state, where the weight of the particle and the spring force are in equilibrium.

2.4 Demo 4: Gimme ya Energy!

This demo does not offer too many new insights, compared to the previous one; it is just a nice two-masses-two-springs experiment, where a big particle sucks energy from a smaller one (just like in real life. . . ). There is no gravity nor drag in this environment

Particle_System = particle_system ([0, 0, 0], 0, 3);

The first particle is just a nail in the middle of a wall, where you can later on fix the first spring

Particle_1 = particle ...

(Particle_System, 1, [0, 0, 0], [0, 0, 0], true, inf);

The second particle is free and has a bigger mass (10)

Particle_2 = particle ...

(Particle_System, 10, [1, 0, 0], [0, 0, 0], false, inf);

which can also be graphically represented

set (Particle_2.graphics_handle, ’markersize’, 60);

The third particle is smaller (1) than the second one, but has an initial velocity (5)

driving it upwards

Particle_3 = particle ...

(Particle_System, 1, [2, 0, 0], [0, 0, 5], false, inf);

Finally we create two identical springs with a strength of 10 and a damping factor of 1;

one spring between the (fixed) first particle and the second particle

Spring_1 = spring ...

(Particle_System, Particle_1, Particle_2, 1, 10, 1);

and another spring between the second and the third particle

Spring_2 = spring ...

(Particle_System, Particle_2, Particle_3, 1, 10, 1);

Ready to rumble

step_time = 0.1;

for i = 1 : inf

Particle_System.advance_time (step_time);

end

The snapshot in Figure 11 a shows the initial positions of both particles and springs. The initial velocity of the small particle that will move it upwards in the next simulation steps is not directly visible. On its way up, the small particle exerts a force through the spring on the big particle; but due to the inertia of the big particle, at first, the small particle rotates about the big one (Figure 11 b). During its rotations, the small particle transfers more and more energy to the big particle, until finally the latter is rotating so fast about the fixed particle that the small particle is not able to do full rotations about the big one any more (Figure 11 c).

illustration not visible in this excerpt

Figure 11: Gimme ya energy!

2.5 Demo 5: Magic Chain Cable

If you take a bunch of particles, arrange them in one long row and connect them with springs, you end up with something like a worm, a rope, a rubber band, a pearl necklace, or a chain cable - depending on the number of particles and the parameters you choose for mass, spring strength, damping, . . . Define some gravity (-1), a bit of drag (0.1) and an ample axes (10)

Particle_System = particle_system ([0, 0, -1], 0.1, 10);

The number of particles (and the corresponding number of springs) directly determine the speed of the simulation. A single core 3 GHz machine running Matlab object oriented programming would be hopelessly overextended even with a hundred particles. Therefore

n_particles = 20;

Next, the 20 particles are created and incorporated into the particle system

for i = 1 : n_particles

Particles{i} = particle ...

(Particle_System, 0.2, [0, 0,0], [0, 0, 0], false, inf);

if i == 1 || i == n_particles Particles{i}.fixed = true;

end

end

Additionally, the first and the last particle of the chain are declared as ”fixed”. This is useful, because their motion should not be influenced by the other particles, but will be externally modified in the simulation loop. The particles are connected by pretty stiff (100) springs with unit length (1) and a reasonable damping factor (4)

for i = 1 : n_particles - 1

Spring{i} = spring ...

(Particle_System, Particles{i}, Particles{i + 1}, 1, 100, 4);

end

Let the games begin

step_time = 0.05;

for i = 1 : inf

Particles{1}.position = ...

[10*sin(0.01*i), 10*sin(0.011*i), 10*sin(0.012*i)];

Particles{n_particles}.position = ...

[-10*sin(0.013*i), 10*sin(0.014*i), -10*sin(0.015*i)];

Particle_System.advance_time (step_time);

end

In the simulation loop the first (Particles{1}) and last (Particles{n_particles}) particles are permanently moved by altering their position property. Both end particles are driven by slow, slightly distorted, harmonic oscillators with an amplitude of 10 in all three spacial dimensions. As a result, the particles seem to float randomly through the axes cube.

Since all particles have their initial positions at the origin, all you can see in Figure 12 a is one red dot. The springs are not visible at all; they are all compressed to zero length and therefore exert expanding forces on the particles immediately. In Figure 12 b the sine generators have moved the end particles apart and the springs have reached their rest length. Since the distance between the end particles is much smaller than the number of (unit) springs, the chain is not tense but relaxed in a random pattern. After the end particles have been forced into different vertices of the axes cube (Figure 12 c), the length of the chain is greater than the sum of the spring rest lengths, straightening the chain and extending every single spring.

illustration not visible in this excerpt

Figure 12: Magic chain cable

2.6 Demo 6: Heavy Chain Mail

This demo is an attempt¹⁴ to simulate the motion of a medieval knight’s heavy chain mail. With adapted particle and spring parameters, other types of cloth can be simulated as well. To make things look realistic, there has to be some gravity (-0.4) and a little bit of drag (0.1)

Particle_System = particle_system ([0, 0, -0.4], 0.1, 5);

view (0, 0);

For a more detailed view the axes limits have to be readjusted

axis ([0, 6, -1, 1, 0, 6]);

Again, we are talking about Matlab object oriented programming; so let’s not over- strain the CPU and restrict ourselves to 25 particles, arranged in a square matrix-like

structure of 5 rows by 5 columns

for col = 1 : 5

for row = 1 : 5

Particles{row, col} = particle ...

(Particle_System, 1, [col, 0, row], if row == 5 && col == 1 || row == 5

Particles{row, col}.fixed = true;

end

end

end

After their creation, the upper left (5th row, 1st column) and upper right (5th row, 5th column) particles are nailed (fixed) to the wall. Next, you have to create the two-dimensional mesh by connecting every particle to its nearest neighbors with unit length springs. This can be done in two similar steps; one double loop for the horizontal springs

for col = 1 : 4

for row = 1 : 5

Springs_horizontal{row, col} = spring (Particle_System, ...

Particles{row, col}, Particles{row, col + 1}, 1, 100, 1);

end

end

and another one for the vertical springs

for col = 1 : 5

for row = 1 : 4

Springs_vertical{row, col} = spring Particles{row, col}, Particles{row

end

end

The simulation loop looks pretty unspectacular

step_time = 0.1; for i = 1 : inf

(Particle_System, ... + 1, col}, 1, 100, 1);

Particle_System.advance_time (step_time);

end

[...]

---

¹ As a convention throughout this toolbox, Matlab objects names begin with uppercase letters.

² Actually, the command does not directly move the graphical particle representation, but only sets the position property of the particle object. The corresponding graphics object is then automatically updated during the next simulation cycle.

³ If you use a faster machine, you might want to reduce the step time (and v. v.).

⁴ The junior high proof is up to the reader . . .

⁵ Can you imagine an environment (or a movie title) where this experiment could take place?

⁶ Actually, the bullet will not really stop before Judgment Day, but its velocity will rapidly become too small to be visually detected.

⁷ Besides, this overcomes the problem that a spring (but not a bungee chord) produces a pressing force if its length is less than its rest length. A similar real-world effect that cannot easily be modeled by a simple spring is the fact that the first part of the jump is actually a free fall where the chord does not produce any force at all.

⁸ If you feel a strong urge to act out your sadistic touch, you could reduce the life span of the first particle to e. g. 42 and see what happens to the jumper. . .

⁹ Nevertheless, in order to avoid ”Divide by zero” warnings, the mass of a particle should never be zero.

¹⁰ Yes, 70 is a reasonable mass for an adult in the SI metric system.

¹¹ Rest length: Length of a spring producing no force

¹² Strength: Deflection-dependent force coefficient

¹³ Damping: Velocity-dependent force coefficient

¹⁴ Traer Bernstein[2] can do that much better. . .