```  *     â (x^2 * y + z^2) / (x^2 + z^2),   x,   x * z * (y - 1) / (x^2 + z^2) â
* M = â  x * (y^2 - 1)  / (x^2 + z^2),   y,    z * (y^2 - 1)  / (x^2 + z^2) â
*     â x * z * (y - 1) / (x^2 + z^2),   z,   (x^2 + z^2 * y) / (x^2 + z^2) â
+ * ```
+ * + * This is stable as long as v (the bone) is not too much aligned with +/-Y + * (i.e. x and z components are not too close to 0). * - * This is stable as long as v (the bone) is not too much aligned with +/-Y (i.e. x and z components - * are not too close to 0). + * Since v is normalized, we have `x^2 + y^2 + z^2 = 1`, + * hence `x^2 + z^2 = 1 - y^2 = (1 - y)(1 + y)`. * - * Since v is normalized, we have x^2 + y^2 + z^2 = 1, hence x^2 + z^2 = 1 - y^2 = (1 - y)(1 + y). * This allows to simplifies M like this: + *
```  *     â 1 - x^2 / (1 + y),   x,     -x * z / (1 + y) â
* M = â                -x,   y,                   -z â
*     â  -x * z / (1 + y),   z,    1 - z^2 / (1 + y) â
+ * ```
* - * Written this way, we see the case v = +Y is no more a singularity. The only one remaining is the bone being - * aligned with -Y. + * Written this way, we see the case v = +Y is no more a singularity. + * The only one + * remaining is the bone being aligned with -Y. * - * Let's handle the asymptotic behavior when bone vector is reaching the limit of y = -1. Each of the four corner - * elements can vary from -1 to 1, depending on the axis a chosen for doing the rotation. And the "rotation" here - * is in fact established by mirroring XZ plane by that given axis, then inversing the Y-axis. - * For sufficiently small x and z, and with y approaching -1, all elements but the four corner ones of M - * will degenerate. So let's now focus on these corner elements. + * Let's handle + * the asymptotic behavior when bone vector is reaching the limit of y = -1. + * Each of the four corner elements can vary from -1 to 1, + * depending on the axis a chosen for doing the rotation. + * And the "rotation" here is in fact established by mirroring XZ plane by that given axis, + * then inversing the Y-axis. + * For sufficiently small x and z, and with y approaching -1, + * all elements but the four corner ones of M will degenerate. + * So let's now focus on these corner elements. * - * We rewrite M so that it only contains its four corner elements, and combine the 1 / (1 + y) factor: + * We rewrite M so that it only contains its four corner elements, + * and combine the `1 / (1 + y)` factor: + *
```  *                    â 1 + y - x^2,        -x * z â
* M* = 1 / (1 + y) * â                            â
*                    â      -x * z,   1 + y - z^2 â
+ * ```
+ * + * When y is close to -1, computing 1 / (1 + y) will cause severe numerical instability, + * so we ignore it and normalize M instead. + * We know `y^2 = 1 - (x^2 + z^2)`, and `y < 0`, hence `y = -sqrt(1 - (x^2 + z^2))`. * - * When y is close to -1, computing 1 / (1 + y) will cause severe numerical instability, so we ignore it and - * normalize M instead. We know y^2 = 1 - (x^2 + z^2), and y < 0, hence y = -sqrt(1 - (x^2 + z^2)). * Since x and z are both close to 0, we apply the binomial expansion to the first order: - * y = -sqrt(1 - (x^2 + z^2)) = -1 + (x^2 + z^2) / 2. Which gives: + * `y = -sqrt(1 - (x^2 + z^2)) = -1 + (x^2 + z^2) / 2`. Which gives: + *
```  *                        â  z^2 - x^2,  -2 * x * z â
* M* = 1 / (x^2 + z^2) * â                         â
*                        â -2 * x * z,   x^2 - z^2 â
+ * ```
```+       * (v)' = a(cceleration) =
```+       * v(t + dt) = v(t) + a(t) * dt