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.

190 Responses to Physics engines for dummies

  1. Tim says:

    I’m using the Angry Birds Part 2. I’ve probably posted in the wrong spot. Sorry about that.

    The reason for the COR is just to test. I figured if I could get a perfectly elastic collision looking right, then reducing it to something more real would probably be correct as well.

    I’ll keep plugging away at it. If I find a solution, I’ll post it and maybe it will help someone else.

    Thanks.

    • Paul Firth says:

      Working on Little Big Planet PSP we found that the entire game was possible with a COR of 0 – the reason being is that most objects in real life have a plastic (near 0) COR :)

      Cheers, Paul.

  2. Pingback: How a physics engine actually simulates physics? - Programmers Goodies

  3. Pingback: Video Game Collision Detection | Thomas Bergamini

  4. Is there any other way to buy it except pay pal?
    alertpay ?moneybookers?

    • Paul Firth says:

      Hi, right now there is only PayPal, very sorry we tried Google checkout but it didn’t work out for us. Remember, if you have a credit card you can use PayPal to pay without needing to sign up to PayPal itself :)

  5. Pingback: Find point of collision for two sprites - Programmers Goodies

  6. Pingback: Physics dummies | Checkinout

  7. Pingback: help about ipad

  8. tarik kaya says:

    Thank you for this amazingly insightful tutorial.

    But… :)

    I think I found a mistake. You say:

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

    to calculate the new velocities.

    So here is a counter-example: if a ball sits still with zero velocity, there can be no impact possible to get it to some velocity other than zero, according to these equations.
    because: suppose: Va = 0 then: (Va • N) = 0 then: Va=0 no mater what Vb is..

    If I am getting something wrong, I would very much appreciate to be enlightened.

    • Paul Firth says:

      Hi Tarik,

      I think you found a mistake that I thought I’d corrected already :)

      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))

      Cheers, Paul.

  9. Pingback: Your Questions About – Mods For Minecraft | Mod For Minecraft

  10. ryan20fun says:

    nice tutorial, but some of the code blocks have the > and < cymbles translated to HTML tags.

  11. Khada Kuraki says:

    Thanks very much for this wonderful article.

    Khada Kuraki.

  12. Paul says:

    Hi Paul,

    Thank you very much for writing this article, really appreciated :)
    but could you please explain how you find the distance from the line to the origin ? is there some kind of formula? im having a really hard time understanding how to get the distance of a arbitary line segment to the origin.

    Thank you
    Paul

    • Paul Firth says:

      Hi Paul,

      Its a case of projecting the origin onto the line and then getting the distance between the origin and the projected point; this website has the details: http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html

      Cheers, Paul.

    • Aleks says:

      I am digesting this article myself and it took me a while to understand some parts of it. So to make it easier to remember I started a blog in which in I wrote my understanding of how the distance worked. If you are interested you can take a look:

      http://nindzadza.blogspot.com/2011/11/setting-up-proper-screen-boundaries.html

      • TasuLife says:

        Paul, thanks for this. Its easy to see the love you put into it and its very helpful. I know you don’t want to spend the rest of your life maintaining this document, but I also got stuck on this part. I think that Alek’s explanation is much, much easier to digest, but it’s way down here in the comments (I missed it last night). I reccomend adding Alek’s explanation up to your own document when you say:
        “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:”

        Anyway, I absolutely love this document, thank you so much for writing it. It’s a great primer.

        • Paul Firth says:

          Hi there,

          Really glad you enjoyed it :)

          I’ll look up Alek’s explanation and see if I can fit it in…

          Cheers, Paul.

        • Paul Firth says:

          I had a quick look at Alek’s explaination, which I have seen before actually :)

          I can’t really update my explanation to match his, because I think his is technically incorrect – the planes I defined are not just the normals of the planes, they are the actual planes – the first two components are the normal x and y, the final one is the distance to the origin.

          Does that make any sense?

          Cheers, Paul.

          • TasuLife says:

            Yes, that tidbit is helpful. I was mostly confused what the numbers meant (to the layman they look like xyz coordinates). I’ve never seen a plane described like that, but I say blame the parents.

  13. Doug Milford says:

    Hi Paul,

    I’ve enjoyed your tutorials… the first one I tripped upon was with Speculative contacts which I found very interesting.

    I recently purchased the code for Physics for Dummies as I was interested in the section regarding Constraints (the string of connected balls). When I tried to open it in Visual Studio 10, though, nothing came up. Peaking at the code via a Notepad, it looks like it’s a Flash game. Is this correct? Is there a Silverlight version? That’s ok if not as I was only interested in a very specific piece and can most likely translate to Silverlight very easily.

    Thank you very much for providing these demos! The speculative code has already saved me several hours of banging my head against the wall. Well worth the price I paid!

    Doug

    • Paul Firth says:

      Hi Doug,

      Glad you enjoyed the tutorials :)

      I switched to using Flash after that speculative contacts article, because it has a much better install base and I wanted the demos to be visible to as many people as possible…

      I use a visual studio 2008 plug-in called Amethyst which allows flash to work with VS, which is why there is a .sln file :)

      Cheers, Paul.

      • Alex A. says:

        Hi Paul,
        Thank you for awesome tutorials/
        I’m wondering, would you make tuts based on Matthias Muller’s papers like “position-based dynamics” and “Meshless Deformations Based on Shape Matching”?
        Cuz these techniques are very flexible and fits well for 2D, but math explanation is much harder to read than code :(

        • Paul Firth says:

          Hi Alex,

          No problem, glad you enjoyed them :)

          I have to confess, I’ve not read either of those two papers – but thanks for the suggestion, I’ll consider it…

          Cheers, Paul.

          • Alex A. says:

            thank you for answer, Paul.
            I forgot to mention that the main reason why i’m asking about that techniques is because as far as i know developers of original LBP was using those techniques(they was mention about it on siggraph presentation)

            p.s. your articles are awesome, and further direction seems very promising

  14. Pingback: How to make a 2d platform game – part 2 collision detection | Paul's blog@Wildbunny

  15. joaquin says:

    I really like the tutorial, its really hard to find a “step by step” introduction to collision response. In most books they dive into the hardcore and its easy to miss the whole picture, speciality if you want to implement a simple 2D game with some collision.

    I’m so glad Google exists =)

  16. Robert says:

    Hi,
    I just purchased this tutorial. Thank you.
    I would be interested how could one implement bounciness on the balls. Right now if two balls hit each other they kind of slide smoothly without any bounce. If it is more than meets the eye maybe you could make another tutorial.

    Anyway, you have some great content here. Keep it up.

    • Paul Firth says:

      Hi Robert,

      Thanks for the purchase :)

      This is one of the limitations of the technique – it only supports totally plastic collisions (i.e. no bounce). I haven’t looked into trying to add the bounce back in; I think it should be possible. When I get a chance I will have a go and write a post if its successful :)

      Cheers, Paul.

  17. Tomas says:

    Hi!

    I really enjoy your articles, thanks a lot! There is one thing I can’t quite work out though, and internet doesn’t seem to know either:

    if (-P[i].y 1)

    How is it that you can put an assignment operator between two if-statements?

    Cheers!

    • Paul Firth says:

      Hi Tomas,

      Are you referring to this code block?

      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)

      If so, I was just saying that you can simplify the two if conditions in the middle and the result is on the right.

      So if (-P[i].y < -1) can be rewritten as if (P[i].y > 1) :)

      Cheers, Paul.

  18. Tomas says:

    Hey!

    Yeah, the code was cut off somehow! I see what you mean now, I suppose the code block threw me off. Thanks a lot for clearing it up for me! :)

    Cheers!

  19. Raul says:

    Hi. I really liked your tutorial as i’m looking for a physics engine, but a very simple one, this will help me get what I want.
    Altough I haven’t finished reading, I can’t understand about the planes. Because you’re creating a 2D engine but use 3D vectors for the planes? I mean, how do I get a normal for a line, that would be the plane in 2D?
    Hope you understand the question.
    Thanks for the attention.

    • Paul Firth says:

      Hi Raul,

      The plane equation is two parts – the normal to the plane (which is two components) and the distance to the origin (which is the final component).

      Hope that makes sense!

      Cheers, Paul.

      • Raul says:

        oh, I get it.
        So a plane (a, b, c) has the normal that is the ‘vector’ (a, b) and a distance c from the origin.
        Thanks for the info.

  20. Raul says:

    Hello, it’s me again.
    As I finished reading the tutorial, I must say it’s very well written.
    But I have a question about the ‘Software engineering’ part.
    I didn’t get how you handle the collisions, I see you created the class Contact but I don’t get how you use it.
    Is this kind of thing available only in the full code?

    • Paul Firth says:

      Hi Raul,

      A contact just contains the data relevant to solving one collision – so it stores both objects, the contact normal and the distance between objects. These values are then used in the equations you’ve already seen :)

      Cheers, Paul.

  21. Chris says:

    Thank you for the very well written tutorial!

    This is the first one that I’ve seen that actually expexplains how to implement the math in code.
    That was my main problem.

    If I end up making any money from what I’ve learned, I’ll make sure to come back and buy the source code. :)

  22. Alan says:

    Dear Paul,

    When I receive my diploma in software engineering, I’ll remember you.
    Thanks a lot for the lessons!

    Sincerely,
    Alan

    P.S. Your source code is a bit much for a college student’s budget, but I’ll try and contribute one day.

    • Paul Firth says:

      Hi Alan,

      Glad I could help! Don’t worry about the cost, I know what its like being a student! :)

      Cheers, Paul.

      • Alan says:

        May I add, actually, that I’m applying the same principles in this tutorial to 3D. Since everything is in vectors I’m able to do that with only minor modifications (except for the actual collision detection). So, if anyone prefers models over sprites, they will benefit from this tutorial as well. Just putting it out there. :)

  23. Kamil T. says:

    Hi, great article and very useful source code. I’m trying to really *get* everything you’re writing here with all the underlying math and physics.

    And I can’t get past one thing… When you write “above two ratios will do the trick” (regarding conservation of momentum) do you only mean that they are “fine” or that they are precise? Wikipedia has somehow different equations (at the very end of that article):
    http://en.wikipedia.org/wiki/Elastic_collision#Two-_and_three-dimensional

    I know that at that point you already solved part of what those equations are doing (namely you’ve got your angles right by that point) and what I’m asking is if m1/(m1 + m2) has any basis in physics of perfectly elastic rigid bodies or is it just some approximation.

    • Paul Firth says:

      Hi Kamil,

      Yes the two ratios I mention are exact and are fundamental to conservation of linear momentum.

      If you look in that wikkipedia link you can see the same equations hiding along with the other parts specific to angular momentum :)

      Cheers, Paul.

  24. mr5 says:

    Hi there…I would also like to create my own physics engine in 2D.
    I have already done some things. the only thing that bothers me most is: how to check all the coordinates scope of the (x, y) axis , which is the origin + its radius (a circular object) . what I mean is like rotating a needle in the middle through its scope. is it like using some sin/cosine functions there? hope you understand what I mean. like the bouncing balls in your article.

    Too lazy to read.

    thanks!

  25. Rich says:

    Hi, forgive me if this is a dumb question but could you please explain to me what exactly N, the normal, is? I see it in multiple equations. Is it a property of a certain class or is it derived in someway? Just bought your source btw so i’ll look at that when I get a chance.

  26. Joe says:

    Hey Paul ,

    Thank you very much for these tutorials! :)
    I learned a lot about how to design classes and how to apply linear algebra in a practical way. Up to now university was all about theory and your tutorials really help me to close the gap to the practical.

    Right now Im trying to figure out your Contact Design:
    1)where do you store Contacts? Are they stored all in one list somewhere, or has each object a list?
    2)how do i resolve them only knowing “m_normal” “m_dist” later in the update loop.

  27. Tim says:

    Hi Paul,

    This article is just awesome. Thanks.
    I was wondering how much more it would take to implement circles’s rotation (is it a torque in that case) as part of the collision reaction ?
    Would you mind giving some tips about it ? equations, simple way to integrate, so forth.
    Thanks again.

  28. Miles says:

    First, thank you SOSOSOSO much for publishing this blog. Very helpful. I’m not particularly amenable to vector maths, but I’m trying to learn. I’ve got a few basic concepts under my belt now, but I’ve had trouble implementing and understanding your planes representation (newb, I know, but bear with me).

    My little sandbox project here: ( Look in – [MAPlaneManager buildPlanes:] )

    I haven’t been able to manipulate the below items to properly match up with my screen edges. Two sides are off by ~20-30px. Is there some wisdom you can point me to to better help understand this?

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

    How exactly do these vectors determine planes? They’re meant to represent directions and not positions, right?

    Thanks in advance,
    Miles

    • Paul Firth says:

      Those planes are specified in world-space where each plane is one metre from the centre of the world (the last parameter of each vector in the plane).

      Keep all your objects in world space, but when rendering them, convert into screen space. Here is one conversion into screen space:

      render.x = p.m_x*Constants.kWidth*0.5 + Constants.kWidth*0.5;
      render.y = p.m_y*Constants.kHeight*0.5 + Constants.kHeight*0.5;

      Hope that helps!

  29. Benjamin Britten says:

    I take it this is Euler integration, if so you should remake the demos with the Runge-Kutta 4 Method as it performs much better.

    • Paul Firth says:

      Actually not in my experience. In fact, I don’t know of any commercial physics engine which uses anything higher than Euler. The integration is never usually the problem unless you’re trying to simulate a pendulum.

  30. Chris Gregg says:

    Hi there — great tutorial, thanks. I know this page is fairly old, but I thought I’d comment on your integration to determine the effect of gravity on the particles. You use Euler integration: V[i] += G*dt; P[i] += V[i]*dt. This introduces a systematic error into the calculation for position, because the value of the position during an entire time-step is not based on the final result of V[i] during that time-step, rather it is based on the average value: (V[i-1]+V[i])/2. The value of P[i] should then be: P[i]+=(V[i-1]+V[i])/2. See the following for more information: http://lol.zoy.org/blog/2011/12/14/understanding-motion-in-games

    As a concrete example, if you dropped a ball from 10 meters above the ground, after one second (neglecting air resistance) the velocity would be -9.8m/s, but the distance the ball fell would be 4.9m, not -9.8m like an Euler integration would calculate.

    I’m guessing you already know all this and wanted to simplify it for the tutorial, but I thought it was worth mentioning to your readers in case any of them (like me) scratched their heads if they have a bit of physics background and see that the numbers don’t quite work out. Cheers!

    • Chris Gregg says:

      Whoops — that should be P[i]+=(V[i-1]+V[i])/2 * dt.

    • Paul Firth says:

      Hi Chris,

      Yes, Euler is an approximation, but it’s one most games are perfectly happy to live with. I’ve never used anything else in any of the games I’ve worked on because it was simply not necessary :)

      Cheers, Paul.

      • Chris Gregg says:

        Fair enough — if it works well enough for the simulation, then that’s cool. I’m building a physics simulator for teaching physics, so I’m trying to make it as true-to-life as I can (esp. when I tell the students that a ball dropped from rest should fall 4.9m!). I’m not looking forward to trying to wrap my head around simulating friction. :)

        • Paul Firth says:

          Ahhh ok, I understand – I can see why you’d want something more true-to-life than Euler if you’re teaching physics. If you’re teaching simulation physics, I’d actually advice steering well clear of higher order integration methods because they make every single thing much more complicated to understand and to implement.

          Semi-implicit Euler does fairly well in most situations.

  31. Elias says:

    Will you accept bitcoins for the source code?

  32. Hector Vasquez says:

    Hi Paul.

    I know that ‘e’ is the coefficient of restitution, but each body have its own coefficient
    so i want to know which coefficient i will use in the ‘I’ ecuation that is
    I = (1+e)*N*(Vr • N) / (1/Ma + 1/Mb);

    Thanks.

    • Paul Firth says:

      Hi Hector,

      This is a common problem. There is no ‘right’ answer to this question – you’ll have to make a function which picks one single value from the two objects in contact using some kind of heuristic.

      The trouble is trying to understand what happens in cases like a perfectly rubber ball touching a perfectly slippery surface :)

      Cheers, Paul.

      • Hector Vasquez says:

        Ook i will try to understand that, but i have another question.
        Some times when they collide, one of the circles stuck to the other one, i dont know why this is happening

        Thanks

        • Paul Firth says:

          Hi Hector,

          It sounds like you need to make sure relative normal velocity is negative before resolving the impulse. This type of impulse should only ‘push’ never ‘pull’.

          Cheers, Paul.

          • Hector Vasquez says:

            Hi Paul.

            Thanks a lot that worked!! relative normal velocity was positive!

  33. Hector Vasquez says:

    Hi Paul.

    How about friction?
    How can we apply friction to this simulation?

    Thanks.

    • Paul Firth says:

      Hi Hector,

      Yes. Friction is always applied in the tangential direction. So, if you take your collision normal, and perp() it (-y, x) that gives you your tangent vector. Then, project relative velocity onto the tangent to get relative tangential velocity. Once you have this, you can use the same impulse equation to remove some portion of the tangential velocity, and that’s what gives you friction. :)

      Cheers, Paul.

      • Hector Vasquez says:

        Hi Paul.

        So we dont need any cof ?

        • Paul Firth says:

          Coefficient of restitution will be 0, and setting it to 0 will remove some of that equation. The result will be what you need to use for friction – although the amount of tangential velocity you chose to remove could be called the coefficient of friction :)

  34. Felipe Obregon says:

    I’m not sure if you noticed but you’re link to the dot product no longer works, it redirects to the home page of http://www.pfirth.co.uk/

    Excellent tutorial might I add!

  35. Tim says:

    Hey Paul, I’m confused about how the planes are positioned. Are they inside the window like a cross or are they the surrounding the window like a box.

    Thanks

  36. Mathias Chunnoo says:

    Hi paul,

    Thank you so much for writhing this article, it have helped me so much.

    But i wondered if you had any tips on how to make a physic simulation with other objects then circles, like squares and triangles?

  37. Hasan says:

    Hi Paul,
    I am going through a critical problem.You seems a pretty good programmer so I hope you may help.I graduated in CS from University of Windsor in the year 2011. I was pretty good in web programming that time.Due my work permit problem in Canada , no company helped me to extended my stay in Canada. So I started odd jobs.They helped my stay in Canada and now PR is on the way I guess.Recent days, I am planning to be game developer.I started with Unity 3D. I never learned C# before but it is familiar to me as I know Java.but still game programming is far away from web programming. I learned C in school but never learned C++.Do you think it will be a good decision to start learning game programming from the scratch using C++.I tried using Unity it seems nice but when comes critical movement I feel like I am a dumb ass.

    Thanks

    • Paul Firth says:

      Hi Hasan,

      Learning to be a game programmer is a very long process, there are many things to cover. C++ is one of the most difficult languages to grasp, so if you chose C++ you will be starting down an even longer road to learning. On the other hand, all the major games studios currently use C++ so it will stand you in good stead if that’s what you want to do.

      To be honest, though you will probably end up getting paid more to do web development than you would game programming anyway, so that might be a better choice for you in the short term.

      Cheers, Paul.

  38. big_fan_for_phy_engine says:

    Hi Paul,

    Thanks for this article. It’s so great!

    But I am having some difficulties when calculating distance from a point to a plane. If I am not misunderstanding, the plane definition you provide is ‘parsed’ like this: (planeNormal.x, planeNormal.y, distanceFromOrigin). I can see the good point of this definition is to simplify the distance calculation from (a*p.x + b*p.y + c)/sqrt(a*a + b*b) to p*planeNormal + distance . But my problem is: (normal.x, normal.y, distance) actually defines two planes(symmetrical to the origin, as the distance here is a scalar) rather than one.

    So for a given plane and point, (p*planeNormal + distance) gives out different results if the normal is flipped. The wrong result is actullay the distance between the point and the “symmetrical plane”.

    As a result, my collision detector behaves really funny(ball goes through certain planes and collides with ‘invisible’ planes out of canvas). Can you shed some light on this?

    Thanks!

    • Paul Firth says:

      Hi there,

      Even when the plane distance is non unit, the definition only defines one plane because the ‘facing’ direction is encoded into the plane normal. Planes under this definition are in fact one sided, so if you need a double sided plane, you’d need two of these single sided planes back-to-back as it were.

      Hope that helps!

      Cheers, Paul.

      • big_fan_for_phy_engine says:

        Hi Paul,

        Thanks for your fast reply!

        By saying the definition defines two planes, i mean (1, 0, 1) can be seen as a plane facing the direction of X axis and intersecting with X axis at x = 1, but it also defines a plane with the same facing direction but intersects with X axis at x = -1.

        I think if the ‘distance’ in plane definition is a unsigned value, it will be difficult to pick the right plane to perform point to plane distance calculation.

        • Paul Firth says:

          I think you’re describing two different planes in your example:

          (1,0,1) is located at x=-1 because the distance is always in the direction of the plane normal (which is 1,0), so to hit the origin after traveling 1 unit, it has to be positioned at x=-1;
          (-1,0,1) would be the opposite plane, located at x=+1, for a similar but opposite reason.

          Does that make sense?

          Cheers, Paul.

          • big_fan_for_phy_engine says:

            “the distance is always in the direction of the plane normal” makes so much sense to me. My code does not restrict this…

            Thanks for the explanation!

  39. Le Minh Duc says:

    Hi Paul,

    Thanks for a great article.

    I implemented distance constraint, and try a step further to completely lock rotation of two object, but I has no idea. I think I must remove some angular velocities like in the distance constraint, but don’t know how to find them, can you help me?

    • Paul Firth says:

      In the article I don’t consider rotation at all – the objects are connected at their centre of mass, so they’re free to rotate in any direction. If you want to constrain the angular velocity of both objects to be the same, you need to construct an equation to remove all relative angular velocity and then apply an angular impulse separately.

      You’ll need to consider moment of inertia instead of mass in the equation, but it looks quite similar (divide by sum of inverse moments of inertia).

  40. Garan says:

    Hi guys, I just implement the first bounching ball demo in your article by myself. But , in my demo, there is always some ball that goes out of the bounding area when i rotate the world. How did you do about this ? Is there something i missed ?

    P.S This is a awesome article, i will study it very carefully.By the way, english is not my native language, so i am sorry if you do not understand what i mean here.

    • Paul Firth says:

      The code in the article actually corrects for this problem by pushing particles back inside the planes by the amount they are out. It’s quite easy to determine the distance to the plane and therefore correct any negative component in the direction of the plane normal.

      • Garan says:

        Thank you for your reply. I can do it right now. It looks pretty good.

        And in the following toturial, i read about the new reflection vector formula, you said it is like this :

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

        You called the Vb – Va is the relative velocity and it is true in the second formula, but is that right ,when it comes up to the first formula ? Is that vb – va is the relative velocity of the body a to body b?How about va – vb ?

        I just get confus about this .

        • Paul Firth says:

          Both objects must share the same (relative) velocity vector in a collision, only the sign of the normal changes which dictates that the objects should be reflected away from each other post collision.

          • Garan says:

            Sorry to bother you again. But in my demo, when i use the same relative velocity, it doesn’t work. When i use the different one, it works.

            How about these test data?
            Va = (0, 10), Vb = (0, -10)
            and N=(0,-1) e = 1 , Ma = 1, Mb = 1

            In your formula, we can get the result as?
            va = (0, 10) – (1 + 1)* (0, -1) * [(0, -10) - (0,10)] * (0, -1) * 1 / (1 + 1)
            = (0, 10) – (0, -1) * (0, -20) * (0, -1) = (0, 10) – (0, -20)
            = (0, 30)
            vb = (0, -10) – (1 + 1) * (0, 1) * [(0, -10) - (0, 10)]
            * (0, 1) *1 / (1 + 1)
            = (0 , -10) – (0, 1) * (0, -20) * (0, 1)
            = (0, -10) – (0, -20)
            = (0, 10)

            Is this right or am i missed something ???

            By the way, thank for your patient .

  41. Nathan says:

    For the part that is calculating to distance between the particles and planes, is P[i].x the particles x position or its velocity?

  42. Roman Veysman says:

    Hello Paul,

    As I know, constraint must solve contacts created by integrator.
    I have tried to found how you joined generated contacts with constraints in the article, but I have not. Help me to understand, please.

    With best wishes, Roman.

    • Paul Firth says:

      Hi Roman,

      In practice contacts and constraints have to be solved sequentially. So, solve constraints first and then solve for contacts, or visa versa.

      Cheers, Paul.

  43. Hoy says:

    Thanks a lot for this very nice and detailed tutorial, I managed to implement all of it and it’s running fine.

    But there’s one issue that remains: when several balls are stacked on each other, there is a constant jitter. I can see this on the flash example which implements the separation of balls if they interpenetrates each other.
    From what I’ve read, it can involve the time-step of the simulation but I still don’t really understand what can be done to fix the jitter.

    Any thoughts on this ?

  44. Hoy says:

    I took a look at it but i’ve found no references to the jittering issue nor an interactive flash demo with stacked balls like the one this page.
    Do you mean that continuous collision detection helps to deal with the jittering ? I think I’ve read this somewhere long ago…
    But anyway, the jittering I had was on a discrete collision detection i wrote like 2 years ago on a spare-time project. It was only really annoying when I tried to put hundreds of balls in a tiny volume. I moved on to other things, so I probably won’t change it to continuous one. I just randomly landed on your article so I thought I could drop a question about this particular problem I had.

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