EDIT:
This series ended up on Tuts+, and these are sort of deprecated. Please visit Tuts+ to see the finalized articles.
Personal Update
Beyond C++ reflection I’ve taken a huge liking to physics programming for games. Reflection in C++ is a largely “solved” area of study. Though resources on the internet may be few and far between for creating custom introspection, people that know how to create such systems have explored much of what is to be accomplished.
However physics is another story entirely. Up until just a few years ago physics for games was in a highly primitive state; techniques just hadn’t been created due to the inability of hardware to perform necessary computation. There are a lot of new and exciting problems to solve.
So this next year while attending DigiPen I’ll be focusing my studies around: engine architecture, data oriented design and optimization, C++ introspection and game physics. This last topic being the topic I’m going to start an article series on.
Game Physics High Level Detail
A game physics engine performs some functionality for ensuring that shapes behave in a particular way. To state this one could say: a physics engine removes degrees of freedom from moving bodies, and enforces these imposed rules.
From here on out, until specified otherwise I’ll be talking specifically about convex rigid bodies within a physics engine. A rigid body is just a shape that is implicitly rigid; the engine does not naturally support deformation of a set of points that make up a shape. Using this implicit definition the mathematics behind manipulating such shapes can be known intuitively and implemented efficiently.
Here’s a short feature list of some features a custom physics system may provide, assuming the engine uses rigid bodies:
 Collision detection
 Collision resolution
 Linear movement, angular rotation and integration
 Raycasting
 Islanding/sleeping/portals
 Friction simulation
 Spatial partitioning/bounding volumes
 Fracturing/splitting
 Multipart bodies
 Various types of shapes
 Advanced mathematical constraints (joints, motors, contact constraints, etc.)
Raycasting from Box2D.
In this article series I will attempt to talk about all of the above topics in my own time. I’ve learned a lot about physics and physics engines in a short amount of time and by documenting what I have learned I hope to solidify my understanding, as well as help out others to achieve what I have. I would not have learned nearly as much as I currently know without help from others so it is only natural to want to do the same.
The best example of an open source physics engine that employs much of the feature list shown above would be Box2D by Erin Catto. The rest of this article series will be detailing the physics engine that I myself have written. There are of course other methods I choose not to talk about, or just don’t know about.
Architecture
There are two main objects that make up a physics engine: shapes and bodies. A body is a small package of data that defines intrinsic properties of a physics object. Properties such as density, restitution (elasticity), friction coefficients, inertia tensors, along with any other flags or settings should be stored within the body. These bits of data are properties that can be isolated away from the shape definition. The shape itself is contained with the body through a pointer.
The shape definition defines what sort of shape a physics object represents. Some physics engines can attach multiple shapes to a body, which are referred to as fixtures. A shape stores any data required to define the shape itself, and provides various functions to operate on the shape, such as testing for line or point intersection, or generating a bounding volume.
Together the body and shape represent a physics object, which by the end of this article series will be able to perform a lot of interesting interactions with other physics objects.
A tower of stacked oriented boxes within a custom physics engine called iiD.
All bodies should be contained within what is known as a scene or world. I refer to this object as a scene. The scene should contain a list of all live objects, as well as functionality inserting or removing bodies from the scene. The scene should also have callbacks for processing shape or ray queries. A query just checks to see if any bodies intersect with something like a point or ray.
The scene has one particular function called step, which steps the scene forward in time by a single delta time (dt). This step function steps all objects forward in time by integration. The integration step just moves the objects forward by using their velocity, position and acceleration to determine their next position.
During the step collisions are detected and then resolved. Often times a broadphase of some sort is used to detect possible collisions, and expensive collisions operations are only used when really needed.
All of this organization allows the user of your physics engine to worry about three main operations: creating and removing bodies, and stepping the scene. The rest is automated and taken care of within the physic system’s internals.
The last major isolated system would be the broadphase. There are two major phases in collision detection: the broad and narrow phases. The narrow phase is the one which determines if two shapes intersect or not. The broad phase determines if two shapes can possibly be intersecting or not. The broadphase should be constructed such that intersection queries are very very fast and simple. An axisaligned bounding box (AABB) will suffice perfectly.
Once the broadphase chooses which objects could perhaps be colliding, it sends them off to the narrow phase. The narrow phase performs much more intensive operations to determine if two objects are colliding or not. The whole point of this is to reduce the amount of times the narrow phase has to be used, as it is expensive to operate with.
Once the narrow phase finds a pair of bodies to be colliding information about the collision is gathered and placed into a small container called the manifold. Do not ask why it is called a manifold, for I have no idea and neither does anyone else seem to! The manifold contains three important pieces of information:
 Collision penetration depth
 Direction to resolve in
 Points of contact
These pieces of information are used to fill in formulas that are used to resolve the penetration and solve for a single unknown: the magnitude of the resolution vector. Here’s a small diagram:
c – Collision contact
n – Resolution normal
d – Penetration distance
It is also useful to to store pointers or handles to the two objects that formed this collision info. This allows some useful functions to placed into the manifold object directly: solve and resolve. Solve would be the function to collect the three pieces of collision information. If no contact points are found, then there is no collision. If contact points are found, then resolving performs a resolution operation to force both objects to not be intersecting after integration.
Velocity
Complex physics manipulation is performed on the velocities of objects. Trying to manipulate the positions of objects directly is very difficult, as it poses a problem that isn’t linear. By using derived position equations for velocity, the problem is thus simplified. Most of the time we will be only dealing with velocity manipulation.
Impulse Resolution in 2D (No Rotation or Friction)
The act of resolving collisions is something that isn’t covered very well on the internet. There are some scattered documents and details, though the information isn’t always easy to find. Most of what I know I learned by talking with other people, but I know most people will not have such an opportunity. Hopefully this article series can provide a strong resource for understanding and constructing a simple physics engine.
The toughest place to code is in my opinion resolution of collision. There exists tons of information on collision detection and broadphase, and thus creating these portions of a physics engine is in my opinion not too difficult. Some resources for collision detection are: RealTime Collision Detection by Christer Ericson, and Game Physics Engine Developement by Ian Millington. I have both of these books right next to me as I write this :)
Generating contact manifolds and resolving such manifolds are what most programmers get caught up in. So lets hit the ground running and tackle a portion of code that will bring your entire physics system to life: contact resolution.
The best type of contact resolution to start with is impulse resolution. The idea behind impulse resolution is to take your contact manifold and solve for a single velocity vector that can be used to add into an object’s current velocity, in order to make it not penetrating the other object during the next frame. The reason for starting with impulse resolution is that it’s quite simple and easy to get up and running, and more complicated and advanced techniques require you to understand impulse resolution anyway.
2D impulse simulation.
Now the contact manifold is assumed to contain the direction of our velocity vector we are solving for. I will cover how to generate the contact manifold in another article in this series. The only unknown left to solve for is the magnitude of this vector. It so happens that it’s possible to solve for this magnitude in relation to both objects. If this magnitude is known, you add velocity along the normal scaled by the magnitude for one object, and you do the same operation to the other object in the direction opposite to the manifold normal. Lets start from the top.
We have two objects moving along in our scene, and there exists a relative velocity between the two, with object A as the reference object at the origin:

Eq 1: VelocityRelative = VelocityA  VelocityB 
The relative velocity can be expressed in terms of the collision normal (from the collision manifold) with a dot product:

Eq 2: VelocityRelative dot ManifoldNormal = (VelocityA  VelocityB) dot ManifoldNormal 
This can be thought of as the relative normal velocity between the two objects. The next thing to include in this derivation is the coefficient of restitution (elasticity factor). Newton’s Law of Restitution states that:

Eq 3: VelocityBefore * restitution = VelocityAfter 
Often times the restitution identifier is specified by an e, or epsilon symbol. Knowing this it’s fairly simple to include it within our current equation:

Eq 4: VelocityRelativeAfter dot ManifoldNormal = resitution(VelocityA  VelocityB) dot ManifoldNormal 
Now we need to go to another topic and model an impulse. I have said the term “impulse” quite a few times without defining it, so here is the definition I use: an impulse is an instantaneous change in velocity. We will use an impulse to change the velocity of an object. Here’s how you could use an impulse to modify the velocity of a single object:

Eq 5: VelocityNew = VelocityOld * Impulse 
Here Impulse would be a scalar floating point value. This isn’t too useful however, as it only scales an object’s current velocity, and thusly makes the object move slower or faster along the positive or negative direction of the vector.
What is needed is a way to do componentwise modification of the vector, so we can make it point slightly in one direction or another, allowing objects to make turns and change directions.

Eq 6: VelocityNew = VelocityOld + Impulse(Direction) 
In the above equation we can take a direction vector with a magnitude of 1, and scale it by our impulse. By adding this new vector to our velocity we can then modify the velocity in any way we wish. Our end goal is to solve for our Impulse scalar that will separate two objects from collision, so in the end we’ll need to distribute this scalar across two equations in terms of velocity.
Lets start with a simple momentum equation:

Eq 7: Momentum = Mass * Velocity 
An impulse is defined to be a change in momentum. Thus we get:

Eq 8: Impulse = MomentumAfter  MomentumBefore Impulse = Mass * VelocityAfter  Mass * VelocityBefore 
To isolate our velocity after we can rearrange into:

Eq 9: VelocityNew = VelocityOld + Impulse(Direction) / Mass 
Now lets change equation 4 into one that contains velocities under the influence of impulses. However we’ll want to express our VelocityNew as one that is acted upon by impulse, and substitute in equation 9:

Eq 10: VelocityA + Impulse(Direction) / MassA  VelocityB + Impulse(Direction) / MassB = Restitution(VelocityRelativeAtoB) * Direction Simplification: VelocityRelativeAtoB * Direction + Impulse( Direction / MassA + Direction / MassB ) + Restitution(VelocityRelativeAtoB) * Direction = 0 Isolate Impulse: Impulse = (1 + Restitution)(VelocityRelativeAtoB dot Direction)  1 + 1   MassA MassB 
Remember that the impulse is a scalar. Also note that all values on the right hand side of the equation are all known, including the Direction which was solved for by the collision detection.
All that is left here is to distribute this scalar impulse over each object’s velocity vector proportional to the total mass of the system (system being collision between both objects).
The total mass is MassA + MassB, so to get an even distribution you do: impulse * Mass / TotalMass. To simplify this one could use the following:

Eq 11 (applying an impulse proportional to mass of object): VelocityNew = VelocityOld + (1 / Mass) * Impulse(Direction) 
This can be done twice, once per object. The total impulse will be applied, except only a portion of the impulse will be applied to each object. This ensures smaller objects move more than larger ones during impact.
One thing you must ensure is that if the velocities are separating (objects moving away from one another) that you do nothing. Here’s a sample version of finalized code for impulse resolution in 2D without rotation or friction:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

// j = (1 + epsilon) * N dot Vrel //  // invM_1 + invM_2 void Manifold::ApplyImpulse( void ) { Vec2 rv = B>m_velocity  A>m_velocity; real contactVel = Vec2::DotProduct( rv, normal ); // Do not resolve if velocities are separating if(contactVel > 0) return; // Calculate restitution real e = std::min(A>m_material.restitution, B>m_material.restitution); // Calculate impulse scalar real j = (1.0f + e) * contactVel; j /= A>m_massdata.inv_mass + B>m_massdata.inv_mass; // Apply impulse Vec2 impulse = j * normal; A>m_velocity = A>m_massdata.inv_mass * impulse; B>m_velocity += B>m_massdata.inv_mass * impulse; } 
Rotation
Now that we have covered resolution without these two factors adding them in to our final equation (equation 10) will be quite a bit simpler. We will need to understand the concept of “mass” in terms of rotations. The inertia tensor of an object describes how “hard it is to turn” along some axis. In 2D there’s only a single axis of rotation (along the z) so a single tensor is all that’s needed. Here’s a new version of equation 11.

Eq 12: RotationalVelocityNew = RotationalVelocityOld + (1 / InertiaTensor) * PointRelativeToCenterofMass cross Direction * Impulse Shortened: V' += I^1 * r x (n * j) 
For the sake of brevity here’s the final equation, where r is the vector from center of mass to a point on the body (contact point). The velocity of a point on the body relative to it’s center is given by:

Eq 13: Vpoint = VCenterMass + AngularVelocity cross r 
Final equation (Direction substituted for n):

Eq 14: Impulse = (1 + Restitution) * (VelocityRelativeAtoB dot n)  1 1 (rA cross n)^2 (rB cross n)^2  +  +  +  MassA MassB InertiaTensorA InertiaTensorB 
And there we have it! This will solve for a separating impulse given a specific contact point. Now you might have noticed there are a couple cross products that are kinda strange in 2D. I won’t cover what they are here, but they still hold true. I’ll just give them to you:

Cross vector and scalar in 2D: V' = x * V.y, x * V.x Cross scalar and vector in 2D: V' = x * V.y, x * V.x 

Cross two vectors in 2D: A.x * B.y  A.y * B.x 
Friction
Friction is the simplest thing to do in this entire resolution article, assuming we’re dealing with 2D! I actually haven’t derived this portion myself, so again I’m going to just throw the equations at you. Assuming we have our resolution direction in the manifold (the collision normal), we can calculate a tangent vector, which is a vector perpendicular to the collision normal. Just replace all instances of n in equation 14 with this tangent vector.
To solve for the tangent vector you do:

TangentVelocity = VRelAtoB  (VRelAtoB dot n) * n TangentVector = Normalize( TangentVelocity ) 
Again, just take the above TangentVector and replace it in equation 14 for n.
There is something I have missed. When solving for the force of friction there’s a max, which is the manifold normal * the coefficient of friction (from the body definition). This is due to the Coloumb friction law. Treat the coefficient of friction * normal as your “friction cap” when solving for your tangential impulse. This is just a simple capping of the final impulse vector.
Penetration Resolution
Due to floating point error energy will always be lost from the system over time. As such, when objects start coming to a rest they will sink into each other due to gravity. How to solve? Well the solution is to actually just handle the leftover penetration manually by just moving the objects apart directly. This will prevent objects from sinking into one another, though doesn’t add any energy into the system.
To resolve penetration, simply move each object along the manifold normal in opposite directions. However there’s an art to doing so; you need to be gentle. If you just resolve 100% of the penetration each frame objects underneath other objects will jitter violently due to the naive penetration resolution scheme I am presenting. To avoid this, try only resolving a percentage of the leftover penetration, perhaps 20 to 80%.
Additionally you only want to resolve penetration if the penetration depth after the impulse is applied is above some constant arbitrary (but small) threshold, (try 0.01 to 0.1). If the penetration is below this, then don’t move either object.
This method of penetration resolution is called linear projection. Here’s a snippet of C++ code demonstrating this:

void Manifold::PositionalCorrection( void ) { const real k_slop = 0.01f; // Penetration allowance const real percent = 0.2f; // Penetration percentage to correct Vec2 correction = (std::max( penetration  k_slop, 0.0f ) / (A>m_massdata.inv_mass + B>m_massdata.inv_mass)) * percent * normal; A>m_tx.position = A>m_massdata.inv_mass * correction; B>m_tx.position += B>m_massdata.inv_mass * correction; } 
Iterative Solving
There is one one additional tweak you can do to increase the believability of your physics simulation: iterate over all contacts and solve + resolve the impulses many times (perhaps 5 to 20 iterations). Since the large equation 14 has relative velocity within it, each iteration will feed in the previous result to come up with a new one.
This will allow the step to propagate energy throughout multiple objects contacting one another within a single timestep. This is essential for allowing the “billiards balls” effect to ensue.
This simple iteration is a very easy way to vastly improve the results of resolution.
Resources/Links
 http://chrishecker.com/images/e/e7/Gdmphys3.pdf
 http://www.cs.cmu.edu/~baraff/sigcourse/notesd2.pdf