I have heard that a trip to a sleep doctor entails hooking up a bunch of sensors to you and they tell you things like how many times you woke up in the night. I don’t have a bunch of sweet EEG nodes, but I realized that I could easily create something that measures my relative movement throughout the night. Looking at an over-arching chart would probably give me a good idea of how I was really sleeping.
The Parts
I am currently working on a system which involves control theory and I purchased one of these IMUs from Sparkfun.
Sparkfun Razor 6DOF IMU
This is a 6 degree of freedom inertial measurement unit. It uses 2 gyroscopes and an accelerometer to give you information about your hardware’s orientation. I am waiting for some parts to come in so I just had it hooked up to one of my development boards just hanging out. I realized I could just use the accelerometer chip to relatively measure “how much” I was moving.
Accelerometers are really simple devices and they have become completely ubiquitous. If you are reading this on a mac or an iphone, then you are using one right now. A 3 axis accelerometer, such as the one on this board, can measure the orientation of the chip in relation to the earth’s magnetic field. Due to it’s sensitivity, it is also really good at measuring general vibration.
The Build
Hooking up to the accelerometer is pretty simple. You only need 5 pins: a 3.3V VCC and GND pin, then 3 analog pins for x, y, and z.
Boarduino + IMU wired
Since most accelerometers run on 3.3V, you will likely need to account for that in some way if you are running your arduino at a higher voltage. The x, y, and z pins output a voltage b/w 0 and 3.3V for each axis. For instance, if you are holding the chip completely level, you may see something like 1.65V on the X pin, as you turn it towards 90 degrees, it will increase. You can measure this using analogRead(int pin); .
The Code
The code for this project was so simple it took me under 1 hour to write the firmware and the client software. First the Arduino code:
#define DEBUG 0 //set to 1 for debug mode
#define AX 0 //accelerometer pin X
#define AY 1 //accelerometer pin Y
#define AZ 2 //accelerometer pin Z
#define DELAY 1000
#define SENSITIVITY 3 //movement required to send a movement event
//storage for x,y,z values
int x, y, z;
//storing the old values to compare relation movement
int old_x, old_y, old_z;
//our 'relative motion'
int movementFactor;
void setup() {
Serial.begin(9600);
readSensors();
saveValues();//initialize
}
void loop() {
readSensors();
determineMovementFactor();
if (movementFactor > SENSITIVITY) {
if (DEBUG) {
Serial.print("x = ");
Serial.print(x);
Serial.print(" | y = ");
Serial.print(y);
Serial.print(" | z = ");
Serial.print(z);
Serial.print(" | movement = ");
Serial.println(movementFactor);
} else {
Serial.println(movementFactor);
}
saveValues();
delay(DELAY);
}
}
/**
* All we are doing here is finding the greatest derivative from
* the 3 axises and choosing that as our relative movement factor
*/
void determineMovementFactor() {
movementFactor = abs(old_x - x);
int temp = abs(old_y - y); //due to impl of max(), no inline functions
movementFactor = max(movementFactor, temp);
temp = abs(old_z - z);
movementFactor = max(movementFactor, temp);
}
void readSensors() {
x = analogRead(AX);
y = analogRead(AY);
z = analogRead(AZ);
}
void saveValues() {
old_x = x;
old_y = y;
old_z = z;
}
Assuming you don’t have it in DEBUG mode, it constantly reads the values off the accelerometer and finds the greatest change b/w all the axises. If that is greater than your SENSITIVITY, it sends it to the computer. I chose to do it this way rather than sending the movement factor every second so I could exploit the sparsity of the data and just use a scatter plot to stretch it out over the time domain. Now for the client side…
The client code is just a single python script. To get it working you need to install pyserial of course.
#!/usr/bin/env python
#
# @file: monitor.py
# @author: Benjamin Eckel
#
# To Run -->
# $python monitor.py > output_file.csv
#
import serial
import datetime
FTDI_USB = "/dev/tty.usbserial-A3000RBH"
if __name__ == "__main__":
print "Time,Movement"
ser = serial.Serial(FTDI_USB, 9600)
while 1:
movement = ser.readline()
print "%s,%s" % (datetime.datetime.now().strftime("%H:%M:%S"), movement.rstrip())
This code is self explanatory, you run it and it prints out a comma-separated file with the hours, minutes, seconds, and the relative motion. You can redirect the output to whatever file or program you choose. Here is an example of some output:
All you need to change is your FTDI_USB location. If you don’t know yours, plug in your arduino and run this command:
ls /dev/usb.tty.usb*
This is assuming you are running OS X or Linux of course. Windows uses COM ports. Refer to pyserial for help.
The Installation
To get good readings, it is best to place the sensor somewhere on your mattress that experiences the maximum displacement of the surface (the anti-node) caused by your movements. You can program it in DEBUG mode and try different locations out. I just wrapped the sensor in a sock and stuck it in one of my pillows.
Sock wrapped sensor
The Results
To graph the results, I imported the resulting csv file into excel and created a scatter plot. This is my sleep from June 20 1:30AM to 11:00AM [ took a while to get out of bed :) ]:
Click to view full chart
From what little I know about the subject, I think I should be seeing spikes separated by intervals of about 90 minutes (the time it takes to get a full cycle of sleep). There are definitely some abnormalities in there. I am going attempt to do this for about a week and maybe send the results to a doctor friend.
I just posted up an article on the Gumbolabs blog explaining how to use the Radioshack infrared receiver module. It contains Arduino code along with the explanations.
This article is one in a series of articles about particle systems in Processing. To view the first of the series, click here.
In the last article, we discussed forces and left off with questions about our technique for determining the particle’s position at each frame. It turns out that it was wildly inaccurate. In this article, I am first going to try to explain the mathematics behind what we are doing and use that knowledge to expose the source of our problem. Then I am going to try to present a solution.
Where Does the Problem Come From?
Recall the equation I presented in the last article which gives us distance given the state of a body and a time period:
This equation gives us an exact distance for the given acceleration, initial velocity, and time period and we used it to measure the error of our technique. So the question remains, why not use this equation to determine the particle position at each frame? Well, the problem is unfortunately too complicated for me to explain easily but I am going to try. From my understanding, that equation only works when acceleration is constant over time. If you take a hard look at it, this property will become apparent. The equation gives you a distance for a given step of time but is not actually a continuous function with respect to time. In other words, you cannot use it to solve acceleration for any given point in time. In a particle system, there may be one to many forces acting on a given particle at any time. On top of that, those forces are changing the acceleration of the particle in seemingly erratic ways. In fact, it is said that the function which determines the acceleration with respect to time is well-defined but is actually unknown. Therefore, we cannot get an exact value at any given point, we can only approximate it. Given that our calculations are approximations, we have to assume there will always be an error. All we can do is find techniques to make the error manageable.
Solving Systems Governed by Unknown Functions
So, how exactly were we solving this unknown function before? Let’s step away from our particle system for a very brief moment. You may remember from high school maths that you can solve the integral of a function by integrating that function over an interval when given boundary conditions. You may also remember that in certain situations, such as when f(t, y) is unknown or cannot be solved analytically, you can solve these problems by taking an initial value and solving in steps. This is called an initial value problem and follows this form:
Some tricky manipulating of this formula results in a series which usually looks something like this:
So, I am pretty bad at all this maths o_0 but hopefully you are starting to see something here. In this series, we are determining the next value [y_n+1] using the current value [y_n] and adding the differential function. Just like we did with our particles, we found velocity by adding the last velocity plus the difference in velocity[acceleration]. Then we did the same for position. Acceleration is a second order derivative so when using this method which solves first order diffEQs, we must break it down into two first order derivatives:
where v is velocity and p is position.
Now let’s try and make sense out of this mess.
Euler’s Method for Solving First Order Differential Equatons
You may not have known it, but the method we were using before in our code is known as Euler’s Method and is the defacto, inherent approach to this problem. The only thing I left out is the time step. The actual process should look something like this:
where h is the time step in the simulation. As an astute reader, you may realize that the h was actually accounted for before. It was just assumed to be ’1′ ;) The common way to make the Euler method more accurate is to reduce the time-step [h]. I will explain why in a second, let’s just quickly alter the code to see what happens. Go check it out in action then come back and examine the source:
ArrayList particles;
Particle a;
float h = 0.2; //new time step < 1
void setup() {
size(600, 600);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
particles = new ArrayList();
a = new Particle();
a.location = new PVector(width/2f, height/2f);
a.fixed = true;
particles.add(a);
}
void draw() {
background(255);
for (int i = particles.size()-1; i >= 0; i--) {
Particle p = (Particle)particles.get(i);
// experiment w/ different friction coefficients
applyDissipativeForce(p, 0.01);
applyAttractiveForce(a, p, 5000f, 50f);
if(!p.exist()) {
particles.remove(i);
}
}
}
void mouseDragged() {
particles.add(new Particle());
}
void applyDissipativeForce(Particle p, float friction) {
PVector f = PVector.mult(p.velocity, -friction);
p.applyForce(f);
}
void applyAttractiveForce(Particle a, Particle b, float strength, float minDistance) {
PVector dir = PVector.sub(a.location, b.location);
float d = dir.mag();
if (d < minDistance) d = minDistance;
dir.normalize();
float force = (strength * a.mass * b.mass) / (d * d);
dir.mult(force);
if (!b.fixed) b.applyForce(dir);
if (!a.fixed) {
dir.mult(-1f);
a.applyForce(dir);
}
}
public class Particle {
PVector location;
PVector velocity;
PVector acceleration;
float mass;
boolean fixed;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
acceleration = new PVector(0, 0, 0);
mass = 1;
fixed = false;
}
public boolean exist() {
if (fixed) {
fill(255,0,0);
} else {
velocity.add(PVector.mult(acceleration, h));
location.add(PVector.mult(velocity, h));
acceleration.mult(0);
fill(150);
}
ellipse(location.x, location.y, 10, 10);
return true;
}
void applyForce(PVector force) {
acceleration.add(PVector.div(force, mass));
}
}
All I did was alter the integration part, added the h float at 0.2, and increased the strength to 5000 to speed things up a little. What you should have noticed is that it has gotten a lot slower and it should be obvious why: Our framerate is the same but we are stepping ahead in the simulation in smaller steps of time.
What you probably did not notice is the huge increase in accuracy. To prove this, let’s re-implement our python script experiment from the last article and compare it with the exact results again. This time we will have a time step-size of 0.1:
t = 0
h = 0.1
vel = 0
pos = 0
acc = 10
while (t <= 10):
print "t = %s" % t
print "pos = %s" % pos
print "vel = %s" % vel
print "-------------------"
vel += acc * h
pos += vel * h
t += h
When we run this script we get [middle cut out for brevity]:
t = 0
pos = 0
vel = 0
——————-
t = 0.1
pos = 0.1
vel = 1.0
——————-
t = 0.2
pos = 0.3
vel = 2.0
——————-
t = 0.3
pos = 0.6
vel = 3.0
——————-
t = 0.4
pos = 1.0
vel = 4.0
——————-
t = 0.5
pos = 1.5
vel = 5.0
——————-
.
.
.
——————-
t = 9.6
pos = 465.6
vel = 96.0
——————-
t = 9.7
pos = 475.3
vel = 97.0
——————-
t = 9.8
pos = 485.1
vel = 98.0
——————-
t = 9.9
pos = 495.0
vel = 99.0
——————-
t = 10.0
pos = 505.0
vel = 100.0
Position = 505. You may remember that using these initial values with our equation:
yielded an exact end position of 500. Now we are only off by 5 whereas with a time-step of 1 we were off by 50! This should show you that the error for this method is close to [i think?]:
There is a mathematical way to determine the real error but we aren’t going to go into it because we will be abandoning this method soon and it is too complicated for the little value it provides us. If you really want to know, check out the Wikipedia article on Euler’s method.
Why does a Smaller h Yield More Accurate results?
The question arises, why is it more accurate? The answer may help you better understand what is happening. Let’s look at a graphic representation of the Euler method in action:
Graphical representation of Euler's method (Source: Wikipedia)
Pretend that the blue line is our unknown but well-defined function in which the acceleration is the dependent variable and that A0, A1, A2, A3, and A4 are the acceleration values we calculate. Essentially, at each step, we know the derivative [the tangent to the curve] and multiply it by the time step, then add it to the current A value to guess the next value. At that point, we assume we are at the correct spot and do this again and again and again. You can now see why the error accumulates. We are basically walking a dark path using a slow strobe light. Now imagine that the A values were taken in smaller intervals [the strobe light is blinking faster] and you can see that we would be closer to the blue line [exact solutions]. But if you took a bigger step size, you would be farther away from the exact values. The following pictures illustrate the point perfectly. The red line is the exact solution and the blue is the approximation:
So as h gets smaller, our approximated curve inches closer to the exact curve.
Picking a New Method
Where do we go from here? Well, we have determined that the Euler method is just not really feasible. Sure, it may have worked out for our little examples but we can achieve more stability and accuracy with some other methods. Plus we can’t pump up the framerate enough to get a realistic speed. It just doesn’t look believable. In the next article, I am going to show you how to implement the Runge-Kutta methods, particularly the Fourth Order Runge-Kutta method. I am also going to try to explain the differences between implicit and explicit integration.
This article is one in a series of articles about particle systems in Processing. To view the first of the series, click here.
To recap, in my last article, I went over the basics of particles and Newtonian motion in Processing. In this article, I am going to explain using forces and Newton’s Second Law of Motion to describe the effects of forces on our particles. I am also going to point out that our current integration technique is flawed as segway into the next article.
Forces
When we last left our particle physics sandbox, it looked liked this:
ArrayList particles;
int LIFESPAN = 120;
//global gravity
PVector acceleration = new PVector(0f, 0.025);
void setup() {
size(400, 400);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
particles = new ArrayList();
}
void draw() {
background(255);
//only create when mouse moves
if (abs(mouseX-pmouseX) > 0.0001) {
particles.add(new Particle());
}
for (int i = particles.size()-1; i >= 0; i--) {
Particle p = (Particle)particles.get(i);
if(!p.exist()) {
particles.remove(i);
}
}
}
public class Particle {
PVector location;
PVector velocity;
int age;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
age = 0;
}
public boolean exist() {
velocity.add(acceleration);
location.add(velocity);
ellipse(location.x, location.y, 10, 10);
if (age > LIFESPAN) {
return false;
}
age++;
return true;
}
}
The PVector acceleration represents the acceleration resulting from a gravitational force. Gravity is a force and results in an acceleration of an object. There is an equation to determine gravity but we usually take it to result in a constant acceleration of 9.8 m/s^2 when viewing things on earth. This is because our measurements on earth usually involve relatively small and close things [relative to earth's mass and space] and the difference b/w earth’s effects on these objects is negligible.
Anyways, like all other units of motion, a force can be represented as a Vector. And different forces have different equations to describe their effects. The result a force has on an object is described by Newton’s Second Law of Motion:
where f is the force, m is the object’s mass, and a is acceleration. It is usually more helpful to represent the equation in this alternate form:
Knowing this, we can now calculate the acceleration resulting from a force on an object. First we will get rid of the global PVector acceleration then we will alter the Particle class to look like this:
public class Particle {
PVector location;
PVector velocity;
PVector acceleration;
float mass;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
acceleration = new PVector(0, 0, 0);
mass = 1;
}
public boolean exist() {
velocity.add(acceleration);
location.add(velocity);
ellipse(location.x, location.y, 10, 10);
//when it comes to a rest, get rid of it
if (abs(velocity.x) < 0.001 && abs(velocity.y) < 0.001) {
return false;
}
acceleration.mult(0); //need to clear acceleration
return true;
}
void applyForce(PVector force) {
acceleration.add(PVector.div(force, mass));
}
}
First, note the changes to the Particle class structure:
The particle now keeps track of it’s own acceleration
Each particle has a mass, which we are just setting to 1 for convenience
We added the applyForce(PVector) method
We are now removing a particle from the system when it’s velocity falls below a certain level
Take a look at the applyForce method. There are a few things to note about this method in particular:
Notice that this is basically Newton’s equation in the second form we presented: a = f/m
We are calling add() on acceleration b/c we may be applying other forces to this particle. It may already have an acceleration value set. But we are also clearing the acceleration after each iteration [line 23]. This is important b/c we don’t want the force to have an accumulative effect. See what happens when you comment out that operation.
Not that this is a Java lesson, but notice that I used the static PVector.div(PVector, float) method rather than using an instance method. The static methods return the results of the operations as new instances whereas the instance methods modify the PVector that the reference points to. Pay attention to this detail, when to use either should become obvious if not already.
Applying a Force: Friction
Now we want to actually apply a force to a particle. There are many types of forces and many equations to describe their behavior. By far, one of the easiest to understand is a dissipative force. It most commonly manifests itself in day to day life in the form of friction but can also be used to simulate other phenomena. The equation is:
where f is the resulting force, C is a mathematical constant known as the “coefficient of friction“, and v is the object’s velocity. Obviously, the larger C is, the more friction is applied and the quicker your object will come to a state of rest. It is usually between 0.0 and 1.0 . To implement this, I added this global function:
Pretty self-explanatory. We take the velocity of the particle that we are applying this force to and we multiply it by the coefficient * -1. Let’s put all this together. Check out the code in action, then come back and examine the whole source:
ArrayList particles;
void setup() {
size(400, 400);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
particles = new ArrayList();
}
void draw() {
background(255);
//only create when mouse moves
if (abs(mouseX-pmouseX) > 0.0001) {
particles.add(new Particle());
}
for (int i = particles.size()-1; i >= 0; i--) {
Particle p = (Particle)particles.get(i);
// experiment w/ different friction coefficients
applyDissipativeForce(p, 0.08);
if(!p.exist()) {
particles.remove(i);
}
}
}
void applyDissipativeForce(Particle p, float friction) {
PVector f = PVector.mult(p.velocity, -friction);
p.applyForce(f);
}
public class Particle {
PVector location;
PVector velocity;
PVector acceleration;
float mass;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
acceleration = new PVector(0, 0, 0);
mass = 1;
}
public boolean exist() {
//when it comes to a rest, get rid of it
if (abs(velocity.x) < 0.001 && abs(velocity.y) < 0.001) {
return false;
}
velocity.add(acceleration);
location.add(velocity);
ellipse(location.x, location.y, 10, 10);
acceleration.mult(0);
return true;
}
void applyForce(PVector force) {
acceleration.add(PVector.div(force, mass));
}
}
Applying Another Force: Attraction
Now we are going to implement a more complicated force, attraction. This force could be used to simulate gravitational force, without all that Einstein, space-time continuum craziness of course. The equation we are going to use is:
where f is the force, G is the gravitational constant [defines the apparent 'strength'], m1 and m2 are the masses of the two bodies in effect, and d is the distance between the two bodies. The problem is that this equation computes the magnitude, or the power, of the force; but not the direction. With the dissipative force, the direction was the opposite of the velocity. This was obvious in the equation. The force in this situation however is between the 2 bodies, so naturally, the equation will involve subtracting the positions. Here it is:
where d is direction and L1 and L2 are the position vectors of the bodies.
Now we will implement the global function to apply an attractive force on two particles. We will also have to modify the Particle class:
void applyAttractiveForce(Particle a, Particle b, float strength, float minDistance) {
PVector dir = PVector.sub(a.location, b.location);
float d = dir.mag();
if (d < minDistance) d = minDistance;
dir.normalize();
float force = (strength * a.mass * b.mass) / (d * d);
dir.mult(force);
if (!b.fixed) b.applyForce(dir);
if (!a.fixed) {
dir.mult(-1f);
a.applyForce(dir);
}
}
public class Particle {
PVector location;
PVector velocity;
PVector acceleration;
float mass;
boolean fixed;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
acceleration = new PVector(0, 0, 0);
mass = 1;
fixed = false;
}
public boolean exist() {
if (fixed) {
fill(255,0,0);
} else {
velocity.add(acceleration);
location.add(velocity);
acceleration.mult(0);
fill(150);
}
ellipse(location.x, location.y, 10, 10);
return true;
}
void applyForce(PVector force) {
acceleration.add(PVector.div(force, mass));
}
}
I added a fixed property. If the Particle is fixed, it is not affected by outside forces. They also render red.
As for the applyAttractiveForce(), for the most part, it follows our two equations above. The first thing to note is that we are constraining the distance:
if (d < minDistance) d = minDistance;
This is because too small a distance value will create too large a force and the particles will go flying off. Also notice that we are applying “equal but opposite” forces on each particle:
if (!b.fixed) b.applyForce(dir);
if (!a.fixed) {
dir.mult(-1f);
a.applyForce(dir);
}
ArrayList particles;
Particle a;
void setup() {
size(800, 600);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
particles = new ArrayList();
a = new Particle();
a.location = new PVector(width/2f, height/2f);
a.fixed = true;
a.mass = 5.0;
particles.add(a);
}
void draw() {
background(255);
for (int i = particles.size()-1; i >= 0; i--) {
Particle p = (Particle)particles.get(i);
// experiment w/ different friction coefficients
applyDissipativeForce(p, 0.01);
//experiment with strength and minDistance, try a minDistance below 30
applyAttractiveForce(a, p, 500f, 50f);
if(!p.exist()) {
particles.remove(i);
}
}
}
void mouseDragged() {
particles.add(new Particle());
}
void applyDissipativeForce(Particle p, float friction) {
PVector f = PVector.mult(p.velocity, -friction);
p.applyForce(f);
}
void applyAttractiveForce(Particle a, Particle b, float strength, float minDistance) {
PVector dir = PVector.sub(a.location, b.location);
float d = dir.mag();
if (d < minDistance) d = minDistance;
dir.normalize();
float force = (strength * a.mass * b.mass) / (d * d);
dir.mult(force);
if (!b.fixed) b.applyForce(dir);
if (!a.fixed) {
dir.mult(-1f);
a.applyForce(dir);
}
}
public class Particle {
PVector location;
PVector velocity;
PVector acceleration;
float mass;
boolean fixed;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
acceleration = new PVector(0, 0, 0);
mass = 1;
fixed = false;
}
public boolean exist() {
if (fixed) {
fill(255,0,0);
} else {
velocity.add(acceleration);
location.add(velocity);
acceleration.mult(0);
fill(150);
}
ellipse(location.x, location.y, 10, 10);
return true;
}
void applyForce(PVector force) {
acceleration.add(PVector.div(force, mass));
}
}
Hopefully this code is pretty easy to understand. Particle a is the center red particle. We set it’s location to the center of the screen, set it as fixed, then give it a bigger mass. For every particle, we are applying both the dissipative force and the attraction to particle a:
If you took the time to closely examine our technique for determining a particle’s position at each frame, you may have noticed some inaccuracies. Let’s use an analogy. Switch over your thought process from 2d down to 1d. Let’s say we are in a car travelling down a road [is commonly described as if it is one dimension]. If we know the car’s acceleration, and we assume that it is constant, we can theoretically calculate the position after a given time with the equation:
where s is the distance, u is the initial velocity, a is the acceleration, and t is time. So if we have an initial velocity of 0, an acceleration of 10 meters per second every second and we travel for 10 seconds, we end up with a distance of 500 meters. Now let’s try this with our methodology. Here is a python script to simulate our technique:
t = 0
vel = 0
pos = 0
acc = 10
while (t <= 10):
print "t = %s" % t
print "pos = %s" % pos
print "vel = %s" % vel
print "-------------------"
vel += acc
pos += vel
t += 1
If we run this, we get:
t = 0
pos = 0
vel = 0
——————-
t = 1
pos = 10
vel = 10
——————-
t = 2
pos = 30
vel = 20
——————-
t = 3
pos = 60
vel = 30
——————-
t = 4
pos = 100
vel = 40
——————-
t = 5
pos = 150
vel = 50
——————-
t = 6
pos = 210
vel = 60
——————-
t = 7
pos = 280
vel = 70
——————-
t = 8
pos = 360
vel = 80
——————-
t = 9
pos = 450
vel = 90
——————-
t = 10
pos = 550
vel = 100
——————-
Position = 550!?!? Why is it so far off? After 10 seconds in the simulation we have already drifted off 50 meters from where we should be. Another problem we see is that this error is accumulating [or compounding]. So as time goes on, our calculations of the particle’s positions will become farther and farther away from the truth. The solution to this problem requires a closer examination of the mathematics behind our process and Newton’s Second Law in general. In the next article, I am going to attempt to explain how this error occurs and how we can reduce it to a manageable level.
I have recently decided to take on the task of creating a particle system from scratch designed to work specifically with Processing. I have always been interested in the subject and I wanted to challenge myself and share something with the Processing community. I am going to write this series of articles to document my learning process as well as the evolution of this library.
For this first article, I am going to start with the most basic concepts of dealing with particles in Processing. No differential equations, no numerical methods, just simple linear algebra and Java. This one may be pretty boring to most Processing users, but I want to start off simple to establish a good base of knowledge because it will get pretty hairy later. If you are a complete beginner, I suggest you start here and come back when you have a decent understanding of Processing.
The “Hello World!” of particle physics
The first thing we want to examine is how to draw and move a single particle. First check out this code in action here then come back and closely examine the source:
//particle's screen coordinates
float x, y;
//particle's velocity
//velocity is change in postion
float vx, vy;
void setup() {
size(400, 400);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
//set position and velocity
resetParticle();
}
void draw() {
background(255);
//draw the particle with diameter of 10
ellipse(x, y, 10, 10);
/*
* update the position by adding on the velocity
* remember that velocity means incremental change
* in position.
*/
x += vx;
y += vy;
}
// on a mouse click
void mousePressed() {
resetParticle();
}
void resetParticle() {
//reset the particle to the middle of the screen
x = y = width/2f;
//set new random velocity
// b/w - and + numbers
vx = random(-0.39, 0.1);
vy = random(-0.1, 0.21);
}
Let’s quickly dissect what is going on here. First we establish the two floats x and y to be the center point of the particle in screen coordinates. Hopefully you have messed with Processing or computer animation enough to understand screen coordinates. Then we establish vx and vy as the velocity of our particle. As we remember in high school physics, velocity is defined as the change in position over time. In computer animation, our general unit of time is a “frame”. As you should know, the draw() function is called at the framerate that you run the program at [Processing default is 60 times a second I think]. So every time that function is called, we redraw a white background drawing over any previous frames:
background(255);
then we draw the particle as an ellipse at it’s current x, y position:
ellipse(x, y, 10, 10);
then we update the position one time step:
x += vx;
y += vy;
and this is done 60 times a second creating the illusion of motion due to persistence of vision.
Adding Acceleration
You may also remember from your childhood science classes that acceleration is the change in velocity over time. Quickly view this code in action and take note of the now non-uniform speed of the particle:
float x, y;
float vx, vy;
//acceleration
float ax, ay;
void setup() {
size(400, 400);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
//set position and velocity
resetParticle();
}
void draw() {
background(255);
//draw the particle with diameter of 10
ellipse(x, y, 10, 10);
/*
* update the velocity by adding on the acceleration
* remember that acceleration means incremental change
* in velocity.
*/
vx += ax;
vy += ay;
/*
* update the position by adding on the velocity
* remember that velocity means incremental change
* in position.
*/
x += vx;
y += vy;
}
// on a mouse click
void mousePressed() {
resetParticle();
}
void resetParticle() {
//reset the particle to the middle of the screen
x = y = width/2f;
//set velocity to 0
vx = vy = 0f;
ax = random(-0.01, 0.01);
ay = random(0.003, 0.01);//downward
}
Now that we are adding on the acceleration to the velocity each frame, the velocity increases as time goes on and the particle moves faster and faster. Kind of like gravity (an acceleration of 9.8 m/s2) or a car accelerating faster and faster. Congratulations! We now have the ability to crudely describe Newtonian motion!
Abstraction and Multiple Particles
Now that we can create one particle with ease, let’s create many. Basically what we need to do is create an abstraction of what a particle is then turn that into a Java class. This isn’t a Java lesson and I hope you know what that means so let’s jump straight into it. Check out this program then come back and learn about the code.
Particle[] particles;
int NUM_PARTICLES = 100;
void setup() {
size(400, 400);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
}
void draw() {
background(255);
if (particles != null) {
for (int i = 0; i < NUM_PARTICLES; i++) {
particles[i].exist();
}
}
}
void mousePressed() {
particles = new Particle[NUM_PARTICLES];
for (int i = 0; i < NUM_PARTICLES; i++) {
particles[i] = new Particle();
}
}
public class Particle {
float x, y;
float vx, vy;
float ax, ay;
public Particle() {
x = mouseX;
y = mouseY;
vx = vy = 0f;
ax = random(-0.2, 0.2);
ay = random(-0.08, 0.05);
}
public void exist() {
vx += ax;
vy += ay;
x += vx;
y += vy;
ellipse(x, y, 10, 10);
}
}
Of course this is too crazy to resemble anything real in the our physical world, but it shows the point. Basically what we have done is created a Particle class which contains references to it’s own state information. Then we can have an array of Particles and update those in a nice loop. There is still some more improvement that can be done to this though. Let’s talk about vectors.
Vectors
A Vector is essentially a 1 row N column matrix. In computer graphics, we use the Vector data structure to represent multi-dimensional state information such as position, velocity, and acceleration. Because they are matrices and we set them to have the same dimensions, we can use standard matrix operations and concepts from linear algebra to help us solve some of the systems of linear equations we face. Processing contains it’s own implementation of a vector called a “PVector”. It is a 3-dimensional vector so we can represent things in 2d or 3d space. Check out this re-write of the last program in action and I will explain how it works after the code:
ArrayList particles;
int LIFESPAN = 120;
//global gravity
PVector acceleration = new PVector(0f, 0.025);
void setup() {
size(400, 400);
stroke(0);
strokeWeight(3);
fill(150);
smooth();
particles = new ArrayList();
}
void draw() {
background(255);
//only create when mouse moves
if (abs(mouseX-pmouseX) > 0.0001) {
particles.add(new Particle());
}
for (int i = particles.size()-1; i >= 0; i--) {
Particle p = (Particle)particles.get(i);
if(!p.exist()) {
particles.remove(i);
}
}
}
public class Particle {
PVector location;
PVector velocity;
int age;
public Particle() {
location = new PVector(mouseX, mouseY);
//get velocity from direction and speed of mouse movement
velocity = new PVector(mouseX-pmouseX, mouseY-pmouseY);
age = 0;
}
public boolean exist() {
velocity.add(acceleration);
location.add(velocity);
ellipse(location.x, location.y, 10, 10);
if (age > LIFESPAN) {
return false;
}
age++;
return true;
}
}
Notice that I changed over the x and y variables to the PVectors and this gives us operations like add, subtract, multiply, magnitude, normalize, and etc. Check out this article on matrix math if you don’t remember it from grade school. In this example, all you need to know is that when you add two PVectors, it does a matrix add. I do, however, highly suggest that you get real familiar with matrix math and systems of equations if you want to move on:
velocity.add(acceleration);
//is equivalent to
velocity.set( acceleration.x+velocity.x,
acceleration.y+velocity.y,
acceleration.z+velocity.z );
I also moved the acceleration vector out to have global static scope. And I am getting the velocity from the changing mouse positions. Remember that velocity is change in position, so we can get the velocity vector by subtracting the new cursor position from the old.
I don’t really want to go into Java stuff but notice that I also changed the array to an ArrayList. And if you are wondering why I didn’t use a parameterized ArrayList, it is b/c Processing uses a light-weight open source java compiler that doesn’t support it.
Until next time
That is it for now. I hope that wasn’t too boring. In the next article, I am going to talk more about vectors, different kinds of forces. I will also briefly point out the inaccuracies in our integration technique and open up a huge can of maths worms.
