May 26th, 2010 — Code, Java, Processing, Software
I have finally released p5particles which is my particle systems library for Processing. You can view the project here:
http://code.google.com/p/p5particles/
There is currently limited documentation and support. I will start improving it if others find it useful.
March 18th, 2010 — Code, Java, Processing, Software
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.
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.
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.
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:
velocity.add(PVector.mult(acceleration, h)); location.add(PVector.mult(velocity, h));
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.
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:
big step of h
medium sized step of h
small step of h
tiny step of h
Here is the source of these photos.
So as h gets smaller, our approximated curve inches closer to the exact curve.
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.
Gaffer on Games : Integration Basics
March 17th, 2010 — Code, Java, Processing, Software
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.
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:
Take a look at the applyForce method. There are a few things to note about this method in particular:
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:
void applyDissipativeForce(Particle p, float friction) {
PVector f = PVector.mult(p.velocity, -friction);
p.applyForce(f);
}
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));
}
}
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);
}
Now let’s look at the whole thing. Check out this code in action then come back and view the source:
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:
applyDissipativeForce(p, 0.01);
applyAttractiveForce(a, p, 500f, 50f);
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.
Click here to view the next article of this series —>
Daniel Shiffman : The Nature of Code
March 12th, 2010 — Code, Java, Processing, Software
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 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.
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!
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.
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.
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.
May 2nd, 2009 — Arduino, Code, Hardware, Java, Leopard, Processing, Python, Software, hack
A while back, I wrote a small python script to act as a proxy between IRC and an Arduino. This video illustrates some of the details.
Here is a zip file with everything to download.
The basic idea was this, the python IRC bot sits in the chat room and waits until someone addresses it’s nick like so:
{NICK}: {command} {arg1} {arg2} ... {etc}
The bot parses the response into a command and a series of arguments. Originally, all this was hardcoded, but through python’s powers of introspection, I realized I could make this way more dynamic. The end result is a framework of sorts. The user can edit the Commands.py file adding functions and settings to his/her specification. The example in the video would be written like this:
#
# module for configuration and commands
#
from Arduinos import Arduino
# Program Settings
DEBUG = True
# IRC Settings
NICK = "arduino"
SERVER = "irc.paraphysics.net"
CHANNEL = "#arduinoroom"
PORT = 6667
# Arduino Settings
USB_PATH = '/dev/tty.usbserial-A7006Qe8'
BAUD = 9600
# define methods
def lightLed(arduino, args):
arduino.send('~') # header
arduino.send(args[0])
arduino.send('~') # terminating
return arduino.read(4) # read 4 bytes from arduino
def readPot(arduino):
arduino.send('}}') # header
return arduino.read(4) # read 4 bytes from arduino
Obviously, you can change the constants to what you need and the daemon picks up on these when it is started. The nice thing about this ‘framework’ is the second half of the script. When you define a method, the IRC bot automagically “understands” it. This is because every time it receives a message directed to it, it calls
reload(Commands)
and reads the function names. So you don’t have to restart the server while developing your functions [An idea stolen from Rails]! Then it parses the message coming in and tries to call the first word as a method and parses the rest of the message into a list. There is no need to worry about the details, what you do need to know to use it is that there are two types of functions, ones with and ones without arguments. The above example illustrates both. Let’s say we define 2 functions:
# define methods
def function1(arduino, args):
print "function1"
print args
def function2(arduino):
print "function 2"
print "no args"
Now, we go into the IRC chat room and issue a command like this:
arduino: function2
The output would be this:
function 2 no args
Then lets say we send this:
arduino: function1 hello world!
The bot would parse it and you would see this output:
function1 ['hello', 'world!']
so args becomes a list by splitting the rest of the statement between whitespaces. Your function always needs to have the ‘arduino’ argument. This argument is a custom class I created which looks like this:
import serial class Arduino(): def __init__(self, path='/dev/tty.usbserial', baud=9600): self.ser = serial.Serial(path, baud) def send(self, data): self.ser.write(data) def read(self, bytes): while (1): if (self.ser.inWaiting() > bytes-1): return self.ser.read(bytes) def flush(): self.ser.flushInput()
It is pretty simple, send() sends a string. read() waits for the defined number of bytes to come in and returns the results. To better understand how the example Commands.py script works, take a look at the arduino sketch:
#define LED 13
void setup() {
pinMode(LED, OUTPUT);
Serial.begin(9600);
}
void loop() {
if (nextByte() == 126) { // header byte ('~' character) led command
char args[] = {0,0,0,0,0,0,0,0,0,0};
char charIn = 0;
byte i = 0;
while (charIn != 126) { // wait for header byte again
charIn = nextByte();
args[i] = charIn;
i += 1;
}
if ((args[0] == 'o') && (args[1] == 'n')) {
digitalWrite(LED, HIGH);
Serial.print("on ");
}
else if ((args[0] == 'o') && (args[1] == 'f')) {
digitalWrite(LED, LOW);
Serial.print("off ");
}
delay(10);
Serial.flush();
}
else if (nextByte() == 125) { // header byte pot command
int val = analogRead(0);
if (val < 10) {
Serial.print(val);
Serial.print(" ");
}
else if (val < 100) {
Serial.print(val);
Serial.print(" ");
}
else if (val < 1000) {
Serial.print(val);
Serial.print(" ");
}
else {
Serial.print(val);
}
}
delay(10);
//if(Serial.available() > 0) {
Serial.flush();
//}
}
byte nextByte() {
while(1) {
if(Serial.available() > 0) {
byte b = Serial.read();
return b;
}
}
}
Yeah, not the cleanest code but hopefully you get the idea. I am not releasing a whole lot of detail on how to use this because I figure If you are using it, then you already know enough about python and Arduino to get by. I have created this to allow for support for firmata but have yet to implement it and probably never will until/if people beg. To run, first get pyserial and irclib.py. Upload the sketch to your arduino and run:
python Main.py
May 1st, 2009 — Code, Java, Processing, Software
Here is a simple sketch I wrote as a joke for a friend. Don’t make fun of the code, I wrote it quick!! Download it here. Download the zip for everything. Probably app only works for os x with iSight. srry.
November 25th, 2007 — Code, DSP, Java, Processing, Software
If you are looking use the FFT algorithm in a java application, you should probably take a look around for a well optimized and tested library. But if you want to get into it yourself, here is a start. I ported this code from C and I got that from “Numerical Recipes in C“. It was originally written by N. M. Brenner in C. This book is essential for anyone interested in DSP programming. I got it for a Numerical Methods class I took last year and I have still been finding uses for it.
For those of you who don’t know, the FFT algorithm allows you to view a discrete function (such as a sound wave) in the frequency domain. This means that you can take a buffer of sampled sound (or function) and retrieve the level of each frequency spectrum. This particular function returns a float array of the levels at each frequency bin b/w 0Hz and the Nyquist frequency. The first half of the array is the real part and the second is the imaginary. Here is a visual representation, the x axis is frequency bins, and the y axis is pressure (or energy):
This is not a very good picture, but, this is a 1 kHz sine wave. Here is the function:
public float[] four1(float data[], int nn, int isign) {
int i, j, n, mmax, m, istep;
float wtemp, wr, wpr, wpi, wi, theta, tempr, tempi;
n = nn << 1;
j = 1;
for(i=1; i
if(j > i) {
float temp;
temp = data[j];
data[j] = data[i];
data[i] = temp;
temp = data[j+1];
data[j+1] = data[i+1];
data[i+1] = temp;
}
m = n >> 1;
while(m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
}
mmax = 2;
while(n > mmax) {
istep = (mmax << 1);
theta = isign*(6.28318530717959/mmax);
wtemp = sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
for(m=1; m
for(i=m;i<=n;i+=istep) {
j = i+mmax;
tempr = wr*data[j]-wi*data[j+1];
tempi = wr*data[j+1]+wi*data[j];
data[j] = data[i] - tempr;
data[j+1] = data[i+1] - tempi;
data[i] += tempr;
data[i+1] += tempi;
}
wr = (wtemp=wr)*wpr-wi*wpi+wr;
wi = wi*wpr+wtemp*wpi+wi;
}
mmax = istep;
}
return data;
}
data[] is the discrete function to be analyzed, nn is the number of complex points. This is generally data[].length/2, and isign is for inversion. Here is some code you can us to do this in Processing:
import krister.Ess.*;
AudioInput myInput;
float gain = 10.;
int bufferSize;
void setup() {
size(512, 200);
Ess.start(this);
bufferSize=512;
myInput = new AudioInput(bufferSize);
myInput.gain(gain);
frameRate(25);
myInput.start();
noFill();
stroke(255, 255, 255);
smooth();
}
void draw() {
background(0, 0, 255);
float[] buffer=four1(myInput.buffer2,
myInput.buffer2.length/4,1);
stroke(255, 255, 255);
for(int i=0; i
line(i, height, i, height-(buffer[i]*1.4));
}
stroke(255, 0, 0);
line(width/2, 0, width/2, height);
}
Of course, you need to drop the four1 function somewhere in there. This is a little crude for sound, but it will allow you to distinguish the difference in pitch of sounds with some accuracy.