Implement new gyro motion integration method

[unpairedbracket]

Replace the "jolt" method for rotational integration with a new one that has better energy preservation (exact for bodies with two or three identical principal moments of inertia)

Objective

The existing "gyroscopic" modification to angular velocity integration has good properties (angular momentum magnitude conservation, seemingly fairly good energy conservation), but it causes the energy of freely rotating bodies to approach $\frac{L^2}{2 I_2}$ - corresponding to the "intermediate" principal moment of inertia. This means that a freely rotating asymmetric body will approach the "tennis racket effect" behaviour over time, even if initially rotating around one of the physically stable axes.

The objective here is to replace the existing solve_gyroscopic_torque with an implementation that exhibits better energy conservation than this

Solution

My proposed solution is described in (excessive) detail in a blog post here. It exhibits better energy-conservation properties than any other method I have tested against, including the "Second order" and "Enhanced second-order" methods published by S. Buss here.

Its energy conservation is probably not as good as that of the discrete Moser-Veselov method (described e.g. here that claims exact energy conservation, but I haven't managed to figure out how I would actually implement that method to compare.

Testing

This PR adds a new example, symmetric_top, that shows off the failure mode of the current gyroscopic integrator. A symmetric body (two distinct principal moments, with one of them repeated twice) like a bar should move by spinning around its axis while its axis rotates around the angular momentum vector. With the existing method the object's axis gradually moves away from the angular momentum vector until it is perpendicular to it; the new method preserves the angle between the two.

The new example logs the rotational energy and the magnitude of the angular momentum at each time step, showing that the new method preserves energy very closely, whereas the old method causes energy loss for the long, thin bar object and energy gain for the thin, wide plate object.

There were (rightly) performance questions about this new method, so I have done some performance testing (not full benchmarking but I can go there if desired). I tiled the bodies in the symmetric_top example in a 100x100 grid, meaning 20,000 bodies total, and turned down the substep count to 1 to minimise per-substep overhead. I watched "Integrate Velocities" on the Diagnostics UI, and saw roughly:

Current method: 0.60ms (30ns per body)
Proposed method 0.75ms (37.5ns per body, +25%)
Proposed method, with fast rotation approximation: 0.66ms (33ns per body, +10%)

For testing, -F new_gyro will activate the proposed method, and additionally activating -F approx_rotation will swap from Quat::from_scaled_axis to the approximate rotation method I describe in the same blog post, which looks to be a few nanoseconds faster per call than creating the quaternion, using it and then throwing it away.

If this is approved I will of course remove the unused functionality and the feature flags, they're just there to make comparisons easier for review purposes