generalrelativity

This article walks through my implementation of modeling orbit with FOAM. The content lays a groundwork which could be used to simulate (a ridiculously simplified version of) our solar system for instance. It has 2 main goals:

1. Give the reader a solid and quick understanding of how one might use FOAM specifically.
2. Exemplify the difference between differential equation solvers.

Orbit Demo
FOAM source

When wanting to create physically realistic animation using FOAM, the first step will often be deciding how to develop the forces required to create the desired motion. Forces in FOAM are primarily handled via modular implementations of IForceGenerator- which simply dictates that the method generate generates a force for its passed element. Gravity in the context of our atmosphere and friction are 2 force generators included in the source.

But how would we go about modeling the force an extremely massive object exerts on another?

Firstly, we need to define an equation which ultimately yields a force that we can apply. Newtonian gravity offers the simple force equation:

• F is the force we’re solving for
• G is a gravitational constant
• m1/m2 are the masses of the bodies in question
• r is the distance between the bodies

It should be clear from looking at this equation that the larger the mass, the shorter the distance, the greater the force- which is about as complex as we need to go to get realistic looking results. One simplification we’ll make is assuming that the body exerting the gravitational force is so much larger than the body being acted on, that we can neglect the force acting back on it.

This is simple enough, but we need to translate Newton’s equation into a form our physics engine can execute- We need to create a new force generator. We’ll look at what’s going on in its constructor and generate method; the class in its entirety can be found here.

First, the constructor takes the element we want to exert the force and gravitationalConstant- this is a magic number that should be tweaked per simulation:

1
2
3
4
5

public function GravitationalForceGenerator( source:ISimulatable, gravitationalConstant:Number = 1.2 )
{
   this.source = source;
   g = gravitationalConstant;
}

I’ve left out some error checking for the sake of brevity here- reference the complete class. Next we implement generate:

1
2
3
4
5
6
7

public function generate( element:ISimulatable ) : void
{
   //find the difference vector
   var diff:Vector = source.position.minus( element.position );
   //add our solved force along the difference vector to the supplied element
   element.addForce( diff.getUnit().times( g * source.mass * element.mass / diff.dot( diff ) ) );
}

Line 6 contains our equation for universal gravitation as discussed above. The change worth noting here is that the force is applied in the direction of the exerting body (GravitationalForceGenerator.source). We do this by normalizing the bodies’ difference vector, and scaling it by the force found in Newton’s equation. Note also that the dot product of a vector and itself is its squared magnitude- this is how we’re finding the square of the bodies’ distance- the equation’s right-hand side denominator.

Now let’s look at how I’ve used this force generator in the context of a FOAM application. Remember that we’re also going to examine the difference between IODESolvers.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

//create a source of the gravitational pull
var source:Circle = new Circle( 400, 300, 60, 10000 );
 
//create an element influenced by source's gravitational pull to be integrated with the Euler IODESolver
var eulerOrbital:Circle = new Circle( 700, 300, 30, 100, 0, -7 );
//create an element influenced by source's gravitational pull to be integrated with the RK4 IODESolver
var rungaKuttaOrbital:Circle = new Circle( 700, 300, 30, 100, 0, -7 );
 
//create the gravitational force generator
var gravity:GravitationalForceGenerator = new GravitationalForceGenerator( source );
 
//add this to each orbital
eulerOrbital.addForceGenerator( gravity );
rungaKuttaOrbital.addForceGenerator( gravity );
 
//add each element to FOAM (no collisions)
foam.addElement( source, false, true, { color:SimpleOrbit.SOURCE_COLOR, alpha:1 } );
//specify the use of our Euler solver
foam.addElement( eulerOrbital, false, true, { color:SimpleOrbit.EULER_COLOR }, new Euler( eulerOrbital ) );
//RK4 is default, no need to explicitly set
foam.addElement( rungaKuttaOrbital, false, true, { color:SimpleOrbit.RK4_COLOR } );

The important aspects here are:

• we create a new GravitationalForceGenerator and add it as a force generator to each of our orbiting bodies, setting our source Circle as its source.
• the orbitals are created with identical masses, positions and initial velocities (tangent to desired orbit at location)
• one orbital gets added with FOAM’s default solver, RK4- the other with an explicitly defined Euler solver.

The whole SimpleOrbit class can be found here.

If you take a look at the demo, you’ll notice that slowly but surely the orbital being integrated with our Euler solver drifts out of desirable orbit. It will continue to spiral outward. If we set the number of iterations we allow the engine to use per frame to 1, you would see that this error accumulates faster. If instead, we set the iterations to 100 for instance, its orbit would match the more precise RK4-integrated orbital much longer.

So why is this happening? To better explain, I’ll steal some images and description from an earlier post of mine- reference the Euler and RK4 classes also.

At every step in time in which we evaluate state and derivative, the forces acting on our orbiting bodies is different. If we were able to do this evaluation an infinite number of times every frame, Euler would provide the correct answer. This is obviously impossible and so even though the forces acting on the body continually change, we’re limited to apply the force found at a single moment in time across the entire interval over which we’re integrating.

The line segment piercing the starting point represents the derivative. Because Euler simply advances the equation via this derivative, it should be easy to see why it is so error prone. Compare this with how RK4 works:

Because it uses the slope at multiple points over the interval to evaluate the ultimate derivative, our approximation using RK4 is entirely more accurate. This accuracy comes at a cost- RK4 is very much more expensive to calculate- although if you tried the different iteration suggestion I made, you’d see it is more accurate at 3 iterations than Euler is at 100. This is one of FOAM’s strong points- it lets you choose the solver best suited for a specific task within a simulation.

Hopefully this has been helpful in understanding both how to solve a problem with FOAM and what the different integrators are useful for. Future releases will include many more solvers.

Sorry for the super lame demo graphics- another early article will detail the implementation of a decent FOAM renderer.

Re-zipped up the source after making numerous bug fixes and adding a new example. This fixes the clumsy, problematic way I was removing elements across the board. Check it out here.

I also put the documentation up on my server and will keep it up-to-date.

• googlecode project hosting with svn repository
• zipped source

I’ve been working on FOAM off and on in my spare time for the past few months. Recently I’ve been working on solving constraints on a global level- which is no easy task. I decided to work laterally for a bit and get a release out.

FOAM is primarily intended as a resource for developers interested in simulating physics. It has a carefully thought out OOP structure and modular design. A savvy developer should have no problem extending and repurposing FOAM to his own ends. The Foam datatype is in fact not a physics engine but an interface for simulating physics. It offers a simple means to create, control and run a simulation- it purposefully keeps the more nitty gritty, behind-the-scenes operations shielded from the casual developer. A physics engine is simply part of its composition.

I think the best way to use FOAM as a resource for creating physics simulations is to look through the source. I did a decent job commenting and use fairly descriptive variable names. When I set out working on simulating rigid body physics, I didn’t find one comprehensive resource- an article would be helpful here, a tutorial there, but then the matching source would be a comment-less single-letter-variable mess. Also- not coming from a mathematical background, the often obtuse-sounding descriptions in articles can be daunting. In practice, seemingly complex equations can translate to very simple algorithms in the world of computers. I hope that FOAM bridges some of these gaps.

Over the next few weeks, I’ll be posting tutorials and examples- Please let me know of anything that interests or confuses you and what you think in general!

simple demos included in source:

• Perpetual Fall
• Zero G Toy Chest

I updated my tutorial on dynamic circle collision detection in Flash to demonstrate a more physically realistic reaction. Reaction is given only to illustrate how to implement the results of the detection scheme, but is physically viable.

I started work on a 2D rigid body physics engine that I’ll be releasing as open-source. There are a few AS3 physics engines out there- Fisix, APE, the unreleased Motor… and while I’m sure this will fill a few holes, I’m most interested in it being a resource for developers wanting to simulate physics, as I’ve had trouble finding a comprehensive source of the hows and whys.

Here’s a demo. Use left and right arrow keys to spin the wheel. You can grab the movable objects and toss them, too.

As it gets further into development I’ll write some tutorials to explain some of the engine’s inner-workings etc.

Current Features:

• Rigid body simulation
• Arbitrary convex polygons
• Circles
• Cubic Bezier curves
• Lines
• Constraints
• Springs
• Bungees
• Easily swappable numerical integrators
• RK4
• Euler
• Midpoint
• Separation Axis Theorem based collision detection
• Modular force generation

Next up is splitting collision detection into coarse and fine phases via spatial partitioning.

Oh- and I’m calling it FOAM.

(Updated 09.30.07)

collision detection demo click and drag to move a circle or its velocity vector
source

INTRODUCTION
Determining whether 2 circles are colliding is one of the simplest types of collision detection. You simply evaluate the sum of the circles’ radii and the distance between them. If that sum is greater than the distance, the circles in question are coincident. This call is made at each time step in a simulation.

This works well in most cases and is most often optimized to evaluate the square of the sum of the radii and the square of the distance between circles- which offers the same response without an expensive square root call.

THE PROBLEM
But what if our circles are moving really fast? At the beginning of the interval they may not be colliding and likewise at the end; but it’s quite probable that at some point within that interval, the circles do collide.

THE SOLUTION
We can evaluate as many as 2 possible collisions (given constant velocity)- their initial collision and the moment at which they are adjacent again (if they passed through each other without effecting each other’s movement). There aren’t necessarily 2 solutions; and the second solution does little good for us here anyway. So we’re left wanting to derive:
(A) if the circles will collide in a given interval and if so,
(B) when they first collide in relation to that interval.

This solution is presented in Real-Time Collision Detection by Christer Ericson. He does a much better job of explaining the following derivation than I can, but here we go anyway!

DERIVATION
Because we’re presuming constant velocity over the interval 0 ≤ t ≤ 1, the position of a circle is equal to:

Where X is its current position and v, velocity. So at the beginning of the interval, the circle’s position is equal to X + v * 0 = X, and X + some portion of v for any t value greater than 0. What we need is a function of t that evaluates the distance between the 2 circles according to this parameterization:

Where

and 

When f is equal to the sum of the 2 circles’ radii, the circles are coincident:

Where r is respective radius. Evaluating distance means we need to compute a square root, which as mentioned above is undesirable. We can square both sides of the equation instead:

The square of the sum of vectors expands into dot products just like a scalar polynomial expands into multiplications:

Moving all terms to the left side, and naming each constant based on its relationship with t in agreement with the equation’s quadratic nature, we have:

Where

and

Which is a quadratic equation where b is known to be a multiple of 2. Its roots are:

EFFICIENCY CONSIDERATIONS
Every calculation requires processing and ActionScript is certainly not the exception! We want to organize our detection method such that nothing is defined/calculated before it is necessary. If we can eliminate the possibility of collision before defining a variable- we should! It’s a common mistake to solve without consideration of the algorithm’s real job- to offer an answer to our collision query as quickly as possible. And that’s the driving force here.

IMPLEMENTATION
This example utilizes CircleParticles with public properties representing radius, position and velocity. The Collision type takes 2 of these CircleParticles at instantiation. Its method detect evaluates collision over the interval implied by these CircleParticle’s velocities.

Given collision, it sets a time variable (t) to the point at which in that interval collision first occurs. A basic understanding of vector arithmetic and quadratic equations should be sufficient to wholly understand this function:

/**
 * Detects whether 2 {@link CircleParticle} instances will collide within an implicit timestep of 1
 *
 
 * 		Given collision, {@link #t} is set to the percentage of the interval at which collision occurs.
 * 		ie: 0 if collision occurs at the beginning of the interval and 1.0 at the end of the interval
 *
 * 		The event of collision is given by the quadratic equation ax^2 + 2bx + c = 0 offers the resolution:
 * 		x = ( -b - Math.sqrt( b^2 - ac ) ) / a, the smaller of 2 possible roots
 *
 * 		Note that we return from the algorithm before declaring/defining any unnecessary variables
 * 
 
 *
 * @return {@code true} if collision is detected, {@code false} otherwise
 * */
public function detect() : Boolean
{
 
	//difference vector of the 2 circles' positions
	var cDiff:Vector = c2.position.getDifference( c1.position );
 
	/**
	 * {@code c} the squared distance between circles minus the squared sum of their radii
	 * this offers a means to check for intersection at the start of the interval
	 * without making an expensive square root call
	 *
	 * ie: if the sum of the circles' radii squared > the distance between them squared,
	 * they are overlapping
	 * */
	var c:Number = cDiff.dot( cDiff ) - ( c1.radius + c2.radius ) * ( c1.radius + c2.radius );
	if( c < 0 )
	{
		//initial overlap condition- return with time 0
		t = 0;
		return true;
	}
 
	//difference between circles' velocities
	var vDiff:Vector = c2.velocity.getDifference( c1.velocity );
 
	var a:Number = vDiff.dot( vDiff );
	if( a < Collision.OMEGA ) return false; //circles not moving relative each other
 
	var b:Number = vDiff.dot( cDiff );
	if( b >= 0 ) return false; //circles moving apart
 
	var d:Number = b * b - a * c;
	if( d < 0 ) return false; //no solution to the quadratic equation- circles don't intersect
 
	//evaluate the time of collision
	t = ( -b - Math.sqrt( d ) ) / a;
	return t <= 1;
 
}

NOTE: While there is not a single trigonometric function call in this detection method, some of its assertions are based on the law of cosines.

REACTION
So, assuming collision at exactly halfway through the interval, our Collision’s t property is set to 0.5. We can now move on to react to the collision.

The idea or our Collision’s react method is to take our knowledge of the moment of collision within the interval and apply that to our CircleParticles. We first move the circles along their velocity vectors scaled by the collision interval. We then calculate the transfer of momentum along the axis of collision and move the circles from their position at collision along their new velocity vectors by the remainder of the interval:

/**
 * Reacts to a collision between {@link #c1} and {@link #c2} at specified time {@link #t}
 *
 *
 
 * 		We project the ghost particles along the velocity vector of c1 and c2 by
 *		t * circle.velocity.
 *         (Description is outdated- needs updating)
 * 		Then we normalize the collision vector and amplify by the magnitude of the
 * 		opposing circle's velocity to obtain each ghost's velocity. Their velocities
 * 		are then multiplied by 1 - t to give us each reaction vector in reference to the
 * 		remaining span of the interval.
 * 
 
 *
 * @param g1 ghost particle representative of {@link c1} at moment of collision
 * @param g2 ghost particle representative of {@link c2} at moment of collision
 * */
public function react( g1:CircleParticle, g2:CircleParticle ) : void
{
 
	//assign the position at percentage t along each velocity vector to each ghost
	g1.position = c1.position.getSum( c1.velocity.getProduct( t ) );
	g2.position = c2.position.getSum( c2.velocity.getProduct( t ) );
 
	//difference vector between circles at time of collision
	var cDiff:Vector = g2.position.getDifference( g1.position );
	//normalized
	var collisionNormal:Vector = cDiff.getUnit();
 
	//relative velocity of circles
	var relativeVelocity:Vector = c1.velocity.getDifference( c2.velocity );
	//coefficient of restitution range of 0 - 1.0
	//0 is completely inelastic (lump of clay) and 1 is purely elastic (superball)
	var restitution:Number = Math.sqrt( c1.elasticity * c2.elasticity );
 
	//the impulse created at collision-
	//this is given by the following equation (j is used to represent impulse)
	//j = -( 1 + e )vAB . n / n . n( 1 / mA + 1 / mB )
	//where e is the coefficient of restitution, vAB is the relative velocity
	//n is the collision normal and mA/mB are the respective masses of the different circles
	// "."s represent the vector dot product
	var impulse:Number = collisionNormal.dot( relativeVelocity.getProduct( -( 1 + restitution ) ) );
	impulse /= collisionNormal.dot( collisionNormal.getProduct( 1 / c1.mass + 1 / c2.mass ) );
 
	//scale each impulse along the collision normal by the circles' respective masses
	var reactionA:Vector = collisionNormal.getProduct( impulse / c1.mass );
	var reactionB:Vector = collisionNormal.getProduct( -impulse / c2.mass );
 
	//Multiply this by the the remainder of the interval (1 - t)
	g1.velocity = c1.velocity.getSum( reactionA ).getProduct( 1 - t );
	g2.velocity = c2.velocity.getSum( reactionB ).getProduct( 1 - t );
 
}