rigid body dynamics - How can I implement revolute (hinge) joints in a 2d physics system?

Further to my previous question on fixed joints, I am now going toe-to-toe with hinge joints in 2D physics. Though the answer to the previous question was excellent, I have come across a few problems in extending it to hinge joints.

My approach was to extend that given here by implementing the following algorithm:

1. Find the 'target' point for both hinge points by taking the midpoint between them (I actually use mass to give a ratio here but assume the bodies are of equal mass -- in my simulation they all are). The target point is labelled T in the diagram below.


2. For body a, find the angle between the hinge point A and the target point T. This is done by transforming T from global space to local space, then taking the difference of Atan2(y, x) for each point - this should return the angle between them.


3. The body a is then rotated by a matrix created from this angle, so that A and T now lie along the same vector. a is then translated by T-A; now A should be at T.


4. Repeat steps 2 and 3 for body b.

The algorithm seems logical to me, but I have found that there is some kind of stability issue whereby the bodies will stay constrained for a short time, then as the angle between the two bodies becomes slightly more acute (in either direction), they will both shoot off at high velocity. The number of constraint iterations is not the issue as far as I can tell - the same problems occur at 1 iteration as at 10.

Pseudocode to my algorithm given below:

// EntityA and EntityB are the two objects
// HingeA and HingeB are represented as vectors relative to their local transforms
Transform tA = EntityA.Transform;
Transform tB = EntityB.Transform;

Vector2 delta = tA.LocalToGlobal(HingeA) - tB.LocalToGlobal(HingeB);

// Locate the target in global space by finding the midpoint of the delta vector A-B
// then translating it by B

Vector2 Target = (delta * 0.5) + tB.LocalToGlobal(HingeB);

// Locate target in Entity A's local space
Vector2 TargetA = tA.GlobalToLocal(Target);

// Avoid undefined behaviour for Atan2
if (TargetA != Vector2.Zero && HingeA != Vector2.Zero)
{
   // Get the angular difference between A and T
   float angleA = Math.Atan2(TargetA.Y, TargetA.X) - Math.Atan2(HingeA.Y, HingeA.X);

   tA.Angle += angleA;
}

// Re-convert the target to local co-ordinates and translate EntityA towards it
tA.Pos += tA.GlobalToLocal(Target) - HingeA;

// Repeat for entity B...

Am I on the right track with this algorithm? Is there some minor detail I have missed that is aggravating the simulation?