Use a Taylor approximation for velocity damping

[Jondolf]

Currently, linear and angular velocity damping are handled with a Padé approximation of exponential decay. This formulation is also used by engines such as Rapier and Box2D. Box2D has a brief derivation for the math.

// Exponential
body.lin_vel.0 *= (-delta_secs * lin_damping.0).exp();

// Pade approximation
body.lin_vel.0 *= 1.0 / (1.0 + delta_secs * lin_damping.0);

This is run at every substep for entities with LinearDamping and AngularDamping, during velocity integration.

The Padé approximation is used to avoid calling exp, which is comparatively expensive and also has non-deterministic precision, meaning that a slower software implementation would likely be required.

An even faster option that avoids a division is the first-order Taylor approximation:

// First-order Taylor approximation
body.lin_vel.0 *= 1.0 - delta_secs * lin_damping.0;

This PR changes velocity damping to use this Taylor approximation.

Accuracy

In my testing, the Taylor approximation produces results that are still extremely close to the exponential function when running at reasonable time steps. Here I am simulating 1,000 objects with damping coefficients ranging from 1 (top) to 0 (bottom) with an initial linear velocity towards the right.

Blue is the exponential function, red is the Palé approximation, and green is the Taylor approximation. Even with just 1 substep and a time step of 10 Hz, the results it settles at are very close:

With a more reasonable but still conservative setup of 4 substeps and a time step of 50 Hz, they are practically indistinguishable:

Based on these results, it seems like the Taylor approximation would provide good enough accuracy while being extremely cheap.

Aside: Caching

We can observe that the damping formula only requires the damping coefficient and delta time for the damping factor. When using a fixed time step (which is recommended for physics), these values typically remain constant. If we cached this factor, and only updated it when the damping coefficient or time step change, we could reduce damping to a single multiplication:

body.lin_vel.0 *= lin_damping.cached_factor;

This is something that should be benchmarked and tested to see if the caching is worth storing the extra float.

[Jondolf]

So this has the problem that the Taylor approximation can become negative if the damping coefficient is large enough. It is rare to hit this, since with a typical 4 substeps and a time step of 60 Hz, you would need a damping coefficient of over 240 for the damping to become negative, while coefficients typically are just between 0 and 1. Still, large coefficients are technically valid, and there is the risk that people get weird behavior because of it, especially with lower time steps.

We could clamp the value to be non-negative, but at that point, just the Palé approximation might make more sense.

So for now I think I will close this. The performance difference is likely negligible. We can consider caching the term if it is measured to improve things at all.