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.

drew ActionScript 3.0, FOAM, Mathematics, Physics