+VG_STATIC void rb_constraint_limits( rigidbody *ra, v3f lca,
+ rigidbody *rb, v3f lcb, v3f limits[2] )
+{
+ v3f ax, ay, az, bx, by, bz;
+ m3x3_mulv( ra->to_world, (v3f){1.0f,0.0f,0.0f}, ax );
+ m3x3_mulv( ra->to_world, (v3f){0.0f,1.0f,0.0f}, ay );
+ m3x3_mulv( ra->to_world, (v3f){0.0f,0.0f,1.0f}, az );
+ m3x3_mulv( rb->to_world, (v3f){1.0f,0.0f,0.0f}, bx );
+ m3x3_mulv( rb->to_world, (v3f){0.0f,1.0f,0.0f}, by );
+ m3x3_mulv( rb->to_world, (v3f){0.0f,0.0f,1.0f}, bz );
+
+ v2f px, py, pz;
+ px[0] = v3_dot( ay, by );
+ px[1] = v3_dot( az, by );
+
+ py[0] = v3_dot( az, bz );
+ py[1] = v3_dot( ax, bz );
+
+ pz[0] = v3_dot( ax, bx );
+ pz[1] = v3_dot( ay, bx );
+
+ float r0 = atan2f( px[1], px[0] ),
+ r1 = atan2f( py[1], py[0] ),
+ r2 = atan2f( pz[1], pz[0] );
+
+ /* calculate angle deltas */
+ float dx = 0.0f, dy = 0.0f, dz = 0.0f;
+
+ if( r0 < limits[0][0] ) dx = limits[0][0] - r0;
+ if( r0 > limits[1][0] ) dx = limits[1][0] - r0;
+ if( r1 < limits[0][1] ) dy = limits[0][1] - r1;
+ if( r1 > limits[1][1] ) dy = limits[1][1] - r1;
+ if( r2 < limits[0][2] ) dz = limits[0][2] - r2;
+ if( r2 > limits[1][2] ) dz = limits[1][2] - r2;
+
+ v3f wca, wcb;
+ m3x3_mulv( ra->to_world, lca, wca );
+ m3x3_mulv( rb->to_world, lcb, wcb );
+
+ rb_limit_cure( ra, rb, ax, dx );
+ rb_limit_cure( ra, rb, ay, dy );
+ rb_limit_cure( ra, rb, az, dz );
+}
+
+VG_STATIC void rb_debug_constraint_position( rigidbody *ra, v3f lca,
+ rigidbody *rb, v3f lcb )
+{
+ v3f wca, wcb;
+ m3x3_mulv( ra->to_world, lca, wca );
+ m3x3_mulv( rb->to_world, lcb, wcb );
+
+ v3f p0, p1;
+ v3_add( wca, ra->co, p0 );
+ v3_add( wcb, rb->co, p1 );
+ vg_line_pt3( p0, 0.005f, 0xffffff00 );
+ vg_line_pt3( p1, 0.005f, 0xffffff00 );
+ vg_line( p0, p1, 0xffffff00 );
+}
+
+VG_STATIC void rb_constraint_position( rigidbody *ra, v3f lca,
+ rigidbody *rb, v3f lcb )
+{
+ /* C = (COa + Ra*LCa) - (COb + Rb*LCb) = 0 */
+ v3f wca, wcb;
+ m3x3_mulv( ra->to_world, lca, wca );
+ m3x3_mulv( rb->to_world, lcb, wcb );
+
+ v3f rcv;
+ v3_sub( ra->v, rb->v, rcv );
+
+ v3f rcv_Ra, rcv_Rb;
+ v3_cross( ra->w, wca, rcv_Ra );
+ v3_cross( rb->w, wcb, rcv_Rb );
+ v3_add( rcv_Ra, rcv, rcv );
+ v3_sub( rcv, rcv_Rb, rcv );
+
+ v3f delta;
+ v3f p0, p1;
+ v3_add( wca, ra->co, p0 );
+ v3_add( wcb, rb->co, p1 );
+ v3_sub( p1, p0, delta );
+
+ float dist2 = v3_length2( delta );
+
+ if( dist2 > 0.00001f )