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!

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?

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*(V*•*N)*

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*(V*•*N)*

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

*R = V - (1+e)*N*(V*•*N)*

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!

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.

*ratio _{a} = M_{b} / (M_{a} + M_{b})*

*ratio*

_{b}= M_{a}/ (M_{a}+ M_{b})The above two ratios will do the trick. Here is an example:

*M _{a} = 1.0
M_{b} = 0.5*

*ratio _{a} = 0.5 / (1.0 + 0.5) = 1/3
ratio_{b} = 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:

*V _{a} = V_{a} - (1+e)*N*((V_{b}-V_{a}) *•

*N)*(M*

V•

_{b}/ (M_{a}+M_{b}))V

_{b}= V_{b}- (1+e)*-N*((V_{b}-V_{a})*-N)*(M*

_{a}/ (M_{a}+M_{b}))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:

*V _{r} = V_{a }- V_{b}*

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

*I = (1+e)*N*(V _{r} • N)*

*V _{a} = -I*(M_{b} / (M_{a}+M_{b}))*

V_{b} = +I*(M_{a} / (M_{a}+M_{b}))

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*(V _{r} • N)*(M_{a}*M_{b})/(_{ }M_{a}+M_{b})*

*V _{a} - = I / M_{a}*

V_{b} + = I / M_{b}

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

*(M _{a}*M_{b})/(_{ }M_{a}+M_{b}) / M_{a }= M_{b}/(_{ }M_{a}+M_{b})*

and

*(M _{a}*M_{b})/(_{ }M_{a}+M_{b}) / M_{b }=_{ }M_{a} / (M_{a}+M_{b})*

Finally we can show that (from fractions)

*(M _{a}+M_{b})/(_{ }M_{a}*M_{b}) = 1/M_{a }+ 1/M_{b}*

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

*I = (1+e)*N*(V _{r} • N) / (1/M_{a }+ 1/M_{b})*

This is handy because we can just store *1/M _{a }*and

*1/M*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).

_{b}**V _{a} - = I * 1/M_{a}**

**V**

_{b}+ = I * 1/M_{b}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 *RigidBody*s.

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;

- Firstly, external forces are applied - i.e. gravity.
- Then, all the collision detection is done to generate a contact which needs to be resolved for each pair of interacting rigid bodies.
- 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*(V _{r} • N) / (1/M_{a }+ 1/M_{b})*

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.

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.

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 *A _{v}* and

*B*on to the length constraint axis.

_{v}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.

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.

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.

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

Pingback: Video Game Collision Detection | Thomas Bergamini

Is there any other way to buy it except pay pal?

alertpay ?moneybookers?

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

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

Pingback: Physics dummies | Checkinout

Pingback: help about ipad

Hi,

I’m fine with you referencing and linking, but I’m not sure about you using large portions of my material

Cheers, Paul.

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.

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.

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

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

Argh, I thought I fixed that – must be WordPress messing me up; thanks for pointing it out

Thanks very much for this wonderful article.

Khada Kuraki.

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

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.

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

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.

Hi there,

Really glad you enjoyed it

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

Cheers, Paul.

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.

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.

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

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.

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

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.

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

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

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

I love google as well; the internet rocks for this kind of info, I can’t even remember what life was like beforehand!

Glad you enjoyed it!

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.

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.

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!

Hi Tomas,

Are you referring to this code block?

If so, I was just saying that you can simplify the two

ifconditions in the middle and the result is on the right.So

if (-P[i].y < -1)can be rewritten asif (P[i].y > 1)Cheers, Paul.

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!

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.

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.

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.

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?

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.

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.

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.

Hi Alan,

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

Cheers, Paul.

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.

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.

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.

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!

Hi, sorry I’m not sure what you mean?

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.

Hi Rich,

N is the contact normal involved in a collision, generated by the collision detection system.

Cheers, Paul.

So is the contact normal considered a distance or a force? Thanks

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.

Hi Joe,

No problem, glad they helped you!

1) Contacts are stored in a big, shared list

2) The distance between objects and normal are all you need if you’re only resolving about the COM of each object. If you need to consider angular effects as well, you’ll need the actual contact points too.

Cheers, Paul.

It works, this is awesome!

Thanks Paul

Great!

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.

Hi Time,

It’s somewhat more complex; I’d recommend following Chris Hecker’s articles here:

http://chrishecker.com/Rigid_Body_Dynamics

These cover the topic thoroughly.

Hope that helps!

Cheers, Paul.

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

basicconcepts 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

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!

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.

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.

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!

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

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.

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.

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.

Will you accept bitcoins for the source code?

Hi Elias,

I’m afraid not… I was tempted but really I could use the cash right now

Cheers, Paul.

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.

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.

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

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.

Hi Paul.

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

Hi Paul.

How about friction?

How can we apply friction to this simulation?

Thanks.

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.

Hi Paul.

So we dont need any cof ?

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

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!

Hi Felipe,

Thanks for pointing that out

Cheers, Paul.

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

Surrounding like a box

I’m also confused about what a normal is, you said that the planes have normals that point inward but what are they?

Unit length vector: http://www.wildbunny.co.uk/blog/vector-maths-a-primer-for-games-programmers/vector/#Unit

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?

Hi Mathias,

There’s a few related articles – have a look through these and see if any fit the bill?

http://www.wildbunny.co.uk/blog/category/physics-2/

Cheers, Paul.

Thank you for answering so fast.

The one about collision detection really helped.

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

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.

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)top*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!

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.

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.

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.

“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!

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?

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

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.

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.

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 .

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.

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 .

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

Hi Nathan,

.x is the x position.

Cheers, Paul.

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.

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.

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 ?

Have you seen my collision detection article? That has a demo of stacked balls without jittering problem:

http://www.wildbunny.co.uk/blog/2011/04/20/collision-detection-for-dummies/

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.

Hey Paul,

I know this is 4 years old now, but I’ve recently been following along and creating my own 2d physics engine to relative success. By that I mean it detects collisions, and resolves them, but there is still some sinking when multiple objects are resting on one another (eg. a bunch of circles gathering at the bottom of a box, the circles on top slowly push the bottom circles into the frame of the box)

This image shows enough circles to line the bottom of the box, and that they are not penetrating the box’s frame (or each other for that matter):

http://s6.postimg.org/s9ww2pco1/Circles_Lining_Box.jpg

This image shows more circles added in order to pile them up, leading to the top circles pushing the bottom circles through the box’s frame (as well as other circles). Note: this happens very slowly:

http://s6.postimg.org/eisf0hnq9/Circles_Penetrating_Box.jpg

I have had a look at your article on Continuous Collision Detection using a speculative approach. However, I’m looking for a way of preventing this penetration without needing to assuming a 0 restitution, as I require ‘bounciness’ in my engine.

My current collision resolution code is as follows:

`// Detect and resolve collisions`

for each (Contact contact in contacts)

{

RigidBody* a = contact.a;

RigidBody* b = contact.b;

Vector2 n = contact.normal;

float distance = contact.dist;

if (distance GetVelocity(), n) GetVelocity(), -n) GetRestitution();

if (b->GetRestitution() GetRestitution();

// Calculate the impulse

Vector2 v = a->GetVelocity() - b->GetVelocity();

Vector2 impulse = (1.f + e) * n * Dot(v, n) / (a->GetInvMass() + b->GetInvMass());

`// Distribute impulse based on mass`

a->ApplyForce(-impulse * a->GetInvMass());

b->ApplyForce(impulse * b->GetInvMass());

}

}

As you can probably see, I’m using your method described in this article for my collision detection and resolution. I am just stumped as to how I can implement penetration prevention/correction.

Any help would be greatly appreciated.

Thanks for the great article by the way.

– David

Sorry, the code block seems to have removed any part of the code that contained a less than symbol followed later by a greater than symbol, and since I’m using C++ and pointers, those symbols come up a fair bit. Here’s the full code (not in a code block):

// Detect and resolve collisions

for each (Contact contact in contacts)

{

RigidBody* a = contact.a;

RigidBody* b = contact.b;

Vector2 n = contact.normal;

float distance = contact.dist;

if (distance GetVelocity(), n) GetVelocity(), -n) GetRestitution();

if (b->GetRestitution() GetRestitution();

// Calculate the impulse

Vector2 v = a->GetVelocity() – b->GetVelocity();

Vector2 impulse = (1.f + e) * n * Dot(v, n) / (a->GetInvMass() + b->GetInvMass());

// Distribute impulse based on mass

a->ApplyForce(-impulse * a->GetInvMass());

b->ApplyForce(impulse * b->GetInvMass());

}

}

Sorry for the multiple posts, it’s done it again, so I’ll just pass a link to a screenshot of the code:

http://s6.postimg.org/iulm38k0x/Collision_Detection_And_Resolution.jpg

There are many different ways to resolve penetration, each with their own trade offs – the simplest is to add a small fraction of the penetration distance as a bias to the relative normal velocity (n * Dot(v,n) in your code. So you add a little bit of the penetration and over the course of a few frames the objects get pushed part.

Thanks for the quick response!

I change the impulse line from:

`Vector2 impulse = (1.f + e) * n * Dot(v, n) / (a->GetInvMass() + b->GetInvMass());`

to

`Vector2 impulse = (1.f + e) * n * (Dot(v, n) + penetration * 4.f) / (a->GetInvMass() + b->GetInvMass());`

And, the addition of penetration fixed my issue, although slower than I’d hoped (there was still visible penetration that lingered while being corrected). However, multiplying the penetration further, seems to have helped with that.

Is this what you meant? Or have I done something that may break something else later on?

Thanks,

– David

That’s about right – the factor you use to scale the penetration will depend on your in game units. The more you iterate over all the contacts, the better the convergence should be.

Thanks, you’ve been a great help!

It’s working quite well and is a bit of a balancing act between number of iterations and penetration scale factor. I’m thinking it’s less performance intensive if you go for a higher penetration scale factor and a lower iteration count.

The issue here is the jitter scaling with it. I’m wondering if there’s a logical way of detecting when the objects are no longer correcting penetration, but continually overcorrecting, so that I can cancel out their velocities.

Another thought is that there may be a way that I can introduce a leniency to the penetration correction, such that a minor penetration is not corrected.

Thanks,

– David

Look into accumulated impulses – this will reduce your need to correct for penetration. You can also try to google split-impulses, which attempts to separate the penetration correction from the velocity resolution.

Hello, Paul!

Thank you for the greate tutorial!

I have a question about constraints. As I understand, the constraint implementation provided in the article doesn’t conserve the full energy system. Even if set COR = 1 and no friction (and assuming no gravity and other forces), ball’s speeds eventually decreasing when balls hitting walls (when flying between, full energy is fine). Am I missing something, or, if not, how to fix it?

Hi Andrey,

It could be related to the integration scheme – are you seeing this drop off very slowly, or more quickly?

Cheers, Paul.

Hello!

I have created several youtube videos to demonstrate this effect:

http://www.youtube.com/watch?v=oWkQrb4q1sw

In this demo two circles are connected with

“equality constaint”, based on your article: the distance

between them is constant.

COR = 1, no friction, no gravity. The string in the low-left

corner is a sum of kinetic energies of linear and rotating

movement of these circles. As you can see, the energy is decreasing each

time some of the circles hit any wall.

When they fly between, energy is not changing.

http://www.youtube.com/watch?v=nXD4nMlDtjc

This is the same demo, but the constraint is disabled, circles are moving independently.

The energy in this case remains constant with very high precision.

The integration scheme is based on another tutorial. In short:

0. find collisions

1. calculate acceleration and increase velocities by acceleration * dt * 0.5

(why 0.5: http://www.niksula.hut.fi/~hkankaan/Homepages/gravity.html. But this doesn’t really matter, as we have no forces, so

acceleration is zero).

2. resolve collisions

3. resolve constraints

4. calculate new coordinates

5. calculate acceleration again and increase velocities by acceleration * dt * 0.5

Collisions are resolved in terms of additional velocities, similiar to constraints.

So for example, lets say we have one circle go strictly down and another circle is just above the first go down too and the equality constraint is set between them.

When the first circle hit the wall, its velocity is changed to the opposite. But another circle have the same velocity as before and

constraint resolution push the first circle again to the wall with slightly less speed, then again and again. And eventually (and very quickly) the circles just stop:

https://youtu.be/B8bgwG2__WQ

If the circles are moving in a random directions and hit the walls randomly, the speed drop is not so big but it takes place too.

The general problem as I see it is how to resolve the constraint if bodies velocities are opposite to each other? Because the method given in the article seems just nullify them.

Yeah, that looks totally wrong, I agree. It might be the positional correction code inside the constraint itself – it biases the resulting velocity of each body by a position correction factor (relDist/dt). This violates the energy conservation principle in exchange for computational speed. Try removing that factor from the constraint to see if it fixes the problem – if so you will probably want to use something like ‘split impulses’ to correct for positional errors without violating the conservation principle.

This is the first time I saw an article about a constraints solver that’s easy to follow and makes sense to a beginner like me. Thanks Paul!

/* Seems you didn’t want the visualization. Then here’s the experimental source code.

Note that force impact happens from two points at once because of aliasing problems in things that happen simultaneously. When force changes direction, therefore a little amount of relative velocity will get lost. First, i thought the loss can be explained with the gravity impact of the entire universe that consists of billions of galaxies on moving objects.

And it would be big on the winning street if it turned out that light is the quick and huge alternating current that appears because of the random oscillation of charges in heated objects.

Anyway, it generates a mesh that is glued together. You could try out and match rotation angles of a corresponding graphical pixel depiction, rotate all the pixels accordingly and then put them on the screen with a simple central projection.

Some elements could also be manipulated with scripts instead with the physics simulation which would then be carried out correspondingly. For example, the motor and the gear shifting would be controlled with an interpreted script, but the resistance from the road would cause a force on a particle that is rotated on an axis and represents the wheel back in the physics simulation.

*/

#include

#include

#define TZAHL 9

#define ELEK 0

#define PROT 1

#define NEUT 2

#define SCHEITEL 3.0

#define LADUNG 17.0

#define MAXKRAFT 4.0 // war beides 57.0

#define LADNULL 10000

#define TRAEGSCHWELLE 0.00

#define TIMESLICE 150.0

#define VERLANGSFAC 1.0 // war 0.9

int main(void)

{

long int runs=0;

struct

{

float x,y,z;

float xv,yv,zv;

float xv_old, yv_old ;

int elec_or_prot;

unsigned char symbol;

} teilchen[TZAHL];

int n=0, n2=0;

unsigned char screen[80][25];

int x,y;

float x_buf, y_buf, power, sign, distance;

teilchen[0].x=15.0, teilchen[0].y=5,

teilchen[0].xv=0.0, teilchen[0].yv=0;

teilchen[0].xv_old=0.0, teilchen[0].yv_old=0.0;

teilchen[0].elec_or_prot=PROT,

teilchen[0].symbol=’+';

teilchen[1].x=18.0, teilchen[1].y=5,

teilchen[1].xv=0.0, teilchen[1].yv=0.0;

teilchen[1].xv_old=0.0, teilchen[1].yv_old=0.0;

teilchen[1].elec_or_prot=ELEK,

teilchen[1].symbol=’-‘;

teilchen[2].x=21.0, teilchen[2].y=5,

teilchen[2].xv=0.0, teilchen[2].yv=0.0;

teilchen[2].xv_old=0.0, teilchen[2].yv_old=0.0;

teilchen[2].elec_or_prot=PROT,

teilchen[2].symbol=’+';

teilchen[3+0].x=15.0, teilchen[3+0].y=8,

teilchen[3+0].xv=0.0, teilchen[3+0].yv=0;

teilchen[3+0].xv_old=0.0, teilchen[3+0].yv_old=0.0;

teilchen[3+0].elec_or_prot=ELEK,

teilchen[3+0].symbol=’-‘;

teilchen[3+1].x=18.0, teilchen[3+1].y=8,

teilchen[3+1].xv=0.0, teilchen[3+1].yv=0.0;

teilchen[3+1].xv_old=0.0, teilchen[3+1].yv_old=0.0;

teilchen[3+1].elec_or_prot=PROT,

teilchen[3+1].symbol=’+';

teilchen[3+2].x=21.0, teilchen[3+2].y=8,

teilchen[3+2].xv=0.0, teilchen[3+2].yv=0.0;

teilchen[3+2].xv_old=0.0, teilchen[3+2].yv_old=0.0;

teilchen[3+2].elec_or_prot=ELEK,

teilchen[3+2].symbol=’-‘;

teilchen[6+0].x=15.0, teilchen[6+0].y=11,

teilchen[6+0].xv=0.0, teilchen[6+0].yv=0;

teilchen[6+0].xv_old=0.0, teilchen[6+0].yv_old=0.0;

teilchen[6+0].elec_or_prot=PROT,

teilchen[6+0].symbol=’+';

teilchen[6+1].x=18.0, teilchen[6+1].y=11,

teilchen[6+1].xv=0.0, teilchen[6+1].yv=0.0;

teilchen[6+1].xv_old=0.0, teilchen[6+1].yv_old=0.0;

teilchen[6+1].elec_or_prot=ELEK,

teilchen[6+1].symbol=’-‘;

teilchen[6+2].x=21.0, teilchen[6+2].y=11,

teilchen[6+2].xv=0.0, teilchen[6+2].yv=0.0;

teilchen[6+2].xv_old=0.0, teilchen[6+2].yv_old=0.0;

teilchen[6+2].elec_or_prot=PROT,

teilchen[6+2].symbol=’+';

/*

teilchen[2].x=55.0, teilchen[2].y=10,

teilchen[2].xv=0.0, teilchen[2].yv=0.0;

teilchen[2].xv_old=0.0, teilchen[2].yv_old=0.0;

teilchen[2].elec_or_prot=ELEK,

teilchen[2].symbol=’-‘;

*/

while(1)

{

n=0;

while ( n < TZAHL)

{

n2=0;

while ( n2 < TZAHL)

{

if ( n2==n) { n2++; continue;}

x_buf=0, y_buf=0;

power=

sqrt (

( teilchen[n].x – teilchen[n2].x ) *

( teilchen[n].x – teilchen[n2].x )+

( teilchen[n].y – teilchen[n2].y) *

( teilchen[n].y – teilchen[n2].y) );

distance=power;

x_buf=(teilchen[n].x-teilchen[n2].x);

y_buf=(teilchen[n].y-teilchen[n2].y);

if ( power-SCHEITEL < 0.0 ) sign= -1.0; else sign= 1.0;

power=power-SCHEITEL;

if ( fabs(power) = LADNULL ) power=0;

if ( distance != 0 )

{

y_buf*=power/distance;

x_buf*=power/distance;

} else x_buf=0, y_buf=0;

if ( (distance=sqrt ( x_buf*x_buf+y_buf*y_buf )) > TRAEGSCHWELLE )

{

power= (distance-TRAEGSCHWELLE)/ distance;

x_buf*=power; y_buf*=power;

} else x_buf=0, y_buf=0 ;

teilchen[n].xv+=x_buf, teilchen[n].yv+=y_buf;

n2++;

}

n++;

}

n=0;

while (n < TZAHL)

{

teilchen[n].x+=teilchen[n].xv/TIMESLICE,

teilchen[n].y+=teilchen[n].yv/TIMESLICE;

teilchen[n].xv*=VERLANGSFAC, teilchen[n].yv*=VERLANGSFAC;

teilchen[n].xv=(teilchen[n].xv_old+teilchen[n].xv)/2.0;

teilchen[n].yv=(teilchen[n].yv_old+teilchen[n].yv)/2.0;

teilchen[n].xv_old=teilchen[n].xv;

teilchen[n].yv_old=teilchen[n].yv;

n++;

}

if ( (runs++)%10==0)

{

x=0;

while ( x < 80)

{

y=0;

while ( y < 25)

{

screen[x][y]=' ';

y++;

}

x++;

}

n=0;

while ( n 0 && teilchen[n].x 0 && teilchen[n].y < 19)

screen[(int)teilchen[n].x][(int)teilchen[n].y]=teilchen[n].symbol;

n++;

}

system("cls\n");

y=0;

while ( y < 23)

{

x=0;

while ( x < 78)

{

printf("%c",screen[x][y] );

x++;

}

printf("\n");

y++;

}

#define MABST(N) sqrt((teilchen[N].x-teilchen[4].x)*(teilchen[N].x-teilchen[4].x)+(teilchen[N].y-teilchen[4].y)*(teilchen[N].y-teilchen[4].y))

printf("%f %f %f\n"

"%f %f %f\n"

"%f %f %f\n",MABST(0), MABST(1), MABST(2),

MABST(3), MABST(4), MABST(5),

MABST(6), MABST(7), MABST(8));

if ( kbhit()) {

if (getch()=='z' ){ teilchen[3].x-=3, teilchen[5].x+=3;

}getch(); }

}

}

}

Puzzling this together could be a waste of time if the trash source code is already crippled from the formatting CGI tool. Sorry, but i’ve seen it running and with a little fixes and adaptions it could be very useful, at least for collision control.

power=

sqrt (

( teilchen[n].x – teilchen[n2].x ) *

( teilchen[n].x – teilchen[n2].x )+

( teilchen[n].y – teilchen[n2].y) *

( teilchen[n].y – teilchen[n2].y) );

distance=power;

x_buf=(teilchen[n].x-teilchen[n2].x);

y_buf=(teilchen[n].y-teilchen[n2].y);

if ( power-SCHEITEL < 0.0 ) sign= -1.0; else sign= 1.0;

power=power-SCHEITEL;

if ( fabs(power) < 0.001 ) power=sign*0.001;

After having calculated the distance on a foreign particle, it calculates the force of course which then changes the movement vector

power= ( LADUNG / (( power * power)+MAXKRAFT) ) *sign ;

if ( (teilchen[n].elec_or_prot==ELEK && teilchen[n2].elec_or_prot==PROT)||

(teilchen[n].elec_or_prot==PROT && teilchen[n2].elec_or_prot==ELEK)) x_buf*=-1.0, y_buf*=-1.0;