
Has anyone modeled angular constraints in Verlet with the same stability as stick constraints? It seems like angular constraints are prone to a lot of instability, if they are resolved in the same way as sticks  that is, if the angle is not correct between AB and BC then set it to the correct angle by moving A and C each half the total distance between the current angle and the required angle.
When angular constraints are stressed they seem to generate a lot of extra energy. This is helped a bit by softening them, but to some degree the problem persists.
Another problem is leverage is not modeled correctly  e.g., angle ABC with AB length of 10 and BC length of 100, constrainted to say Math.PI / 2, and with particle B 'pinned' in place. BC will not be 'pulled down' lower due to gravity. In fact, the opposite seems to occur.
Below is my code. Forgive the homebrewed the trig there  Im sure it could be tightened up. I used an additional reference line to get angles, which is probably not necessary.
Also, what is the correct way to add inverse mass values to the resolution?
AngularConstraint.prototype.resolve = function() {
// make sure the reference line position gets updated
// lineC is pBpD
this.lineC.p2.x = this.lineC.p1.x + 0;
this.lineC.p2.y = this.lineC.p1.y  1;
abRadius = this.pA.distance(this.pB);
bcRadius = this.pB.distance(this.pC);
thetaABC = this.calcTheta(this.pA, this.pB, this.pC);
thetaABD = this.calcTheta(this.pA, this.pB, this.pD);
thetaCBD = this.calcTheta(this.pC, this.pB, this.pD);
halfTheta = (this.targetTheta  thetaABC) / 2;
paTheta = thetaABD + halfTheta * this.coeffStiff;
pcTheta = thetaCBD  halfTheta * this.coeffStiff;
this .pA .x = abRadius * Math .sin(paTheta ) + this .pB .x ;
this .pA .y = abRadius * Math .cos(paTheta ) + this .pB .y ;
this .pC .x = bcRadius * Math .sin(pcTheta ) + this .pB .x ;
this .pC .y = bcRadius * Math .cos(pcTheta ) + this .pB .y ;
}
AngularConstraint.prototype.calcTheta = function(pa, pb, pc) {
AB = new Vector(pb.x  pa.x, pb.y  pa.y);
BC = new Vector(pc.x  pb.x, pc.y  pb.y);
dotProd = AB.dot(BC);
crossProd = AB.cross(BC);
return Math .atan2(crossProd , dotProd );
}

