Billiards Physics Simulation

You have 2 balls, A and B. Initially, you know the positions of both balls and is represented by the 2 vectors P_A and P_B. The initial velocity of each ball is represented by the vectors V_A and V_B. The friction acting on each ball is represented by an acceleration vector pointing in the opposite direction of the velocity vector, and of a constant magnitude. These friction vectors will be called F_A and F_B.

We can express the position of a ball as a function of time with the following formula:

F(p, v, f, t) = p + vt + 0.5ft²

p is the position of the ball at t=0
v is the velocity of the ball at t=0
f is the vector representing the friction acting on the ball
t is the time

We can use this formula to express the positions of balls A and B as a function of time:

f(P_A, V_A, F_A, t) = P_A + V_At + 0.5F_At² f(P_B, V_B, F_B, t) = P_B + V_Bt + 0.5F_Bt²

Now we can define a function that will express the distance between two balls as a function of time:

distance(P_A, V_A, F_A, P_B, V_B, F_B, t) = |f(P_A, V_A, F_A, t) - f(P_B, V_B, F_B, t)|

To determine when a collision occurs all we have to do is set this equation equal to a balls diameter (2 * radius) and solve for t.

The first step is rewriting these equations so they use real numbers as opposed to vectors.

The position function f(p, v, f, t) can be broken down into two functions that calculate the balls position along the x-axis as a function of time, and the balls position along the y-axis as a function of time:

f_x(p, v, f, t) = p_x + v_xt + 0.5f_xt² f_y(p, v, f, t) = p_y + v_yt + 0.5f_yt²

In order to rewrite the distance formula, recall that |Z| denotes the magnitude of a vector Z. The magnitude can be calculated by using the x and y components like so:

|Z| = sqr(Z_x² + Z_y²)

Using this we can rewrite our distance formula as follows:

distance(P_A, V_A, F_A, P_B, V_B, F_B, t) = sqr((f_x(P_A, V_A, F_A, t) - f_x(P_B, V_B, F_B, t))² + (f_y(P_A, V_A, F_A, t) - f_y(P_B, V_B, F_B, t))²)

Setting this equal to our balls diameter (D) we get:

D = sqr((f_x(P_A, V_A, F_A, t) - f_x(P_B, V_B, F_B, t))² + (f_y(P_A, V_A, F_A, t) - f_y(P_B, V_B, F_B, t))²)

Squaring each side this becomes

D² = (f_x(P_A, V_A, F_A, t) - f_x(P_B, V_B, F_B, t))² + (f_y(P_A, V_A, F_A, t) - f_y(P_B, V_B, F_B, t))²

Substituting in the formulas for f_x and f_y we get:

D² = ((P_A.x + V_A.xt + 0.5F_A.xt²) - (P_B.x + V_B.xt + 0.5F_B.xt²))² + ((P_A.y + V_A.yt + 0.5F_A.yt²) - (P_B.y + V_B.yt + 0.5F_B.yt²))²

Now lets remove some of those brackets

D² = (P_A.x + V_A.xt + 0.5F_A.xt² - P_B.x - V_B.xt - 0.5F_B.xt²)² + (P_A.y + V_A.yt + 0.5F_A.yt² - P_B.y - V_B.yt - 0.5F_B.yt²)²

Now we have to evaluate the squares which will produce a monster of an equation with 72 terms:

D² = P_A.x² + P_A.xV_A.xt + P_A.x0.5F_A.xt² - P_A.xP_B.x - P_A.xV_B.xt - P_A.x0.5F_B.xt² + V_A.xtP_A.x + V_A.x²t² + V_A.xt³0.5F_A.x - V_A.xtP_B.x - V_A.xt²V_B.x - V_A.xt³0.5F_B.x + 0.5F_A.xt²P_A.x + 0.5F_A.xt³V_A.x + 0.25F_A.x²t⁴ - 0.5F_A.xt²P_B.x - 0.5F_A.xt³V_B.x - 0.25F_A.xt⁴F_B.x - P_B.xP_A.x - P_B.xV_A.xt - P_B.x0.5F_A.xt² + P_B.x² + P_B.xV_B.xt + P_B.x0.5F_B.xt² - V_B.xtP_A.x - V_B.xt²V_A.x - V_B.xt³0.5F_A.x + V_B.xtP_B.x + V_B.x²t² + V_B.xt³0.5F_B.x - 0.5F_B.xt²P_A.x - 0.5F_B.xt³V_A.x - 0.25F_B.xt⁴F_A.x + 0.5F_B.xt²P_B.x + 0.5F_B.xt³V_B.x + 0.25F_B.x²t⁴ + P_A.y² + P_A.yV_A.yt + P_A.y0.5F_A.yt² - P_A.yP_B.y - P_A.yV_B.yt - P_A.y0.5F_B.yt² + V_A.ytP_A.y + V_A.y²t² + V_A.yt³0.5F_A.y - V_A.ytP_B.y - V_A.yt²V_B.y - V_A.yt³0.5F_B.y + 0.5F_A.yt²P_A.y + 0.5F_A.yt³V_A.y + 0.25F_A.y²t⁴ - 0.5F_A.ytP_B.y - 0.5F_A.yt²V_B.y - 0.25F_A.yt⁴F_B.y - P_B.yP_A.y - P_B.yV_A.yt - P_B.y0.5F_A.yt² + P_B.y² + P_B.yV_B.yt + P_B.y0.5F_B.yt² - F_B.ytP_A.y - V_B.yt²V_A.y - V_B.yt³0.5F_A.y + V_B.ytP_B.y + V_B.y²t² + V_B.yt³0.5F_B.y - 0.5F_B.yt²P_A.y - 0.5F_B.yt³V_A.y - 0.25F_B.yt⁴F_A.y + 0.5F_B.yt²P_B.y + 0.5F_B.yt³V_B.y + 0.25F_B.y²t⁴

First lets simplify by grouping like terms:

D² = P_A.x² + 2P_A.xV_A.xt + P_A.xF_A.xt² - 2P_A.xP_B.x - 2P_A.xV_B.xt - P_A.xF_B.xt² + V_A.x²t² + V_A.xF_A.xt³ - 2V_A.xP_B.xt - 2V_A.xV_B.xt² - V_A.xF_B.xt³ + 0.25F_A.x²t⁴ - F_A.xP_B.xt² - F_A.xV_B.xt³ - 0.5F_A.xF_B.xt⁴ + P_B.x² + 2P_B.xV_B.xt + P_B.xF_B.xt² + V_B.x²t² + V_B.xF_B.xt³ + 0.25F_B.x²t⁴ + P_A.y² + 2P_A.yV_A.yt + P_A.yF_A.yt² - 2P_A.yP_B.y - 2P_A.yV_B.yt - P_A.yF_B.yt² + V_A.y²t² + V_A.yF_A.yt³ - 2V_A.yP_B.yt - 2V_A.yV_B.yt² - V_A.yF_B.yt³ + 0.25F_A.y²t⁴ - F_A.yP_B.yt - F_A.yV_B.yt² - 0.5F_A.yF_B.yt⁴ + P_B.y² + 2P_B.yV_B.yt + P_B.yF_B.yt² + V_B.y²t² + V_B.yF_B.yt³ + 0.25F_B.y²t⁴

Now lets rewrite in Quartic (biquadratic) form:

(0.25F_A.x² - 0.5F_A.xF_B.x + 0.25F_B.x² + 0.25F_A.y² - 0.5F_A.yF_B.y + 0.25F_B.y²)t⁴ + (V_A.xF_A.x - V_A.xF_B.x - F_A.xV_B.x + V_B.xF_B.x + V_A.yF_A.y - V_A.yF_B.y - F_A.yF_B.y - F_A.yV_B.y + V_B.yF_B.y)t³ + (P_A.xF_A.x - P_A.xF_B.x + V_A.x² - 2V_A.xV_B.x - F_A.xP_B.x + P_B.xF_B.x + V_B.x² + P_A.yF_A.y - P_A.yF_A.y - P_A.yF_B.y + V_A.y² - 2V_A.yV_B.y - F_A.yP_B.y + P_B.yF_B.y + V_B.y²)t² + (2P_A.xV_A.x - 2P_A.xV_B.x - 2V_A.xP_B.x + 2P_B.xV_B.x + 2P_A.yV_A.y - 2P_A.yV_B.y - 2V_A.yP_B.y + 2P_B.yV_B.y)t + (P_A.x² - 2P_A.xP_B.x + P_B.x² + P_A.y² - 2P_A.yP_B.y + P_B.y² - D²) = 0

We can make the following substitutions to pretty up our equation:

a = 0.25F_A.x² - 0.5F_A.xF_B.x + 0.25F_B.x² + 0.25F_A.y² - 0.5F_A.yF_B.y + 0.25F_B.y² b = V_A.xF_A.x - V_A.xF_B.x - F_A.xV_B.x + V_B.xF_B.x + V_A.yF_A.y - V_A.yF_B.y - F_A.yF_B.y - F_A.yV_B.y + V_B.yF_B.y c = P_A.xF_A.x - P_A.xF_B.x + V_A.x² - 2V_A.xV_B.x - F_A.xP_B.x + P_B.xF_B.x + V_B.x² + P_A.yF_A.y - P_A.yF_A.y - P_A.yF_B.y + V_A.y² - 2V_A.yV_B.y - F_A.yP_B.y + P_B.yF_B.y + V_B.y² d = 2P_A.xV_A.x - 2P_A.xV_B.x - 2V_A.xP_B.x + 2P_B.xV_B.x + 2P_A.yV_A.y - 2P_A.yV_B.y - 2V_A.yP_B.y + 2P_B.yV_B.y e = P_A.x² - 2P_A.xP_B.x + P_B.x² + P_A.y² - 2P_A.yP_B.y + P_B.y² - D²

We can divide the co-efficient of t⁴ into all terms to eliminate it. We will do this with some more substitutions:

g = b/a h = c/a i = d/a j = e/a

After the substitutions our equation is in the following form:

t⁴ + gt³ + ht² + it + j = 0

To find the times that the collision occurs we must find the roots of this equation. There is an algorithm for finding the roots of a quartic equation which I will apply here.

First we must rewrite it as its resolvent cubic equation and solve for one of the real roots. We can do this like so:

y³ - hy² + (gi - 4j)y - g²j + 4hj - i² = 0

We can make some substitutions to pretty this up:

k = -h m = gi - 4j n = -g²j + 4hj - i²

This gives us the form:

y³ + ky² + my + n = 0

Now we apply the method to solve cubic equations which requires us to rewrite this equation yet again in the following form:

x³ + px + q = 0 Where x = y + k/3, p = (1/3)(3m - k²), q = (1/27)(2k³ - 9km + 27n)

We now calculate 2 values r and s like so:

The roots of this cubic equation are then:

x = r + s x = -((r + s) / 2) + ((r + s) / 2)sqr(-3) x = -((r + s) / 2) + ((r - s) / 2)sqr(-3)

At least one of these roots must be real. Find any real root and use it to calculate y. Recall:

y = x - (k / 3)

Now we calculate a value u:

u = sqr((g² / 4) - h + y)

Now we calculate 2 values w and z. The forumulas used depend on the value of u.
If u != 0 then

w = sqr(((3g²) / 4) - u² - 2h + ((4gh - 8i - g³) / 4u)) z = sqr(((3g²) / 4) - u² - 2h - ((4gh - 8i - g³) / 4u))

If u = 0 then

w = sqr(((3g²) / 4) - 2h + 2sqr(y² - 4j)) z = sqr(((3g²) / 4) - 2h - 2sqr(y² - 4j))

The roots of our quartic equation are given by:

t = -(g / 4) + (u / 2) +/- (w / 2) t = -(g / 4) - (u / 2) +/- (z / 2)

This will give us between 0 and 4 distinct roots. We have to eliminate some of them. We can calculate the time at which the ball will stop moving due to friction:

t_f = - v / f

We can now throw away any roots that don't satisfy 0 <= t < t_f

We should be left with at most 2 real roots.

If there are 0 or 1 real roots then no collision occurs.
If there are 2 distinct real roots then the collision occurs at the smaller of the 2 roots.