Archive

Archive for the ‘Physics’ Category

Constrained Dynamics 1: Particle Constrained to Curve

December 31st, 2009

This is the first in a series on simulating constraints; in this case a particle is constrained to move along a curve in planar space. To follow along, you’ll need to know a little physics, calculus and some algebra, and hopefully you’ve implemented a simple particle system before!

I’ve chosen a parabolic curve. We’re actually going to solve the problem in 2 different ways, so we’ll be seeing 2 separate representations of the curve, the first being the implicit form:

implicit form of parabola

Which looks like this for x {-200…200}:

parabola graphed from -200 to 200

The first method we’ll use for solving the constraint is described for a circle in this paper by Andrew Witkin. There are some convenient aspects of his example, and so it may not be entirely apparent how to constrain the particle to our parabola. I’ll follow Witkin’s format to show how his formula works in this context.

First, you write a positional constraint which is just our curve equation:

positional constraint

The reason the curve is the constraint is because it implies that when the function’s value at (x,y) is zero, then [x,y] is located on the parabola. Next, we define the velocity constraint, which is the derivative of the positional constraint with respect to time:

velocity constraint

If the velocity constraint is satisfied, then we’ll always end up with a valid position. Finally, we find the derivative again with respect to time. This is the acceleration constraint function:

acceleration constraint

As was the case with velocity → position, if the acceleration constraint is satisfied, so will the velocity. This means that, provided we’re starting with a valid position and velocity, we only ever need to enforce the acceleration constraint on the system.

Let’s do a quick review of Newton’s 2nd law:

Newton's 2nd law

Which states that force is equal to mass times acceleration. We’ll be making use of the fact that F represents the sum of all forces acting on the particle. So:

Newton's component 2nd law

Solving for acceleration gives:

Newton's component 2nd law solved for acceleration

The subscripts ext and c stand for external and constraint respectively. It’s worth pointing out for clarity’s sake that when a term is bold it represents a vector, whereas italicized variables are scalars.

Now we can replace the appearances of x and y’s 2nd derivatives, breaking F into its component parts to allow the constraint force to appear explicitly:

Equation 1


Equation 1

The 2 unknowns now appear, but only in the 1 equation, so we can’t solve yet. If we consider a particle attached to a frictionless wire, it should be apparent that the wire, which is the source of the constraint force, can do no work and so it can not add nor remove energy from the system. The kinetic energy of the system is defined as:

kinetic

To gain anything useful out of this, we have to consider what direction a force can push without changing T. From examining the derivative of this function with respect to time:

kinetic

it should be apparent that the constraint force has an important relationship with the velocity vector, specifically that their dot product should equal zero (or else T would change as a direct result of the constraint force). We can use this knowledge to write the constraint force components into a 2nd equation:

Equation 2


Equation 2

This makes sense not just in terms of energy conservation, but also based on the constraint function chain described above. We already know our particle has a legal velocity, and so we don’t want to add any force that would ruin the circle of trust–the constraint force should only pull the external forces onto the acceleration curve.

I first solved for the x constraint component in terms of y from equation 2:

Equation 3


Equation 3

Solving for y’ in the velocity constraint and substituting for x’ in equation 3 reduces the variable count:

Constraint force x

Plugging the right-hand side in for Fcx in equation 1 and solving for Fcy yields:

Constraint force y

We could explicitly solve for Fcx too, but that’s a non-trivial number of multiplications and a divide. Luckily, we don’t have to because this is for a computational simulation; after Fcy is known, the solution to Fcx is given in equation 3.

Plugged into F=ma, we have the equations of motion for a particle constrained to the parabola x2 + 100y = 0.



Parametrized Form

Now we’ll solve the problem a different way, using a reformulation of Newton’s equations called Lagrangian mechanics. Lagrangian mechanics involves finding the number of degrees of freedom (DOF) in a system, then expressing that system in terms of generalized coordinates, one for each DOF. In our example we have a single DOF, even though our particle exists in a 2-dimensional space. First, we phrase our parabola parametrically:

Equation 4


Equation 4

Here, q is our lone generalized coordinate. Noting that x is a vector, the transformation from the generalized coordinate to our particle’s position is given in the parametric equation (4):

Parametrized position

The Lagrangian equations of motion for a particle experiencing only conservative forces, as is the case with a particle constrained to a wire, is:

Equation 5


Equation 5

David Eberly gives a full derivation of equation 5 from Newton’s 2nd law in his excellent, if often daunting, Game Physics. F, again, represents all externally applied conservative forces and T is the kinetic energy of the system, which is a function of q and q’. Before we can define T for this system, we have to derive dx/dq and x’ because kinetic energy is dependent upon velocity.

velocity

and

dxdq

Now we can define the kinetic energy of the system as a function of our generalized coordinate because x’ can be expressed entirely in terms of q and its first derivative:

kinetic energy

Having determined the kinetic energy function, we can define the rest of the terms in equation 5. We’ll see that its first term includes q’s second derivative, which is the last piece needed to advance the simulation.

derivatives

Lastly, F, in this case, is just the gravitational force, which we define as the vector mg[0,1]:

generalized force

Next we plug all of this back into equation 5, which then defines the Lagrangian equations of motion for a particle constrained to the parametrized parabola (q,-q2/100):

Equation 6


Equation 6

Solving for acceleration yields:

Equation 6

Now that we’ve solved for acceleration, we can advance our 2-dimensional simulation by numerically integrating our 1-dimensional differential equations. Below is a demo, showing both methods running side-by-side. If you watch long enough, you’ll see the implicit form (left) eventually leaves the parabola as error accumulates. In contrast, the parametric form will always be on the parabola, but energy will dissipate. These errors are associated with numerical integration inaccuracies and floating point imprecision.

Roll over to run simulation. Left is the implicit form, right is the parametric


Get Adobe Flash player



Considerations

Each method has advantages and disadvantages. Both offer a relatively easy-to-follow format. Explicitly solving for a constraint force as we did with the implicit curve is more intuitive to me, but the parametric form is more computationally efficient, and the algebra isn’t as messy. A downside to the parametric form is that it isn’t a generic solution–adding new constraints to the system means solving a totally new set of equations of motion.

Below are the separate implementations that each run in a generic dynamics harness. You can see that there is more involved computationally in the implicit form:

Implicit implementation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
var particle:Particle = Particle( p[ 0 ] );
particle.exertForce( G ); //apply gravity
 
var m:Number = particle.mass;
 
var X:Vector3D = particle.position;
var V:Vector3D = particle.velocity;
 
var Fext:Vector3D = particle.force; //generically grab all applied forces
 
var Fc:Vector3D = new Vector3D();
 
//solve for the constraint force
Fc.y = ( -50 * m * V.x * V.x - 100 * X.x * Fext.x - 2500 * Fext.y ) / ( X.x * X.x + 2500 );
Fc.x = ( Fc.y * X.x ) / 50;
 
particle.exertForce( Fc ); //exert the constraint force on the particle
Parametric implementation
1
2
3
4
var q:Number = x[ 0 ];
var qPrime:Number = v[ 0 ];
 
a[ 0 ] = ( -q * qPrime * qPrime - 50 * G * q ) / ( 2500 + q * q );

As mentioned above, both of these are prone to computational error. In the demo, I’m using an RK4 integrator. Using a Euler integration scheme results in visible error within the first second of simulation. Although RK4 does a good job of reducing error, it isn’t perfect. Ideally, an error term is introduced to keep drift in check.

drew ActionScript 3.0, Mathematics, Physics

BallRacer!

April 21st, 2009

This is my first multi-player game! You race a tiny silhouetted hamster in a fashionable, plastic globe! I recommend playing it over at OMGPOP so that you’re not constrained to the embedded dimensions here.

drew ActionScript 3.0, Network, Physics

Dealing with Synchronicity in Multiplayer Simulations

April 5th, 2009

At OMGPOP, I’m working on a physics-based, real-time multiplayer game. Keeping a user-run simulation synchronized is a difficult (read: impossible) challenge. I’m fairly happy with the result, so I thought I’d share the idea and logic behind the networked simulation synchronization pipeline.

hammy

The model on the end user application consists of the current player and a list of all other players in the match. Each player is a subclass of a MovableElement type which defines relevant properties:

MovableElement
{
    //position
    public var x:Number;
    public var y:Number;
 
    //velocity
    public var vx:Number;
    public var vy:Number;
 
    //angular velocity
    public var av:Number;
}

Each user runs a copy of the simulation locally. Because frame rate can vary so drastically from machine to machine, it’s necessary to advance the simulation by a dynamic time-step. This introduced a problem in the collision detection scheme, as a large enough time-step would allow the character to tunnel through the floor without intersection being detected. As per the previous post, collision detection occurs between a circle and a line segment. To guarantee tunneling does not occur (without a true continuous solution, and thus a rewrite), I had to ensure that the time-step was small enough that any change in position was less than the circle’s radius:

max Δt = (r - λ) / ||v||

Where r is the circle’s radius, λ an error term to ensure we’re working with less than the radius, and v is the linear velocity of the circle in question.

Here is the main simulation loop, including solving for the maximum allowed time-step, without any opponent input:

//solve for dt
var dt:Number = synchronizedTime - lastTime;
 
//ActionScript does not throw a divide by zero error- max will be positive infinity if speed = 0
var max:Number = ( currentPlayer.radius - 1e-8 ) / currentPlayer.velocity.magnitude;
 
physicsEngine.step( dt, max );//ensures we never take a sub-step larger than max
 
lastTime += dt;

Now, when a user changes input (through key presses) and thus applies forces to their character, a notification must be sent to everyone with the current physical state of the sending player’s character; alongside a list of keys being pressed and a timestamp. Upon receipt of the character state, the corresponding player is updated on the receiver’s end, and the simulation is rewound to the time the sending player changed inputs, and then re-simulated over the latent time for just that player, according to the new set of inputs (which map to forces applied remotely to all opponents). This looks something like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
private function onUserUpdatedRemotely( event:CharacterUpdateEvent ) : void
{
 
    var s:PlayerState = event.playerState;
    //find the local instance of the remote player who has updated
    var updatedOpponent:MovableElement = playerMap[ event.characterID ];
 
    //copy the state
    updatedOpponent.x = s.x;
    updatedOpponent.y = s.y;
    updatedOpponent.vx = s.vx; 
    updatedOpponent.vy = s.vy;
    updatedOpponent.av = s.av;
 
    //update forces being applied to opponent
    remoteKeyMapper.updateForces( updatedOpponent, event.keyList );
 
    //solve for change in time and advance just this opponent
    var dt:Number = synchronizedTime - event.updateTime;
    physicsEngine.stepSingleElement( updatedOpponent, dt );
 
}

And that’s pretty much it- very simple. All the hard work is writing a clean model and tight logic so that it doesn’t become a heaping pile of updates here and events there, etc.

drew ActionScript 3.0, Network, Physics

SLRBS demo

February 11th, 2009

demo image

Here’s a demo of the auto-gen physics engine (SLRBS- Super Lightweight Rigid Body Simulator) I’ve been working on which allows you to select images to use as rigid bodies. First, the image is analyzed, ignoring transparent pixels, to determine its convex hull. Then mass/center of mass/inertia tensor are solved for. Finally, the rigid body is created and dropped into the world.

I’m sure that in playing with it, you’ll notice all sorts of wonky behavior… but I think it’s fun anyway. Let me know if there’s a class or type of shape/image that tends to be particularly misbehaved!

demo
image pack

drew ActionScript 3.0, Flash Player 10, Geometry, Mathematics, Physics

Fluid Dynamics

February 6th, 2009

This was something I did a while ago and never got around to posting about. It’s a simulation of an incompressible fluid- when you render the density and apply temperature by adding an upward force, it looks like smoke!

It’s fairly computationally intensive. The state of the system has to be solved simultaneously. But because it uses an iterative solver to approximate a solution, it’s easy to tune for performance to realism.

Click and drag to knock the smoke around, release to add density (more smoke).

My implementation is based on this paper by Joe Stam. The leap from the Navier Stokes equations to this compact form is pretty radical. Mr. Stam’s a pretty smart guy.

drew ActionScript 3.0, Mathematics, Physics

SLRBS (Super Lightweight Rigid Body Simulator)

February 1st, 2009

For the past few months I’ve been working on a general purpose geometry library for representing your interactive elements in a more precise geometrical way. There are some pretty cool features like collision detection management, auto-generation of geometry, and an easy way to bind the model to your DisplayObject hierarchy (or not!). 

It will be released under the MIT license alongside SLRBS- a physics engine that uses the library. The absolute last thing the Flash world needs is another physics engine, but this has some really nice features for beginner developers, prototypers etc. And it does a great job of showing how to use the geometry library.

One of the most frequent questions or requests I get/got with FOAM was something along the lines of “how do I make my MovieClip a rigid body?” It can be hard from my position to know how to tell someone who probably uses the Flash authoring tool a proper way to handle the difference between the model and view. So with SLRBS, I set out to address this issue.

With SLRBS, you have a typical API to create physics simulation, but it also has some streamlined auto-creation methods that make setting it up a breeze. To ensure this, I even downloaded a copy of Flash CS4 to play with it!!! AAAAHH!!! wtf happened to Flash!? 

So, all you do is lay out your MovieClips on the stage, and call a single static method and voila! physics simulation! There’s some pretty cool stuff going on that determine how to handle overlapping DisplayObjects (each demo does it differently), and it’s just a single parameter.

Here’s a screenshot of the stage in a demo where a car is being created:

Car Demo

Because 2 of the tires overlap the car, they get pinned to the car at their center of mass.

 

This is the only ActionScript in the .fla:

1
2
3
4
5
6
7
8
9
import org.generalrelativity.slrbs.SLRBSimulator;
 
var simulator:SLRBSimulator = SLRBSimulator.createSLRBSimulatorFromContainer( this );
addEventListener( Event.ENTER_FRAME, onEnterFrame );
 
function onEnterFrame( event:Event ) : void
{
   simulator.step( 1 / 60 );
}

 

Here’s a second demo’s stage in the Flash authoring tool- this time a ragdoll. The overlap this time gets pinned at the center of overlap instead of the center of mass as in the car demo- this works well for creating ragdoll joints.

Ragdoll Demo

drew ActionScript 3.0, Flash Player 10, Mathematics, Physics

Quick Maths 1: Projectiles and Targeting

May 26th, 2008

This is the first in a series of quick tutorials on helpful maths, physics and algorithms, and their ActionScript implementations. Today we’re going to learn how to shoot a projectile out of the sky. Sweet.

DEMO

First, let’s define the algorithm’s structure. We’re going to pass the time to collision. The algorithm will return the muzzle velocity that our projectile requires to hit a moving target out of the air after this specified amount of time.

function solveMuzzleVelocity( timeToImpact:Number ) : Vector

Next, we have to define the equations of motion:

x1 is a particle’s position after t seconds. x0 is the original position and the dots represent the first and second positional derivatives (velocity and acceleration, respectively). This equation allows us to solve for position after any specified time (given constant acceleration). Referring to our algorithm’s structure, we know t (timeToImpact) to be defined. This means that we can determine where the target will be when we hit it, which gives us a point of impact:

//point of impact
var poi:Vector = target.position.plus( target.velocity.times( timeToImpact ) );
poi.plusEquals( new Vector( 0, Projectile.G * timeToImpact * timeToImpact * 0.5 ) );

It’s assumed that the only accelerative force is gravity (Projectile.G here), that we know the starting position of the target and our projectile, and the current velocity of the target.
Plugging this point of impact in for our projectile’s equation of motion (as x1), we eliminate all unknowns except the muzzle velocity. Through basic algebraic manipulation, we can isolate and solve for this unknown:

Finally, divide each side of the equation by t:

The algorithm will return this value. Here’s my implementation, in full:

1
2
3
4
5
6
7
8
9
10
11
12
13
private function solveMuzzleVelocity( timeToImpact:Number ) : Vector
{
 
    //point of impact
    var poi:Vector = target.position.plus( target.velocity.times( timeToImpact ) );
    poi.plusEquals( new Vector( 0, Projectile.G * timeToImpact * timeToImpact * 0.5 ) );
 
    var diff:Vector = poi.minus( projectile.position );
    diff.y -= Projectile.G * timeToImpact * timeToImpact * 0.5;
 
    return diff.dividedBy( timeToImpact );
 
}

Super simple! Now you can integrate dynamically and knock some shit out of the sky.

Here’s a verbose implementation that doesn’t use my Vector class and isn’t dependent on my application:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
private function solveMuzzleVelocity( t:Number, xp:Point, xt:Point, vt:Point, pixelToMeterRatio:Number = 40 ) : Point
{
 
    const g:Number = 9.81;
    var vDotTerm:Number = g * pixelToMeterRatio * t * t * 0.5;
 
    var poiX:Number = xt.x + vt.x * t;
    var poiY:Number = xt.y + vt.y * t + vDotTerm;
 
    var diffX:Number = poiX - xp.x;
    var diffY:Number = poiY - xp.y - vDotTerm;
 
    return new Point( diffX / t, diffY / t );
 
}

t is time to impact, xp holds the x and y coordinates of the projectile’s starting position, xt is the target’s starting position, vt is the target’s starting velocity and pixelToMeterRatio defines how many pixels represent one meter.

drew ActionScript 3.0, Mathematics, Physics

Rendering FOAM

November 27th, 2007

I’m sorry that this took so long- especially those of you who’ve written me and patiently waited for my help with rendering… I simply don’t have as much time as I would like. Anyway- this will offer a solution for rendering your FOAM physics simulation with:


Here’s a very silly demo- Mouse down anywhere to switch from the DisplayObject renderer to the Bitmap renderer.

demo
source

Note that these new renderers are meant as example only- I do think in most cases the developer will need a more specific setup than I could hope to generically offer. So consider this my implementation and at best a nudge in the right direction.

Hopefully both of these will be fairly self-explanatory, but I’ll briefly go over the DisplayObject renderer and what my thinking was when addressing this issue.

I wanted a means to render a DisplayObject (base class of Sprite, MovieClip etc.) that’s the visual representation of a RigidBody being simulated. I knew that I’d like to be able to pass either the DisplayObject’s Class or the instance itself. I also knew that I’d need an easy way to adjust the position of the DisplayObject- so first, I wrote a datatype to hold these properties- DisplayObjectData.

You’ll see that it has a setter for defining its displayObject property. This is what facilitates passing either a DisplayObject Class or a DisplayObject instance. The other properties:

  • offsetX - the amount to offset the DisplayObject from its origin in the x direction
  • offsetY - the amount to offset the DisplayObject from its origin in the y direction
  • autoCenter - Boolean indicating whether to center the contents about the origin- Because of the way a RigidBody rotates about its center of mass, I figured this would be typical.
  • hasBeenDisplayed - a helper Boolean for determining whether to perform the indicated transformations

In Rocket.as, you can see the way I’m embedding my visual assets:

[Embed(source="../../../bin/images/rocket/lunarPod.png")]
private var rocketSpriteClass:Class;
 
[Embed(source="../../../bin/images/rocket/asteroid.png")]
private var asteroidSpriteClass:Class;

If this is unfamiliar with you, simply note that it’s analogous to the classes you’re instantiating as MovieClips in Flash CS3. As you’ll see in a moment, it’s not necessary to pass the Class, we can pass the DisplayObject instance as well.

Here’s how I’m passing the class to FOAM to tie it to to a RigidBody. Note that the instance of DisplayObjectData simply defines Foam.addElement’s renderData argument.

foam.addElement( rocket, true, true, new DisplayObjectData( rocketSpriteClass, 0, -5 ) );

Compare this with how I’m creating the asteroid (I wanted to scale the asteroid, or else I’d have created it in the same way as the rocket- it worked out well for the purpose of demonstration):

var asteroidSprite:DisplayObject = new asteroidSpriteClass();
asteroidSprite.scaleX = asteroidSprite.scaleY = size / 64;
var asteroid:Circle = new Circle( sx, sy, size, size * 4, svx, svy, 0.5, 0.25, 0, -0.1 + Math.random() * 0.2 );
foam.addElement( asteroid, true, true, new DisplayObjectData( asteroidSprite ) );

There’s one last important thing to note in the Rocket example class. Because I don’t want to use the default FOAM renderer, I have to tell it which renderer to use:

foam.setRenderer( new DisplayObjectFoamRenderer() );

Now let’s move on to DisplayObjectFoamRenderer. Note that I’ve chosen to extend SimpleFoamRenderer so that I don’t have to rewrite drawing routines for objects I don’t have graphics for (such as the mouse dragger bungee). You’ll see my use of the DisplayObjectData properties in the addRenderable method.

So- ultimately all of this boils down to 3 lines of code which line up our DisplayObject with our RigidBody:

displayObject.x = ISimulatable( renderable.element ).x;
displayObject.y = ISimulatable( renderable.element ).y;
if( renderable.element is IBody ) displayObject.rotation = IBody( renderable.element ).q * 180 / Math.PI;

It places the DisplayObject at the x and y coordinates of the RigidBody’s center of mass. Because rotation is given in radians, we have to convert to an angle which is what DisplayObject.rotation expects.

And that’s pretty much it; for the most part a FOAM renderer is bookkeeping. It might suit you far better to just maintain reference to your DisplayObject and update it every frame just as I am in the above 3 lines of code… simply set render to false in Foam.addElement.

drew ActionScript 3.0, FOAM, Mathematics, Physics

Modeling Simple Orbit with FOAM

November 11th, 2007

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.

demo

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:

gravitational 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.

drew ActionScript 3.0, FOAM, Mathematics, Physics

FOAM Update

November 10th, 2007

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.

drew ActionScript 3.0, FOAM, Mathematics, Physics