I recently found this wonderful article: "The Mystery of the Bouncing Ball" [Jonathan A. Feucht, 2011]. Feucht's article starts with a simple model of a ball that falls quadratically under gravity following a parabolic trajectory until it hits the ground, at which point its coefficient of restitution $r$ determines how much velocity is lost upon impact before starting the next parabolic trajectory. For $r<1$, the maximium height reached with each parabola decreases, so does the time between impacts. After a finite amount of time the maximum height of the bounces converges to zero. Feucht works out the formulas for $t_\text{final}$ as a function of the ball's initial velocity $v_0$ and the acceleration due to gravity1 $g$: $$ t_\text{final} = \frac{2 v_0}{-g \left(1-\sqrt{r}\right)} $$
Feucht also works out the continuous function which interpolates the bounce apexes as a function of time: $$ h(t) = \frac{-g}{2} \left(\frac{1-\sqrt{r}}{1+\sqrt{r}}\right)^2 \left(t - \frac{2 v_0}{g\left(1-\sqrt{r}\right)}\right)^2. $$
We can use this equation to compute the height of the ball at any time $t$ after it is dropped. For a given value of $t$, we first need to determine which bounce the ball is currently experiencing. This is given by $$ N(t) = \left\lfloor \frac{2 \log\left(1 - \frac{- t g \left(1-\sqrt{r}\right) }{2 v_0}\right)}{\log(r)} \right\rfloor. $$
The ball is mid parabolic flight in the $N(t)$-th bounce, but we need to now its flight time relative to the start of the bounce. The bounce start time as a function of the bounce index $n$ is given by: $$ S(n) = \frac{-2 v_0 \left(1 - r^{n/2}\right)}{g\left(1-\sqrt{r}\right)}. $$ So, we can say that the duration within the $N(t)$-th bounce is: $$ x = t - S(N(t)), $$ where I drop the dependence of $x$ on $t$ for brevity.
We need three points to fully determine a parabola. We know that at time $S(N(t))$ the ball is at height $0$ and at time $S(N(t)+1)$ the ball is back to height $0$. If the midway moment in this bounce is then $x_\text{mid} = \frac{1}{2}(S(N(t)) + S(N(t)+1))$, then Feucht's $h$ function gives us the height at the midway moment: $h_\text{mid} = h(x_\text{mid})$. For convenience of notation let's also introduce the duration of the bounce as: $s = S(N(t)+1) - S(N(t))$.
Thus, the height of the ball at time $t$ is given by the quadratic function: $$ y = \frac{4 h_\text{mid} \left(x \left(s - x\right)\right)}{s^2}. $$ Unless of course, $t≥t_\text{final}$, in which case the ball is on the ground and $y=0$.
Having a directly samplable function is really nice for procedural animation.
Relying on numerica integration is error prone and sensitive to parameters, especially the time-step value / sampling rate.
The $y$ function above is infinitely sampleable. It's also differentiable
(though I leave that as homework for the reader/symbolic toolkit/automatic
differentiation library). It makes creating the animation above very easy. We can also make a super slow-mo version:
We can also make a trippy animation of a bouncing ball where we "zoom-in" in time as the ball reaches $t_\text{final}$.
So far we've been assuming the ball starts at height $0$ at $t=0$ with an initial velocity $v_0$ and coefficient of restitution $r<1$ and acceleration due to gravity $g$. Given some initial height $\tilde{h}$ and velocity $\tilde{v}$ we can also solve for when this height is reached $\tilde{t}$ and what initial velocity $v_0$ would result in it (this really doesn't require any of the bouncing business above, but I found it useful for working with the bounce function as an easing curve for animation): $$ \tilde{t} = \frac{\tilde{v} - \sqrt{\tilde{v}^2 - 2 g \tilde{h}}}{g} \quad \text{and} \quad v_0 = \sqrt{\tilde{v}^2 - 2 g \tilde{h}}, $$ where lack of real values means that $g$ is too strong to reach the given pair.
Suppose you have a certain desired start velocity $v_0$ and you want to be sure that the bouncing is done at a desired $t_\text{final}$. Then you can solve for the coefficient of restitution $r$ as a function of $v_0$, $t_\text{final}$, and $g$: $$ r = \left(1 + \frac{2 v_0}{g t_\text{final}}\right)^2. $$ Or similarly, given $v_0$, $t_\text{final}$ and $r$ you can solve for $g$: $$ g = \frac{-2 v_0}{t_\text{final} \left(1 - \sqrt{r}\right)}. $$
Compared to alternative methods for direct evaluation of bounce-like behavior, such as: $$ y_\text{alt} = \left(e^{-t}\right) \cdot \left|4 \cos(5t)\right|, $$ as seen on GameDev StackExchange, function I derive above is albeit more complicted but has the nice property that it goes to zero at a known finite time and that the ball is always on some perfect parabolic trajectory not slowing down mid-air. Energy is "lost" sharply at bounces. And the frequency of impacts increases (as expected) rather than stays constant (unnatural).
Here are a few implementations. I'm sure you could get ChatGPT to convert this to your language of choice.
function [y,t_final] = bounce(t,v0,g,r)
q = 1-sqrt(r);
t_final = 2.*v0/(-g*(q));
h = @(t,g,r,v0) (-g./2) .* (q./(2-q)).^2 .* (t - ((2.*v0)./(-g.*q))).^2;
S = @(n) (2*v0/-g) .* (1 - r.^(n/2)) ./ q;
N = floor(2.* log(1 - (t.*q.*-g)./(2.*v0)) ./ log(r));
S0 = S(N);
S1 = S(N+1);
x = t-S0;
s = S1-S0;
t_mid = 0.5*(S0+S1);
H_max = h(t_mid,g,r,v0);
y = 4.*H_max./(s.^2) .* (x.*(s-x));
y(t>=t_final) = 0;
end
function bounce(t,v0,g,r)
{
const q = 1 - Math.sqrt(r);
const t_final = 2 * v0 / (-g * q);
function h (t,g,r,v0) {
return (-g / 2) * Math.pow(q / (2 - q), 2) * Math.pow(t - (2 * v0) / (-g * q), 2);
}
function S(n) {
return (2 * v0 / -g) * (1 - Math.pow(r, n / 2)) / q;
}
const N = Math.floor(2 * Math.log(1 - (t * q * -g) / (2 * v0)) / Math.log(r));
const S0 = S(N);
const S1 = S(N + 1);
const x = t - S0;
const s = S1 - S0;
const t_mid = 0.5 * (S0 + S1);
const H_max = h(t_mid, g, r, v0);
const y = t>=t_final ? 0 : 4 * H_max / Math.pow(s, 2) * (x * (s - x));
return {y: y, t_final: t_final};
};
import numpy as np
def bounce(t, v0, g, r):
t = np.asarray(t)
q = 1 - np.sqrt(r)
t_final = 2 * v0 / (-g * q)
def h(t):
return (-g / 2) * (q / (2 - q))**2 * (t - (2 * v0) / (-g * q))**2
def S(n):
return (2 * v0 / -g) * (1 - r**(n / 2)) / q
N = np.floor(2 * np.log(1 - (t * q * -g) / (2 * v0)) / np.log(r)).astype(int)
S0 = S(N)
S1 = S(N + 1)
x = t - S0
s = S1 - S0
t_mid = 0.5 * (S0 + S1)
H_max = h(t_mid)
y = 4 * H_max / s**2 * (x * (s - x))
# Set to zero where t ≥ t_final
y = np.where(t >= t_final, 0, y)
return y, t_final
CSS can't run an abitrary function but you can approximate an idealized bouncing ball very well. Using the functions above we can predefine keyframes for an arbitrary number of bounces before $t_\text{final}$. I found that 97% looks pretty good. We could just use `ease-out` and `ease-in` to approximate the parabolic up and down trajectories, but with `cubic-bezier` we can get a much closer approximation. Through numerical fitting I found that these work great:
cubic-bezier(0.020367,0.020367, 0.44008, 1); /* parabolic up */
cubic-bezier(0.55992, 0, 0.979633, 0.979633); /* parabolic down */
1 Feucht's article uses positive $g$, but I assume $v_0$ is positive and $g$ is negative (as in -9.8m/s²). Many of my equations include sign changes.