Physics engines for dummies

Can you code this?
Submit to StumbleUpon

Hello and welcome back to my blog!

This time i’m going to talk about the basic components that make up a physics engine and how to put them together; this tutorial is aimed at programmers who have a basic grasp of maths and geometry but would like to step into the world of simulation.

It is my hope that if, at the beginning of this article, you are able to code the physics behind the 1972 game Pong, by the end of the article you will be equally happy writing your own constraints to use in your own physics solver!

Can you code this?

I’ve always found the title of those books of from the ‘…for dummies’ series reassuring; after all, if a dummy can learn this stuff you should stand a good chance, right?

...for dummies?

Hence the title of this article.

Plan of action

Ok, so i’m going to cover a few of the things you might want in a physics engine:

  • Rigid bodies
  • Collisions
  • Resting contact
  • Constraints (Joints)

Simulation

Bear with me, this is going to start out really basically, but hopefully later on it will become clear why.

Starting out with particles, which we assume will have a position P and a velocity V. Each frame we advance the position of P by adding the velocity V on to it. This is called integration.

You can use whatever units you like for you simulation; the usual choice is one unit/meter which is what i’m using here – the screen is two metres by two metres (in our world), so our velocities are specified in metres/second. To make this work, we must know how many seconds there are per frame, then we can do P += V*dt in our integration to advance the particles properly, where dt is the number of seconds per frame.

We can simulate multiple particles by storing an array of positions and velocities, and looping over integrating each one. To prevent our particles going off screen we would like to collide them against the screen’s edges.

To bounce the particles off the edges we simply detect collisions against the edges and reverse the velocity in the offending axis, thus:

for all i
  if (P[i].x < -1 && V[i].x < 0)
  {
    V[i].x = -V[i].x;
  }
  if (P[i].y > 1 && V[i].y > 0)
  {
    V[i].y = -V[i].y;
  };
  if (P[i].x > 1 && V[i].x > 0)
  {
    V[i].x = -V[i].x;
  }
  if (P[i].y < -1 && V[i].y < 0)
  {
    V[i].y = -V[i].y;
  }

The first condition of each case checks for intersection of the particle in space and the second checks that the particle is actually heading towards the edge. The second condition is important otherwise we would do the collision response again next frame erronusly should the particles move more than one pixel per frame. This second condition carries though all of impulse based rigid body simulation and is what separates inequality contraints (e.g contacts) and equality constraints (most joints).

This is the kind of example code you will no doubt have written a long time ago when you first started getting interested in simulation. I’m including it here because I think there is a natural progression from this easy example to more complex physics code.

Whats going on in this simple simulation?

Well, without knowing it, we have assumed a physical matieral type for our particles including the coefficient of restituition and followed Newton’s principle of conservation of momentum. We have also chosen that the world have infinite mass such that it doesn’t feel a collision response when the particles impact on it. By choosing to reflect particles on impact we are also preseving their momentum (even though we have ignored their mass in our calculations), indicating that we have chosen a coefficient of restitution of 1, i.e. perfectly elastic – like a super-ball. Additionally, we have chosen an impulse/velocity model for the collision response, rather than a force/acceleration model – we did this simply by changing the velocity instantaneously, rather than over a period of time.

What we have here is actually a highly optimised special case physics simulation. Optimised how, you ask? Well, let me explain:

If we were to code the above “properly”, not taking any short cuts, we would have to assume the following. Our environment is defined by its boundaries which are four axis aligned planes (actually, they’re lines because we are in 2d). Each one has a normal (which points inwards) and a distance to the origin. They look like this:

Planes[4] =
{
   (1,  0, 1),
   (0, -1, 1),
   (-1, 0, 1),
   (0,  1, 1)
}

So, now we have to detect our collisions as we did before, but we can’t take any shortcuts this time. In order to detect if our particles have penetrated the planes we must perform a dot product and add the plane’s distance to origin:

for all particles i
{
  for all planes j
  {
    distance = P[i].x*Planes[j].a + P[i].y*Planes[j].b + Planes[j].c;
 
    if (distance < 0)
    {
      // collision responce
    }
  }
}

What this code is doing is finding, by projection, how much of the vector from the plane to the particle is in the direction of the plane normal. Because our plane’s normals are unit length, this is a measure of the closest distance from the particle to the plane. Obviously, if this distance is less than 0, our particle has penetrated the plane and we need to take action to perform the collision response.

Now, if we take a closer look at that distance check above, including each plane’s coefficents:

plane0dist = P[i].x*1  + P[i].y*0  + 1;
plane1dist = P[i].x*0  + P[i].y*-1 + 1;
plane2dist = P[i].x*-1 + P[i].y*0  + 1;
plane3dist = P[i].x*0  + P[i].y*1  + 1;

If we simplify that down a bit:

plane0dist =   P[i].x + 1;
plane1dist =  -P[i].y + 1;
plane2dist =  -P[i].x + 1;
plane3dist =   P[i].y + 1;

Re-arranging slightly we get these plane tests conditions:

if (P[i].x  < -1)
if (-P[i].y < -1) = if (P[i].y > 1)
if (-P[i].x < -1) = if (P[i].x > 1)
if (P[i].y  < -1)

Which are exactly the same as those which we used in our first simple example. The second condition for each plane must also be satisfied so we only do the collision response once per collision. This can be done by performing a dot product between the particle’s velocity and the normal for each plane;

if (P[i].x < -1 && V[i]•N[i] < 0)

We call this the normal velocity since its the particle’s velocity in the direction of the normal. If this results in a value less than zero then the particle is moving towards the plane and we allow the collision. Once the collision is resolved, the normal velocity will be >= 0 depending on the coefficient of restitution. I could perform the same analysis on the second condition as i did the first to prove that this general solution is the same as the first specific version, but i will leave that as an excercise for the reader.

This shows our two approaches are physically the same.

Ok, so what do we do now for the collision response?

We need a solution which achieves the same result as the original program but is more general. We can do this using the reflection vector from the law of reflection. The reflection vector is calculated thus:

R = V – 2*N*(VN)

Where V is the velocity vector and N the surface normal. We need to do a similar comparison with this operator as we did with the collision checks:

plane0vel x = V.x - 2* 1*(V.x* 1 + V.y* 0) = V.x -2*V.x   =   -V.x
plane0vel y = V.y - 2* 0*(V.x* 1 + V.y* 0) = V.y - 0      =    V.y
plane1vel x = V.x - 2* 0*(V.x* 0 + V.y*-1) = V.x - 0      =    V.x
plane1vel y = V.y - 2*-1*(V.x* 0 + V.y*-1) = V.y - 2*V.y  =   -V.y
plane2vel x = V.x - 2*-1*(V.x*-1 + V.y* 0) = V.x - 2*V.x  =   -V.x
plane2vel y = V.y - 2* 0*(V.x*-1 + V.y* 0)                =    V.y
plane3vel x = V.x - 2* 0*(V.x* 0 + V.y* 1)                =    V.x
plane3vel y = V.y - 2* 1*(V.x* 0 + V.y* 1) = V.y - 2*V.y  =   -V.y

Which, you should be able to see, are the exact same resulting velocities after collision with each plane in the first simple example.

You will have noticed that mass has been ignored in our calculations so far, this is because it just drops right out of the equations when you are dealing with collision against a body of infinite mass (i.e. our simple world).

So what does this new version look like so far in pseudo-code?

for all particles i
{
  for all planes j
  {
    N = {Planes.a, Planes.b};
    distance = P[i]•N + Planes[j].c;
 
    if (distance < 0 && V[i]•N < 0)
    {
      // collision response, reflect particle
      V[i] -= 2*N*N•V[i];
    }
  }
}

Great, so we’ve turned a nice simple example into one that’s far more complicated than it needs to be. Why?

Firstly, to demonstrate that the roots of all our inital assumptions were grounded in physics and secondly, to demonstrate the advantages of the more complex implementation. Because we now have a way to deal with arbitrary 2d planes instead of axis aligned planes, it means we can rotate our world and still have the simulation work correctly:

Move the mouse over the demo to rotate the planes. Note that in this demo, i have done some correction to the points as the planes are rotated to make sure the points still lie within the planes.

Ok, so this is all very well but its not a very real world simulation. There is no gravity and the coefficient of restistution is far too high – in the real world, most everyday things have a near plastic coefficient of restitution, i.e near 0.

So, adding gravity is pretty simple, we just subtract gravity from our particle’s velocities before we do any collision, remembering to account for the fact that gravity is an acceleration:

V[i] += G*dt

In this case G is the vector {0, -9.8} and dt is the number of seconds per frame as before.

We would also like to use a different coefficient of restitution as well because our particles might not be made of super-rubber. Recall the reflection vector equation that we used earlier:

R = V – 2*N*(VN)

We can actually re-write this equation to include the coefficient of restitution:

R = V – (1+e)*N*(VN)

Where e is the coefficient of restitution which varies from 0 (totally plastic) to 1 (totally elastic). You can see that these two equations are equivelent; if we set e to 1, we get the original equation. The particles in this example have coefficients of restitution varying from 0 to 1.

The only thing we can now do to this particle simulation to make it more realistic (without going to a full rigid body simulation with rotation) is to handle masses in our collisions. To accomplish this, we will have to use circles instead of particles, since its quite unlikely that two particles will ever collide what with them having no actual size!

Figure 1

As per Figure 1, two circles A and B intersect if the distance between them is less than the sum of their radii. The collision normal is simply the vector d between the two, but normalised.

In order to deal with the masses of our particles we need to do some thinking regarding the reflection equation that we already have. What we really need is some kind of mass ratio to apply a bias to the equation so that lighter particles get pushed out of the way by heavier particles, but we need it to conserve momentum as well.

ratioa = Mb / (Ma + Mb)
ratiob = Ma / (Ma + Mb)

The above two ratios will do the trick. Here is an example:

Ma = 1.0
Mb = 0.5

ratioa = 0.5 / (1.0 + 0.5) = 1/3
ratiob = 1.0 / (1.0 + 0.5) = 2/3

Obviously 1/3 + 2/3 = 1, so you can see that this conserves momentum correctly, and since body b is half the mass of body a, it should experience twice the reaction upon collision, which it does.

So we can plug these two ratios into our reflection vector equation to obtain two new equations for the reflected velocity of each body after collision:

Va = Va – (1+e)*N*((Vb-Va) N)*(Mb / (Ma+Mb))
Vb = Vb – (1+e)*-N*((Vb-Va)
-N)*(Ma / (Ma+Mb))

Our collision code gives us the normal pointing from body b towards body a, so we must reverse the normal to calculate body b‘s reflected velocity. You will have noticed that there are a lot of duplicated calculations in the two equations above. We would like to simplify this down a bit. What we can do is use the relative velocity of the two objects in question:

Vr = Va - Vb

Now that we have only one velocity to deal with the equation gets simpiler.

I = (1+e)*N*(Vr N)

Va = -I*(Mb / (Ma+Mb))
Vb = +I*(Ma / (Ma+Mb))

What we are doing here is to treat one object as being stationary relative to the other moving object. But there is still more we can do here to simplify these equations:

I = (1+e)*N*(Vr N)*(Ma*Mb)/( Ma+Mb)

Va – = I /  Ma
Vb + = I / Mb

By combining the ratios we can save some more calculation and reduce the final step to a simple divide though by the mass of the objects in question to produce the final velocities for those objects. You can see this works because

(Ma*Mb)/( Ma+Mb) / Ma Mb/( Ma+Mb)

and

(Ma*Mb)/( Ma+Mb) / Mb = Ma / (Ma+Mb)

Finally we can show that (from fractions)

(Ma+Mb)/( Ma*Mb) = 1/Ma + 1/Mb

Which means we can re-write the final linear impulse equation

I = (1+e)*N*(Vr N) / (1/Ma + 1/Mb)

This is handy because we can just store 1/Ma and 1/Mb inside the definition of our circle rigid bodies so we don’t have to re-calculate. It also means that we now have a way to represent infinitly massive objects, such as our world (by storing the inverse mass as 0).

Va – = I * 1/Ma
Vb + = I * 1/Mb

In these demos you can manipulate the circles with the mouse.

Ok, so you probably already noticed we are starting to see problems due to penetration caused by the collision detection system detecting the collision after its already too late. This problem only gets more severe as we add more rigid-bodies, and is compounded by the fact that the system is currently doing nothing to correct the penetration:

To combat this problem, we need to be able to detect collisions before they actually happen and deal with penetration. Luckily i’ve covered this technique already in a previous article which i suggest you read.

Once this technique is applied the problem just goes away; and there are only a few lines of code required. The only caveat is that we have to assume that the coefficient of restitution is never required to be non-0; in the majority of real-world cases this isn’t too much of a limitation, especially given the benefits.

Software engineering

In physics engines, the mathsy part is only half the story; the second half is how you organise your physics engine in terms of classes hierarchy and how you actually code it.

I’m going to present one possible solution here which has worked for me in the past; it may be over-simplified for your requirements but it should give you an idea. The language here is actionscript, but the principles apply to any OO language.

public class RigidBody
{
	protected var m_pos:Vector2;
	protected var m_vel:Vector2;
	protected var m_invMass:Number;
 
        ///
	///
	///
	public function RigidBody( pos:Vector2, vel:Vector2, invMass:Number )
	{
	        m_pos=pos;
		m_vel=vel;
		m_invMass=invMass;
	}
 
        ///
	///
	///
	public function GenerateContact( rb:RigidBody ):Contact
	{
		throw new Error("Not implemented on base class");
	}
 
	///
	///
	///
	public function Integrate( dt:Number ):void
	{
		m_pos=m_pos.Add( m_vel.MulScalar( dt ) );
	}
}

So, i’ve defined a RigidBody base class containing the three parameters that we need so far in our example. There is an integrate function to move foward in time and (what should be a virtual) function called GenerateContact() which will generate a collision normal and distance between any two RigidBodys.

public class Circle extends RigidBody
{
	private var m_radius:Number;
 
        ///
	///
	///
	public function Circle( pos:Vector2, radius:Number, invMass:Number )
	{
		m_radius = radius;
		super(pos, new Vector2(), invMass);
	}
 
        ///
	///
	///
	public override function GenerateContact( rb:RigidBody ) :Contact
	{
		if ( rb is Circle )
		{
			...
		}
                else if (rb is ...)
                {
                        ...
                }
                else
		{
		    throw new Error("unahandled case!");
                }
	}
}

And we have a Circle which derives from RigidBody and extends its functionality by having a radius parameter. Also, we have implemented the GenerateContact() function to return the required information when needed.

By separating functionality out into base class and concrete implementations we will be able to deal with multiple RigidBodys in the same lists, which greatly simplifies the code.

public class Plane extends RigidBody
{
	private var m_n:Vector2;
	private var m_d:Number
 
	public function Plane(n:Vector2, d:Number )
	{
		m_n=n;
		m_d=d;
 
		super(n.MulScalar(-d), new Vector2(), 0);
	}
 
        ///
	///
	///
	public override function GenerateContact( rb:RigidBody ):Contact
	{
		if ( rb is Particle )
		{
			...
		}
		else if ( rb is Circle )
		{
			...
		}
	        else
                {
			throw new Error("unhandled case!");
                }
	}

Above is another concrete implementation, this time representing an infinite plane (our screen edges).

public class Contact
{
	public var m_normal:Vector2;
	public var m_dist:Number;
 
	///
	///
	///
	public function Contact( n:Vector2, dist:Number )
	{
		m_normal = n;
		m_dist = dist;
	}
}

Above is the definition for the Contact class, which represents the information the system will need to resolve one collision.

private function Update( e:GameLoopEvent ):void
{
	const dt:Number = Math.min(e.m_elapsed, 1.0/15.0);
 
	// apply gravity
	for each ( var p:RigidBody in m_rigidBodies )
	{
		if ( p.m_InvMass>0 )
		{
			p.m_vel=p.m_vel.Add( Constants.kGravity.MulScalar( dt ) );
		}
	}
 
	// collide
	for ( var i:int=0; i<0||rbj.m_InvMass>0 )
	{
		const c:Contact=rbi.GenerateContact( rbj );
 
		...
 
		//resolve collision
	}
 
	// integrate
	for each ( p in m_rigidBodies )
	{
		p.Integrate( dt );
	}
}

Above is what our update loop currently looks like. The order here is fairly important;

  1. Firstly, external forces are applied – i.e. gravity.
  2. Then, all the collision detection is done to generate a contact which needs to be resolved for each pair of interacting rigid bodies.
  3. Finally, each RigidBody is integrated forward in time.

The time-step clamping at the beginning is so that as we debug the code, we don’t end up with such a massive time-step the next frame that everything explodes!

Constraints

Ok, so we have a pile of spheres up and running and our physics engine is starting to take shape quite nicely, what about constraints?

Before i talk about this, i feel i should introduce some concepts needed to understand the problem.

Constraint

Loosely, this is something which enforces a condition between two rigid bodies. So, when there is a collision we create a constraint which enforces the rule that the colliding rigid body may not move through the object it collides with but also that it do this by only pushing. We already derrived the code for this above (not taking into account rotation yet):

I = (1+e)*N*(Vr N) / (1/Ma + 1/Mb)

All the constraints i’m going to work with are impulse velocity level constraints. Forces and accellerations will not come into it.

Inequality constraint

This is a constraint which only acts in one direction. So a collision constraint is an inequality constraint because it is only allowed to push and never pull. If it were allowed to pull, it would stick the rigid body to the object it collided with. The way we enforce this inequality is simply the check i mentioned at the beginning:

if (V•N < 0)
{
   // handle constaint
}

Little Big Planet contained a lot of constraints of this type; winch, piston, string are a few examples – they only allowed motion in one direction within the limits of the constraint. So a piston, for example was a length constraint with lower and upper bounds on the allowed length.

Equality constraint

This constraint type is allowed to pull and push equally. An example is a point to point constraint, where two rigid bodies are connected by a single point – one point on each body is forced to lie in the same place. With a hinge constraint the whole hinge axis is constrained to lie in the same place on both rigid bodies. A rod constraint would be a a length constraint where the length is never allowed to change.

Designing a constraint

The way to approach how to design a particular constraint is to think about what restricting effect it has on the two rigid-bodies its connected between.

For example, a distance constraint (a rod in LBP) simply prevents the two end points of the rod getting any closer than (or any further away from) the specified distance. Putting the end points exactly at the centre of mass of each rigid body, as we are about to describe keeps simple circular rigid bodies at a constant distance from each other, while still allowing both of them to continue to move around in 2d space.

Figure 2

In Figure 2, bodies A and B must not get any closer than distance d, but are still allowed to move around independently otherwise – this resolves as each body being able to orbit the other as they fly through space.

Impulse/velocity level

As mentioned previously, the simulator i’m describing is an impulse/velocity level simulator, rather than a force/acceleration level one. This means that all the constraints must only be solved by applying impulses to change velocities of bodies.

What this means when designing a constraint is that once you’ve worked out how it behaves on a position/distance level (i.e. on paper, with a diagram), you need to then think about how that translates into velocities.

So, to recap: our distance constraint stops the end points getting any closer or further away from each other. In velocity terms this is describing a one dimensional constraint; the relative velocity of the bodies is constrained to be zero in one particular axis. That axis is defined by vector between the end points.

Figure 3

In Figure 3, A and B have starting velocities which violate the length constraint that we’ve placed between them (from Figure 2). The first step in solving a constraint is to figure out where that velocity actually is. In this case you can see its the projection of Av and Bv on to the length constraint axis.

Once we have the velocities we want to remove, we convert into a single, relative velocity (as we did before with the contact constraint) and we already derived the maths required to solve one dimensional constraints (from the contact constraint):

I = (1+e)*N*(Vr • N) / (1/Ma + 1/Mb)

Lets re-write that so it becomes less of a mouth-full:

I = RelativeVelocityMagnitudeToRemove / (1/Ma + 1/Mb)

In addition to the velocity part, there will also be some positional fix-up to do, since with more than one constraint acting in the system simultaneously it will often fail to be resolved with one one solver iteration. In this simple example we just add an artificial correction to the velocity (current length – desired length) / time-step.

Above is the result – you can manipulate the circles with the mouse as before, although one is fixed in place to demonstrate the constraint chain.

Software engineering

So, lets cover the additions to the code.

public class Constraint
{
	protected var m_bodyA:RigidBody;
	protected var m_bodyB:RigidBody;
 
	/// 
	///
	/// 
	public function Constraint( bodyA:RigidBody, bodyB:RigidBody )
	{
		m_bodyA = bodyA;
		m_bodyB = bodyB;
 
		Assert( m_bodyA.m_InvMass>0||m_bodyB.m_InvMass>0, "Constraint between two infinite mass bodies not allowed" );
	}
 
	/// 
	///
	/// 
	public function ApplyImpulse( I:Vector2 ):void
	{
		m_bodyA.m_vel = m_bodyA.m_vel.Add( I.MulScalar(m_bodyA.m_InvMass) );
		m_bodyB.m_vel = m_bodyB.m_vel.Sub( I.MulScalar(m_bodyB.m_InvMass) );
	}
 
	/// 
	///
	/// 
	public function Solve( dt:Number ):void
	{
		throw new Error("base class doesn't implement!");
	}
}

Above you can see we’ve added another base class; this one defines the basic part of a constraint but doesn’t actually define any concrete details – as before, we leave this to the derived classes. It does however provide a framework for us – each constraint must implement the Solve() function which does the meat of the work. Additionally, there is an ApplyImpulse() function defined to save us some work in the concrete implementations.

public class DistanceConstraint extends Constraint
{
	private var m_distance:Number;
 
	/// 
	///
	/// 
	public function DistanceConstraint( bodyA:RigidBody, bodyB:RigidBody, distance:Number )
	{
		super(bodyA, bodyB);
 
		m_distance = distance;
	}
 
	/// 
	///
	/// 
	public override function Solve( dt:Number ):void
	{
		// get some information that we need
		const axis:Vector2 = m_bodyB.m_Pos.Sub(m_bodyA.m_Pos);
		const currentDistance:Number = axis.m_Len;
		const unitAxis:Vector2 = axis.MulScalar(1/currentDistance);
 
		// calculate relative velocity in the axis, we want to remove this
		const relVel:Number = m_bodyB.m_vel.Sub(m_bodyA.m_vel).Dot(unitAxis);
 
		const relDist:Number = currentDistance-m_distance;
 
		// calculate impulse to solve
		const remove:Number = relVel+relDist/dt;
		const impulse:Number = remove / (m_bodyA.m_InvMass + m_bodyB.m_InvMass);
 
		// generate impulse vector
		const I:Vector2 = unitAxis.MulScalar(impulse);
 
		// apply
		ApplyImpulse(I);
	}
}

Above is the implementation for the actual distance constraint, it does the work as described above.

private var m_joints:Vector.<Constraint>;

In the solver we define a list of constraints – they are different to contacts in that they’re more permanent, so they get their own list.

const dt:Number = Math.min(e.m_elapsed, 1.0/15.0);
 
// apply gravity
for each ( var p:RigidBody in m_rigidBodies )
{
	...
}
 
// collide
for ( var i:int=0; i<m_rigidBodies.length-1; i++ )
{
	const rbi:RigidBody=m_rigidBodies[i];
	for ( var j:int=i+1; j<m_rigidBodies.length; j++ )
	{
		const rbj:RigidBody=m_rigidBodies[j];
 
		if ( rbi.m_InvMass>0||rbj.m_InvMass>0 )
		{
			...
		}
	}
}
 
// solve constraints
for each(var jt:Constraint in m_joints)
{
	jt.Solve(dt);
} 
 
// integrate
for each ( p in m_rigidBodies )
{
	p.Integrate( dt );
}

The solver loop now looks like the above – we have gained a little loop for solving the constraints – i’ve chosen to do that after solving for contacts which makes the joints ‘stronger’ than the contacts (in that because they’re resolved last, they will be more correct than the contacts which are done first), you can do it the other way around if you like.

Conclusion

This article has uncovered the physics lurking behind the most simple of 2d games (pong), shown that the principles therein are grounded in physics and expanded on the technique used in that game step by step all the way up to the point where designing a constraint in a physics engine can be seen as a simple extension of colliding two objects.

I hope this article has whet your appetite for physics simulation while also being easy to understand. I have run out of time with this article, but depending on its popularity i can write a lot more on the subject. Register your interest by posting in the comments! :)

Source code

As ever, if you enjoyed reading this article and would like to see more, please consider buying the source code which accompanies it; this will allow an indie developer like me to pay the rent and buy more food to eat! :)

You will receive all the framework code that i’ve laid out in this article, including the code behind the demos on this page which should serve as a good base to expand upon. They are written in actionscript 3.0 but the techniques are applicable to all OO languages.

After purchasing, you will be redirected to a page where you can download the source immediately.

USD 14.99

Subscribers can access the source here

Until next time, have fun!

Cheers, Paul.

Submit to StumbleUpon

About Paul Firth

A games industry veteran of ten years, seven of which spent at Sony Computer Entertainment Europe, he has had key technical roles on triple-A titles like the Bafta Award Winning Little Big Planet (PSP), 24: The Game (PS2), special effects work on Heavenly Sword (PS3), some in-show graphics on the BBC’s version of Robot Wars, the TV show, as well as a few more obscure projects.   Now joint CEO of Wildbunny, he is able to give himself hiccups simply by coughing.   1NobNQ88UoYePFi5QbibuRJP3TtLhh65Jp
This entry was posted in AS3, Beginner, Collision Detection, Geometry, Graphics, Physics, Technical and tagged , , , , , , , . Bookmark the permalink.

185 Responses to Physics engines for dummies

  1. Another superb article, Paul!

    Love it! :D

    (Seriously, people – These aren’t the ramblings of some deranged hobo, Paul is one clever muffin indeed!)

    Keep posting this stuff Paul. You will have enough for a book soon enough!

    • nip says:

      You sure? I read:
      “this will allow an indie developer like me to pay the rent and buy more food to eat! ”
      as
      “will code for food”.

      • val says:

        Any indie developer able to live off his own business has to be clever. Working for food so he can follow his passion is a common start, even after experience on big companies.

  2. questor says:

    thanks for the article, it’s a really great introduction!

    I’ve also read your previous article and really like little big planet, but I have the feeling that the physics in lbp are a bit “fluffy”. I can’t really describe it, but I have the feeling that something is a really tiny bit off. May it be because of the speculative contacts and the missing energy?

    • Paul Firth says:

      Hi questor,

      Thanks!

      Regarding the fluffy feeling physics in LBP – this might be down to the character controller (the thing which drives sackboy’s motion); this came in for a bit of criticism after the game was launched. In the PSP version we changed that a bit to make sackboy feel more solid… Could that be what you’re experiencing? :)

      Cheers, Paul.

  3. raigan says:

    Awesome!

    I’ve got a question about the position correction; you’re doing “error = relVel + relDist/dt” which seems like a Baumgarte factor of 1… surely this is way too much and would result in lots of noticeable added velocity?

    • Paul Firth says:

      Hi raigan,

      Yes, this is probably true – it worked ok in the sample on the page though, which was only running 1 iteration. It may need some tweaking when deployed in production code, of course :)

      In LBP PSP we were using split impulses for this, but i wanted to keep the tutorial as simple as possible…

      Cheers, Paul.

  4. Pingback: fizix » Turtlehead Blag

  5. The SkyLife says:

    I happen to be a physicist and have done some programming in my day… excellent write from both fields of study. good job.

  6. kevin says:

    Cheers man, you just made my multi media class so much simpler.

  7. Ralph Dratman says:

    Thanks for this article. It’s very lucid and non-threatening.

    I used to write code like this. Yours are good principles to follow. But at the time, I wanted to make a rigid body, and could never figure out how to do it. Real rigid bodies consist of about 10^23 little masses coupled to each other with springs. Since I never quite grasped the whole Hamiltonian thing, I just stuck with the masses and springs I understood. Unfortunately, that is the long, slow way to create a quasi-rigid body.

    Can the methods you use here be extended to make something 2-D or 3-D and rigid?

    • Paul Firth says:

      Hi Ralph,

      Yes, these methods extend into 3d very easily :)

      The masses and springs model is very like what Jacobsen talked about in his early article about rigid bodies composed of length constrants:

      http://www.gamasutra.com/resource_guide/20030121/jacobson_01.shtml

      Cheers, Paul.

      • Ralph Dratman says:

        That looks like a good article by Jacobsen. I used to build all kinds of little computable objects from masses and springs. You can get interesting effects by varying individual masses and spring constants, and of course gravity is easy to add.

        I used to worry about approximating the square root. Eventually I decided it wasn’t worth it to over-optimize. Especially with today’s processors there is almost no speed penalty for using floating point and square roots. Writing code for fixed point is not a great way to enjoy life.

        Spring models are outstandingly stable, and objects can in principle be made very realistic (with some work), but it’s a lot of computation.

        The simplest spring has zero equilibrium length, but it’s easy to give them non-zero length. I never did very much with collisions. With lots of little masses and springs, fully general collision tests would take far too much time. I will study Jacobsen’s method for making that practical.

        Thanks again for your article.

        • Paul Firth says:

          Hi Ralph,

          I was very excited when i first read that Jacobsen article and indeed the way he pushes the constraint solver onto the geometry that make up the interacting rigid bodies is an amazing idea :) However, when it comes to modelling rigid bodies with this technique i ran into so many problems trying to maintain the rigidness, and with objects flipping inside out on fast collisions that i eventually abandoned the technique in favour of traditional rigid body dynamics but adopting the relaxation technique for solving the constraints :)

          Cheers, Paul.

  8. Pingback: Physics engines for dummies | Paul's blog@Wildbunny

  9. Pingback: Physics engines for dummies | Paul’s blog@Wildbunny : Popular Links : eConsultant

  10. Pingback: Physics Engines for Dummies « Wobbits

  11. Pingback: Physics engines for dummies | Pauls blog@Wildbunny

  12. Pingback: futilefiles » Blog Archive » Physics Engines for Dummies

  13. Pingback: Learning A Foreign Language: Hebrew For Dummies

  14. Pingback: My daily readings 04/08/2011 « Strange Kite

  15. Vikram says:

    Wow! Amazin tuts.. Thanks

  16. Kevin Reid says:

    Nice article! I might pass it on to the next person asking about how to do simulations.

    One suggestion: In “The collision normal is simply the vector d between the two, but normalised” it might be worth explaining what “normalised” means and that it’s different from “normal vector”.

  17. Latika patil says:

    I dont like PHYSICS…!!!! :) :) :)

  18. LoveIt says:

    Hey, lovely write-up!
    I’ve got 2 weeks off, and writing my own one of these is just what I need to keep me entertained!

    One type I noticed: “You can use whatever units you like for you simulation; the usual choice is one unit/meter ” – I think you meant unit/second or unit/frame.

  19. ScreamingMage says:

    Outstanding.

    Thanks you for putting this up!

  20. Pingback: JavaScript Magazine Blog for JSMag » Blog Archive » News roundup: iOS viewport fixes, Mobile Boilerplate, CommunityJS, Ender.js

  21. Pingback: Links for 2011-04-08 « Micha? Piaskowski

  22. m_eiman says:

    Consider adding flattr.com to your sources of revenue, I’d have clicked that button right away!

  23. Pingback: Robert McGhee » April 9th

  24. Florian says:

    Your physics simulations have some issues with oscillations and constraint enforcements, and because of that you made them really “soft”. There is a solution to solve these problems without making everything feel packed in wool, which is verlet integration. I’ve written two articles about that (including examples and code as well). I also discuss a hybrid approach to verlet integration with impulse preservation and stable at-rest acceleration.

    • Paul Firth says:

      Hi Florian,

      The reason things look ‘soft’ as you say, is that the solver is only running 1 iteration :) Nothing to do with the integration method…. Thanks for the links, though :)

      Cheers, Paul.

    • Mike says:

      Hi there Florian, your sim looks very good. How many iterations are you using for the “bridge builder” style steel girder sim?

  25. Wolter says:

    You lost me at the “Planes[4] =” part.

    What do all the 1, 0, and -1 values mean? Why are there 12 numbers there?

    “Each one has a normal (which points inwards)” — What does this mean?

    • Paul Firth says:

      Hi Wolter,

      The planes array is describing the planes which make up the edges of the screen, the values {x,y, d} in each plane are {direction x, direction y, distance to plane}.

      The first two components describe a unit length direction vector (also called a normal, sometimes), the third represents the distance from the plane to the origin :) The normals point inwards so that everything inside the screen is a positive distance to each plane…

      Hope that makes sense,

      Cheers, Paul.

      • Jake Boxer says:

        Thanks for this clarification. There is one point that I’m still a little unclear on. Maybe you can confirm my understanding.

        At first, I thought that the origin point of each of these normals was the middle of the screen. If that were the case, I think Planes[0] would represent the right plane, Planes[1] would represent the bottom plane, and so on (clockwise).

        However, after re-reading “the normals point inwards”, it seems to me that the middle of the screen is actually the endpoint of each normal. If this is the case, then Planes[0] would represent the left plane, Planes[1] would represent the top plane, and so on (still clockwise).

        Is my second interpretation correct? It makes sense, but it’s a bit different from other vector representations I’ve seen, so I want to make sure I’m understanding correctly.

        • Paul Firth says:

          Hi Jake,

          The planes are

          Planes[] =
          {
          left,
          bottom,
          right,
          top
          }

          The coordinate system in flash has down the screen = positive, which is why 0,-1 points up (and is hence the bottom plane).

          Try not to think about the end-points of the normals – they’re just a unit length direction vector, so they would all exist exactly in the middle of the screen (0,0), except that we add the distance parameter which ‘pushes’ them out along their direction vector by the distance to origin. e.g:

          point on plane = direction vector * -distance

          For the left plane (1, 0, 1):

          point on plane = (1,0) * -distance = -1,0

          For the right plane (-1, 0, 1):

          point on plane = (-1,0) * -distance = 1,0

          We use negative distance because the planes are all pointing inwards, so to approach them we must move against their direction vectors :)

          Cheers, Paul.

  26. Pingback: Learn A Language: Copyediting and Proofreading For Dummies

  27. Pingback: Eric Blue’s Blog » Weekly Lifestream for April 10th

  28. Derrick says:

    Paul,
    I’ve been enjoying your tutorials so far.

    I’ve noticed in some of the demos that the circles go outside of the bounding boxes (or overlap, etc). When the change in x,y is greater than 1 pixel, this can happen. What can be done to fix this ?

    Thanks!

  29. Pingback: BertaS' Blog» Blogarchiv » Physik Engine für Dummies

  30. Volgogradetzzz says:

    Hello, Paul. Thanks for tutorials. Sadly that you didn’t mention friction. Or maybe there will be separate tutorial?

  31. Sinan says:

    i have done some correction to the points as the planes are rotated to make sure the points still lie within the planes`

    Can you give us information about corrections you’ve done to make sure the points still lie within the planes?

    I’m trying to code this examples in Pygame(Python), and having difficult times to hold points between planes.

    • Paul Firth says:
      if ( dist<0 )
      {
      	if ( relNv<0 )
      	{
      		p.m_vel = p.m_vel.Sub( n.MulScalar( 2*relNv ) );
      	}
      	p.m_correct = n.MulScalar( -dist );
      }
      

      p.m_correct is added on in the integration phase, then zeroed :) This corrects for penetartion...

      • Sinan says:

        Can you please more explanatory? What is relNv, p.m_vel, Sub(), m_correct, etc.

        Actually it would be better if you just tell me the way you correct the positions of points :) .

        Thanks.

        • Paul Firth says:

          Sorry, I find the penetration distance for each point, multiply by the contact normal, store that value and then add it on to the position in the integration phase… :)

  32. Pingback: Physics engines for dummies | News Ninja

  33. Pingback: Collision detection for dummies | Paul's blog@Wildbunny

  34. Ben says:

    This is a very informative and well written article.

    Thank you for spending the time to educate us!

  35. Pingback: Blog J.Schweiss | Physics engine

  36. Brandon says:

    This is a great article, it really made the math easy enough to implement. However, I although I have it implemented up to the coefficient of restitution, I don’t really understand how the math relates to the macroscopic representation (i.e. the simulation). I’ve never taken linear algebra so although I can manipulate the math I don’t comprehend how it relates. If you are looking for more material for articles in the future I would be extremely grateful if you would cover that topic.

    Here is my (rudimentary) implementation of this in Javascript / HTML 5:
    http://cgi.cs.arizona.edu/~bmf/physical/

    Click to make bunches of balls that float around. Press r to make them random colors. Press g to turn on gravity and watch them bounce.

    • Paul Firth says:

      Hi Brandon,

      When you ask how the maths relates to the simulation, which part are you talking about? Is it just a software engineering question, or more of a mathsy one? :)

      You demo look good, though!

      Cheers, Paul.

  37. Ivan says:

    Hi Paul. Thanks for the article – really inspiring. :)

    One thing I would like to ask you to clarify is why add distance to the origin in “In order to detect if our particles have penetrated the planes we must perform a dot product and add the plane’s distance to origin …” instead of subtracting from it?

    Given that the normals are inverted as in originating at the origin and pointing toward the planes; the dot product, assuming normals are of unit length, gives component of particle vector onto one of the normals. Why not subtract it from plane’s distance to the origin? In this case if a particle is closer to plane A than to plane B, distance to plane A is 1 – comp_An_ p and distance to plane B is 1 + comp_Bn_ p because the component is negative now. Here An, An are normals to planes A and B respectively, p is a particle’s position vector.

    Since your demos work and no-one complained about this part, I’m probably just rushing through the article and it is misunderstanding on my side (hope I didn’t make a fool of myself here :) ), but still, could you (or anyone) explain why we are adding but not subtracting the component.

    Thanks a bunch!

    Cheers.

    • Paul Firth says:

      Hi Ivan,

      The reason the distance is added on to the dot of the particle and the plane is that the planes all have their normals pointing towards the origin (so the world appears to be an inside out box, with all normals pointing inwards):

      Consider the left plane, at distance d=1 from the origin with the normal 1,0. It actually sits at -1 in the x axis, since its on the left.

      Particle P at -0.5, 0.5 is 0.5 units distance from this plane.

      P . N = -0.5
      P . N + 1 = -0.5 + 1 = 0.5
      :)

      Hope that helps,

      Cheers, Paul.

  38. Pingback: Real-Time Rendering · Seven Things for May 4th, 2011

  39. David says:

    This is a great article, it really made the math easy enough to implement. (with some work). Keep posting this kind of material . And i can’t wait to buy your book in the future.

  40. Stephen says:

    Great article Paul, really takes alot of the confusion out of the whole collision & response ideas that I’ve been stuck on.

    One question I have though, (and I don’t know if this would be compatible with this style of engine) would be to have bodies with bounding shapes made up of squares & circles, forming a custom collision shape… how would this be done? or is it something completely different?

    • Paul Firth says:

      Hi Stephen,

      Glad you enjoyed the article :)

      If you want a custom collision shape it should be quite easy to add one – for example you could add another class which extends RigidBody called CompoundObject which would actually contain a list of other shapes – it might be worth separating the direct relationship of Shape->RigidBody so that you can have a list of shapes corresponding to only one RigidBody (which is what you want).

      Then you can just query the individual shapes within the CompoundObject to generate your contact information as before :)

      Cheers, Paul.

  41. Henry says:

    Great article Paul. I have a few friends who are physicists and I think they would really enjoy this. I will have to share this with them.

    Thanks,
    Henry

    cantilever parasol</a

  42. Pingback: Programando la física del movimiento (con ejemplos) | CyberHades

  43. Pingback: iPhone Programming Book and a physics engine - iPhone Dev SDK Forum

  44. Pingback: How to make Angry Birds – part 2 | Paul's blog@Wildbunny

  45. walter says:

    Hi!
    First of all, thanks for the great article. I’m new to game development and this article helped a lot.

    Now for my newbie question, you said:

    “GenerateContact() will generate a collision normal and distance between any two RigidBodys”

    Then you have this code:

    // collide
    for ( var i:int=0; i0 )
    {
    const c:Contact=rbi.GenerateContact( rbj );

    ...

    //resolve collision
    }

    In your 1st example collision of a particle in a plane, what is the collision normal and the distance?

    Is the collision normal equal to the vector perpendicular to the plane? (which is always pointing to the origin (0, 0))?

    Is the distance equal to 1? Is it because you declared the planes to be:

    (1, 0, 1),
    (0, -1, 1),
    (-1, 0, 1),
    (0, 1, 1)

    Again, thanks a lot!

  46. walter says:

    Hi Paul,

    Another question regarding the “GenerateContact” method, say you have a class Plane and a class Circle, in its GenerateContact method, will it have the same logic? Please see code below:


    public class Plane extends RigidBody {
    ...

    public override function GenerateContact( rb:RigidBody ):Contact {
    if ( rb is Plane )
    ...
    else if ( rb is Circle )
    //is the logic here the same as in the class Circle?
    else if ( rb is Particle )
    ...
    else
    throw new Error("unhandled case!");
    }
    }

    public class Circle extends RigidBody {
    ...

    public override function GenerateContact( rb:RigidBody ):Contact {
    if ( rb is Plane )
    //is the logic here the same as in the class Plane?
    else if ( rb is Circle )
    ...
    else if ( rb is Particle )
    ...
    else
    throw new Error("unhandled case!");
    }
    }

    Thanks again!

    • Paul Firth says:

      Hi walter,

      The logic is the same, but the objects ra and rb are swapped around… This is what I do in the 2nd case:

      if (rb is Plane)
      {
      var contact:Contact = rb->GenerateContact(this);
      // flip normal, and A and B references around
      contact.flip();
      }

      So I avoid having to write the same code twice by simply swapping the inputs around :)

      Cheers, Paul.

  47. Liy says:

    Hi Paul,
    Thank you for this great article!

    I have a question about your equations for the reflected velocity of two objects after collision:

    Va = Va – (1+e)*N*(Va • N)*(Mb / (Ma+Mb))
    Vb = Vb – (1+e)*-N*(Vb • -N)*(Ma / (Ma+Mb))

    If Va=[1,0], Vb=[-1,0], N=[1,0], Ma=1, Mb=1, e=1.
    I imagine the two objects’ velocity should be [-1,0] and [1,0] after collision, but your equation always returns 0. Am I doing something wrong?

    Could you also explain a bit more about the equations below. Why there is a minus sign in Va and opposite in Vb?
    I = (1+e)*N*(Vr • N)
    Va = -I*(Mb / (Ma+Mb))
    Vb = +I*(Ma / (Ma+Mb))

    Excuse my poor English…
    Thanks!

    • Paul Firth says:

      Hi Liy,

      I think you found a mistake! :)

      I think it should be:

      Va = Va – (1+e)*N*((Vb-Va) • N)*(Mb / (Ma+Mb))
      Vb = Vb – (1+e)*-N*((Vb-Va) • -N)*(Ma / (Ma+Mb))
      :)

      But I would use the equations just below in the article; they are less complex to follow…

      There is a – + here

      Va = -I*(Mb / (Ma+Mb))
      Vb = +I*(Ma / (Ma+Mb))

      because we’re applying an impulse equal and oppositely to each object :)

      Cheers, Paul.

      • Liy says:

        I’ll use your final impulse equations.
        But I’m trying to understand your two equations:

        Va = Va – (1+e)*N*((Vb-Va) • N)*(Mb / (Ma+Mb))
        Vb = Vb – (1+e)*-N*((Vb-Va) • -N)*(Ma / (Ma+Mb))

        Still quite confusing…

        Can you show me how you plug the mass ratio into original R = V – (1+e)*N*(V•N) and formed up the two equations? Any hints will be much appreciated!
        Thanks!

        • Paul Firth says:

          The mass ratios are Mb / (Ma+Mb) and Ma / (Ma+Mb), and simply scale the response for each body based on the relative masses. That’s why the reflection vector equation is multiplied by the individual ratio for each. I guess they don’t ‘plug in’ exactly because there is no mass-ratio term in the reflection vector equation :)

          Cheers, Paul.

  48. Tim says:

    Hi Paul,

    Great articles. I have a couple of questions.

    1. In modeling physics using constraints, is the reflection of a ball off of a wall or floor considered a ‘constraint’? I wasn’t sure if I simply used the reflection math at the collision or if I should turn it into a constraint and solve it with the other constraints.

    2. In the formula for reflection, there is “(1+e)”. I understand that e=coefficient of restitution of one of the objects. Does that mean “1″ is the COR of the other? So, if I wanted a floor that deadened the reflection some, would I do “…(0.5+e)…”.

    Thanks!

    • Paul Firth says:

      Hi Tim,

      1) Reflection of a ball off a wall would be considered a contact, which is like a transient constraint which only exists for one frame :)
      2) (1+e) is from the original reflection vector equation (well, it was 2 originally)… You need to combine the two CORs into one – how you combine the COR of two objects during a collision is a bit hand-wavey and vague… I *think* in previous engines, I’ve just picked the minimum of the two :)

      Cheers, Paul.

  49. Pingback: Quora

  50. Tim says:

    Hi Paul,

    I’ve been working with the impulse based physics engine, but I’m a little stuck. What I noticed with the “contact” is that things seem to come to a stop instead of bounce (when not rotating).

    I wanted to make a ball bounce, so I created a new type of contact called “bounce”. The challenge I had was that if I simply applied the reflection equation at the point where the motion bounding-box collides, it doesn’t look like the ball hits the ground.

    So, in my “bounce” contact, I modified the velocity based on your original “contact” to get the ball to the ground this frame and then I take the remaining reflection velocity and apply it the next frame.

    95% of the time, it works perfectly. However, every once in a while my “delayed velocity” adds a little extra jump to the ball. I started looking at the original “contact” again and noticed that the ball does bounce, but only if it’s rotating.

    I was hoping you could enlighten me a bit on the best way to get a ball to bounce (whether it’s rotating or not) while still having it look like it impacts the ground.

    I hope this all makes sense. :)

    Thanks,
    Tim

    • Paul Firth says:

      Hi Tim,

      Which source code are you using? :) The physics engine for dummies code didn’t handle rotation at all…

      Having a non 0 coefficient of restitution is hard if you want to use speculative contacts. I’d recommend against trying to use them together because they’re not really suited…

      At a push I would suggest you try and do what you have already attempted – though I’ve never done this myself!

      Sorry I can’t be more help,

      Cheers, Paul.

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>

WP-SpamFree by Pole Position Marketing