As explained previously, the current algorithm does not handle collisions of multiple balls at the same time. I found a solution for a double collision, and I am pretty sure that there is no single solution for more complex cases using only conservation laws and the coefficient of restitution formula – but I still have to prove it. Of course there would be a solution modeling the actual physics of the collisions – a problem for another day.
Double ball collision
..or more formally: elastic collision of three spheres in two dimensions
Ball 0 – with velocity \((v_0,w_0)\), hits at the same time ball 1 and ball 2 – that have velocities \((v_1, w_1)\) and \((v_2, w_2)\) respectively.
Projecting the velocity \((v_0,w_0)\) onto the lines connecting the center of ball 0 with the centers of balls 1 and 2 gives the velocity components \(q_{01}\) and \(q_{02}\). Projecting the velocity \((v_1,w_1)\) onto the line connecting the center of ball 0 with the center of ball 1 and onto its perpendicular through the center of ball 1 gives the velocity components \(q_1\) and \(p_1\). Doing the same thing for ball 2 gives the velocity components \(q_2\) and \(p_2\).
Here is the drawing:
As usual, balls are not spinning and there is no ball-ball friction, so the forces between balls are only along the lines connecting the center of ball 0 with the centers of balls 1 and 2. This means that the velocity components perpendicular to those lines are the same before and after the collision:
$$ p’_1 = p_1 $$
$$ p’_2 = o_2 $$
(indicating with a prime the velocities after the collision), whereas the velocity components along those lines are related by the coefficient of restitution formula:
$$ q’_{01} – q’_1 = c(q_1 – q_{01}) $$
$$ q’_{02} – q’_2 = c(q_2 – q_{02}) $$
Let’s not use trigonometry this time. The projection of a velocity onto the line connecting the center of two ball is the scalar product of the velocity with the unit vector pointing from one center to the other. For balls 0 and 1 this unit vector is:
$$ \left( \frac{x_{01}}{\sqrt{x_{01}^2 + y_{01}^2}}, \frac{y_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} \right)$$
where \(x_{01} = x_1 – x_0\) and \(y_{01}=y_1-y_0\), hence:
$$ q_1 = \frac{v_1x_{01} + w_1y_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} $$
$$ q_01 = \frac{v_0x_{01} + w_0y_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} $$
Similarly, the projection on the perpendicular of that line is the scalar product of the velocity with the perpendicular unit vector:
$$ \left( -\frac{y_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} , \frac{x_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} \right) $$
hence:
$$ p_1 = \frac{-v_1y_{01} + w_1x_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} $$
Replacing these values in the four initial equations gives:
$$ –v’_1y_{01} + w’_1x_{01} = -v_1y_{01} + w_1x_{01} $$
$$ –v’_2y_{02} + w’_2x_{02} = -v_2y_{02} + w_2x_{02} $$
$$ v’_0x_{01} + w’_0y_{01} – v’_1x_{01} – w’_1y_{01} = c(v_1x_{01} + w_1y_{01} - v_0x_{01} - w_0y_{01}) $$
$$ v’_0x_{02} + w’_0y_{02} – v’_1x_{02} – w’_2y_{02} = c(v_2x_{02} + w_2y_{02} - v_0x_{02} - w_0y_{02}) $$
The conservation of momentum gives two additional equations:
$$ m_0v’_0 + m_1v’_1 + m_2v’_2 = m_0v_0 + m_1v_1 + m_2v_2 $$
$$ m_0w’_0 + m_1w’_1 + m_2w’_2 = m_0w_0 + m_1w_1 + m_2w_2 $$
- where \(m_0\), \(m_1\) and \(m_2\) are the masses of the balls, resulting in a system of 6 linear equations for the 6 unknowns \(v’_0, w’_0, v’_1, w’_1, v’_2, w’_2\). Introducing new unknowns \(u_j\) defined as:
$$ u_0 = v’_0 – v_0,\quad u_1 = w’_0 – w_0,\quad u_2 = v’_1 – v_1,\quad u_3 = w’_1 – w_1,\quad u_4 = v’_2 – v_2,\quad u_5 = w’_2 – w_2 $$
the system becomes a bit simpler:
$$ x_{01} u_3 – y_{01} u_2 = 0 $$
$$ x_{02} u_5 – y_{02} u_4 = 0 $$
$$ x_{01} (u_0-u_2) + y_{01}(u_1-u_3) = B_{01} $$
$$ x_{02} (u_0-u_4) + y_{02}(u_1-u_5) = B_{02} $$
$$ m_0u_0 + m_1u_2 + m_2u_4 = 0 $$
$$ m_0u_1 + m_1u_3 + m_2u_5 = 0 $$
with
$$ B_{01} = (1+c)[x_{01}(v_1 – v_0) + y_{01}(w_1-w_0)] $$
$$ B_{02} = (1+c)[x_{02}(v_2 – v_0) + y_{02}(w_2-w_0)] $$
After quite a lot of substitutions and simplification the result is
$$ u_0 = v’_0 – v_0 = –\frac{1}{\Delta}[B_{01}m_1(x_{02}m_2P – x_{01}(m_0+m_2)S_{02}) + B_{02}m_2(x_{01}m_1P-x_{02}(m_0+m_1)S_{01})] $$
$$ u_1= w’_0 – w_0 = –\frac{1}{\Delta}[B_{01}m_1(y_{02}m_2P – y_{01}(m_0+m_2)S_{02}) + B_{02}m_2(y_{01}m_1P-y_{02}(m_0+m_1)S_{01})] $$
$$ u_2 = v’_1 – v_1 = \frac{x_{01}m_0}{\Delta}[B_{02}m_2P-B_{01}(m_0+m_2)S_{02})] $$
$$ u_3= w’_1 – w_1 = \frac{y_{01}m_0}{\Delta}[B_{02}m_2P-B_{01}(m_0+m_2)S_{02})] $$
$$ u_4 = v’_2 – v_2 = \frac{x_{02}m_0}{\Delta}[B_{01}m_1P-B_{02}(m_0+m_1)S_{01})] $$
$$ u_5 = w’_2 – w_2 = \frac{y_{02}m_0}{\Delta}[B_{01}m_1P-B_{02}(m_0+m_1)S_{01})] $$
with
$$ P = x_{01}x_{02}+y_{01}y_{02} $$
$$ S_{01} = x_{01}^2 + y_{01}^2 $$
$$ S_{02} = x_{02}^2 + y_{02}^2 $$
$$ \Delta = (m_0 + m_1)(m_0 + m_2) S_{01}S_{02}-m_1m_2P^2 $$
The corresponding code is:
/**
* Updates the velocities of this ball and two other balls after a double collisions
* with both balls at the exact same time
* The coordinate of the balls must be at the collision point.
* @param otherBall1 second colliding ball
* @param otherBall2 third colliding ball
* @param restitution coefficient of restitution for a ball-ball collision
*/
collide2(otherBall1: Ball, otherBall2: Ball, restitution: number) {
var m0 = Math.pow(this.radius, 3);
var m1 = Math.pow(otherBall1.radius, 3);
var m2 = Math.pow(otherBall2.radius, 3);
var x01 = otherBall1.x - this.x;
var y01 = otherBall1.y - this.y;
var x02 = otherBall2.x - this.x;
var y02 = otherBall2.y - this.y;
var p = x01 * x02 + y01 * y02;
var s01 = x01 * x01 + y01 * y01;
var s02 = x02 * x02 + y02 * y02;
var delta = (m0 + m1) * (m0 + m2) * s01 * s02 - m1 * m2 * p * p;
var v01 = otherBall1.v - this.v;
var w01 = otherBall1.w - this.w;
var b01 = (1 + restitution) * (x01 * v01 + y01 * w01);
var v02 = otherBall2.v - this.v;
var w02 = otherBall2.w - this.w;
var b02 = (1 + restitution) * (x02 * v02 + y02 * w02);
this.v = -(b01 * m1 * (x02 * m2 * p - x01 * (m0 + m2) * s02) + b02 * m2 * (x01 * m1 * p - x02 * (m0 + m1) * s01)) / delta + this.v;
this.w = -(b01 * m1 * (y02 * m2 * p - y01 * (m0 + m2) * s02) + b02 * m2 * (y01 * m1 * p - y02 * (m0 + m1) * s01)) / delta + this.w;
var r1 = m0 * (b02 * m2 * p - b01 * (m0 + m2) * s02) / delta;
otherBall1.v = x01 * r1 + otherBall1.v;
otherBall1.w = y01 * r1 + otherBall1.w;
var r2 = m0 * (b01 * m1 * p - b02 * (m0 + m1) * s01) / delta;
otherBall2.v = x02 * r2 + otherBall2.v;
otherBall2.w = y02 * r2 + otherBall2.w;
} // collide2
As a verification, when \(m_2 = 0\) the first four equations become:
$$ u_0 = v’_0 – v_0 = x_{01} \frac{m_1}{m_0+m_1} \frac{B_{01}}{S_{01}} $$
$$ u_1 = w’_0 – w_0 = y_{01} \frac{m_1}{m_0+m_1} \frac{B_{01}}{S_{01}} $$
$$ u_2 = v’_1 – v_1 = -x_{01} \frac{m_0}{m_0+m_1} \frac{B_{01}}{S_{01}} $$
$$ u_3 = w’_1 – w_1 = -y_{01} \frac{m_0}{m_0+m_1} \frac{B_{01}}{S_{01}} $$
if \(\alpha\) is the angle between the vector from ball 0 to ball 1 and the horizontal, then:
$$ \cos\alpha = \frac{x_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} $$
$$ \sin\alpha = \frac{y_{01}}{\sqrt{x_{01}^2 + y_{01}^2}} $$
hence
$$ \frac{B_{01}}{S_{01}} = -\frac{1+c}{\sqrt{x_{01}^2 + y_{01}^2}}[\cos\alpha(v_0 – v_1) + \sin\alpha(w_0-w_1)] $$
and the result is the same as in the two balls case.