+/*
+ * Solve a list of positional constraints
+ */
+VG_STATIC void rb_solve_position_constraints( rb_constr_pos *buf, int len ){
+ for( int i=0; i<len; i++ ){
+ rb_constr_pos *constr = &buf[i];
+ rigidbody *rba = constr->rba, *rbb = constr->rbb;
+
+ v3f wa, wb;
+ m3x3_mulv( rba->to_world, constr->lca, wa );
+ m3x3_mulv( rbb->to_world, constr->lcb, wb );
+
+ m3x3f ssra, ssrat, ssrb, ssrbt;
+
+ m3x3_skew_symetric( ssrat, wa );
+ m3x3_skew_symetric( ssrbt, wb );
+ m3x3_transpose( ssrat, ssra );
+ m3x3_transpose( ssrbt, ssrb );
+
+ v3f b, b_wa, b_wb, b_a, b_b;
+ m3x3_mulv( ssra, rba->w, b_wa );
+ m3x3_mulv( ssrb, rbb->w, b_wb );
+ v3_add( rba->v, b_wa, b );
+ v3_sub( b, rbb->v, b );
+ v3_sub( b, b_wb, b );
+ v3_muls( b, -1.0f, b );
+
+ m3x3f invMa, invMb;
+ m3x3_diagonal( invMa, rba->inv_mass );
+ m3x3_diagonal( invMb, rbb->inv_mass );
+
+ m3x3f ia, ib;
+ m3x3_mul( ssra, rba->iIw, ia );
+ m3x3_mul( ia, ssrat, ia );
+ m3x3_mul( ssrb, rbb->iIw, ib );
+ m3x3_mul( ib, ssrbt, ib );
+
+ m3x3f cma, cmb;
+ m3x3_add( invMa, ia, cma );
+ m3x3_add( invMb, ib, cmb );
+
+ m3x3f A;
+ m3x3_add( cma, cmb, A );
+
+ /* Solve Ax = b ( A^-1*b = x ) */
+ v3f impulse;
+ m3x3f invA;
+ m3x3_inv( A, invA );
+ m3x3_mulv( invA, b, impulse );
+
+ v3f delta_va, delta_wa, delta_vb, delta_wb;
+ m3x3f iwa, iwb;
+ m3x3_mul( rba->iIw, ssrat, iwa );
+ m3x3_mul( rbb->iIw, ssrbt, iwb );
+
+ m3x3_mulv( invMa, impulse, delta_va );
+ m3x3_mulv( invMb, impulse, delta_vb );
+ m3x3_mulv( iwa, impulse, delta_wa );
+ m3x3_mulv( iwb, impulse, delta_wb );
+
+ v3_add( rba->v, delta_va, rba->v );
+ v3_add( rba->w, delta_wa, rba->w );
+ v3_sub( rbb->v, delta_vb, rbb->v );
+ v3_sub( rbb->w, delta_wb, rbb->w );
+ }