+/*
+ * based on: https://box2d.org/files/ErinCatto_NumericalMethods_GDC2015.pdf,
+ * page 76.
+ */
+static void rb_solve_gyroscopic( rigidbody *rb, m3x3f I, f32 h ){
+ /* convert to body coordinates */
+ v3f w_local;
+ m3x3_mulv( rb->to_local, rb->w, w_local );
+
+ /* Residual vector */
+ v3f f, v0;
+ m3x3_mulv( I, w_local, v0 );
+ v3_cross( w_local, v0, f );
+ v3_muls( f, h, f );
+
+ /* Jacobian */
+ m3x3f iJ, J, skew_w_local, skew_v0, m0;
+ m3x3_skew_symetric( skew_w_local, w_local );
+ m3x3_mul( skew_w_local, I, m0 );
+
+ m3x3_skew_symetric( skew_v0, v0 );
+ m3x3_sub( m0, skew_v0, J );
+ m3x3_scalef( J, h );
+ m3x3_add( I, J, J );
+
+ /* Single Newton-Raphson update */
+ v3f w1;
+ m3x3_inv( J, iJ );
+ m3x3_mulv( iJ, f, w1 );
+ v3_sub( w_local, w1, rb->w );
+
+ m3x3_mulv( rb->to_world, rb->w, rb->w );
+}
+