Years ago, when I was still in college and had lots of free time, I played around with physics simulations of various kinds. One that most programmers try at some point is a classic gravitational simulation. It’s great fun to setup complex arrangements of planets and stars, turn on the simulation, and watch star systems evolve.

## I am become death, the destroyer of worlds

After trying numerous generic simulations, I decided to try something a bit more complex: A collision simulator. I altered the force equations for the particles, so that when they got close enough together, they would repulse instead of attract. This caused particles to stick together in clumps. By simulating small systems, then saving the results, and throwing the clumps into more clumps, I eventually worked up to a planetery collision simulation with over 1,000 particles. I exported the results to the POV-Ray renderer, and entered the animation into the Internet Ray Tracing Competition.

## Is this simulation *ever* going to finish?

Recently, I started thinking about better ways to do animations like that. The major issue I ran into when creating the full planetary collision was slowness. A gravity simulation with *n* particles needed to calculate *n*^{2} forces every frame in order to update the particles’ positions. This becomes slow very quickly when you increase the number of particles. If you’re doing scientific simulations, then this is probably the best you can do since accuracy is more important than speed. However, if your goal is to make a “cool” animation, then you can get away with some “cheating” as long as it’s not too drastic.

## From up here, they all look like little ants

I’ve done some recent work with kd-trees as a way to create a more efficient nearest neighbor algorithm. This data structure is useful for organizing points in space so that you can search more efficiently in multiple dimensions. Another similar data structure is the octree which has a more regular structure than the kd-tree.

Since my goal is aesthetic rather than scientific, I can speed up force calculations by grouping all of the particles into an octree. Since a group of points becomes more and more like a single point the further away they get, I can get the force from a particular octree node by averaging the positions of the particles inside the node and treating the node as a single particle. Thus, the force from a large group of points far away from another point can be calculated with a single operation instead of many.

Although complete accuracy isn’t necessary for a purely visual effect, I don’t want the simulation to be *too* unrealistic. So, before I start implementing an octree-based simulation, I need to know how much “cheating” I can get away with. To do this, I wrote a PHP script (pasted below) that creates 2x2x2 cubes of 10 random masses at various distances from the origin, and then calculates the magnitude and angle error of the averaged force compared to the actual force.

## Results

As shown in the graphs below, the accuracy is acceptable at surprisingly low distances. When the edge of the cube is as close to the origin as it is to the center of the cube, the magnitude error is less than 10%, and the angle error is less than 6 degrees. In order to bring the errors down to 1% and 1 degres respectively, the distance of the cube from the origin needs to be only 2.5 times the side length of the cube! So, even if there are 100,000 points in the scene, a point which is even somewhat isolated from the rest of the points may only need a handful of force calculations to determine its next position.

Errors for a 2x2x2 cube of 10 points at a distance X from the origin

This PHP code snippet is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

class Vector
{
public $x;
public $y;
public $z;
public function __construct($x = null, $y = null, $z = null)
{
if($x instanceof Vector)
{
$z = $x->z;
$y = $x->y;
$x = $x->x;
}
$this->set($x,$y,$z);
}
public function set($x,$y,$z)
{
$this->x = $x;
$this->y = $y;
$this->z = $z;
}
public function length()
{
return sqrt($this->length2());
}
public function length2()
{
return self::dot($this,$this);
}
public function normalize()
{
$len = $this->length();
$v = new Vector($this);
$v->x /= $len;
$v->y /= $len;
$v->z /= $len;
return $v;
}
public function scaleTo($len)
{
$v = $this->normalize();
$v->x *= $len;
$v->y *= $len;
$v->z *= $len;
return $v;
}
public function scaleBy($len)
{
$v = new Vector($this);
$v->x *= $len;
$v->y *= $len;
$v->z *= $len;
return $v;
}
public static function add($v1,$v2)
{
$v = new Vector($v1);
$v->x += $v2->x;
$v->y += $v2->y;
$v->z += $v2->z;
return $v;
}
public static function sub($v1,$v2)
{
$v = new Vector($v1);
$v->x -= $v2->x;
$v->y -= $v2->y;
$v->z -= $v2->z;
return $v;
}
public static function dot($v1,$v2)
{
return
$v1->x * $v2->x +
$v1->y * $v2->y +
$v1->z * $v2->z;
}
public static function cross($v1,$v2)
{
return new Vector($v1->y*$v2->z - $v1->z*$v2->y, $v1->z*$v2->x - $v1->x*$v2->z, $v1->x*$v2->y - $v1->y*$v2->x);
}
public static function angle($v1,$v2)
{
$len2 = $v1->length() * $v2->length();
$cosTheta = self::dot($v1,$v2) / $len2;
$sinTheta = self::cross($v1,$v2)->length() / $len2;
return atan2($sinTheta,$cosTheta) * 180.0 / 3.1415936535897;
}
}
function random_float ($min,$max) {
return ($min+lcg_value()*(abs($max-$min)));
}
class Mass
{
public $mass;
public $location;
const G = 100.0;
public function force()
{
return $this->location->scaleTo(self::G * $this->mass / $this->location->length2());
}
public static function average($massList)
{
$v = new Vector(0,0,0);
$mass = 0.0;
foreach($massList as $item)
{
$v = Vector::add($v,$item->location->scaleBy($item->mass));
$mass += $item->mass;
}
$v = $v->scaleBy(1.0 / $mass);
$m = new Mass;
$m->mass = $mass;
$m->location = $v;
return $m;
}
public static function random($vOffset, $offsetVariance, $minMass, $maxMass)
{
$m = new Mass();
$m->mass = random_float($minMass, $maxMass);
$m->location = new Vector(
random_float(-$offsetVariance,$offsetVariance),
random_float(-$offsetVariance,$offsetVariance),
random_float(-$offsetVariance,$offsetVariance)
);
$m->location = Vector::add($m->location,$vOffset);
return $m;
}
}
function getResults($offset)
{
$offset = new Vector($offset,0,0);
$variance = 1.0;
$min = 1.0;
$max = 10.0;
$numMasses = 10;
$masses = array();
$forceActual = new Vector();
for($i=0;$i<$numMasses;$i++)
{
$mass = Mass::random($offset,$variance,$min,$max);
$masses[] = $mass;
$forceActual = Vector::add($forceActual,$mass->force());
}
$forceAverage = Mass::average($masses)->force();
$magDiff = $forceActual->length() - $forceAverage->length();
$magDiffRatio = 100.0 * $magDiff / $forceActual->length();
$angle = Vector::angle($forceActual,$forceAverage);
return array(
'forceActual' => $forceActual,
'forceAverage' => $forceAverage,
'magDiff' => $magDiff,
'magDiffRatio' => $magDiffRatio,
'angle' => $angle
);
}
echo "<pre>";
$iterations = 20;
$limit = 20;
$increment = 0.05;
$results = array();
for($offset=0.1; $offset<=$limit; $offset+=$increment)
{
$magDiffAvg = 0.0;
$angleAvg = 0.0;
for($i=0; $i<$iterations; $i++)
{
$curResults = getResults($offset);
$magDiffAvg += abs($curResults['magDiffRatio']);
$angleAvg += abs($curResults['angle']);
}
$magDiffAvg /= $iterations;
$angleAvg /= $iterations;
$results["$offset"] = array(
'magDiffAvg' => $magDiffAvg,
'angleAvg' => $angleAvg
);
}
foreach($results as $offset => $data)
{
extract($data);
echo "{$offset}, {$magDiffAvg}, {$angleAvg}\n";
}