Billiard balls collision Given the 2-dimensional positions and velocities of a pair of billiard balls right before impact, calculate their velocities after a perfectly elastic collision. The balls are assumed to be ideal spheres (or equivalently: circles) with the same radius, same mass, uniform density, and no friction.

Input consists of 8 numbers: p0x,p0y,v0x,v0y,p1x,p1y,v1x,v1y where p0x,p0y is the centre of the first ball, v0x,v0y its velocity, and similarly p1x,p1y,v1x,v1y for the second ball. You can accept input in any order and structured in any convenient way, e.g. as a 2x2x2 array, or maybe a 2x2 array for p and two length-2 arrays for v0 and v1. It's also fine to take complex numbers (if your language supports them) instead of xy pairs. However, you should not take input in a coordinate system other than Cartesian, i.e. polar is not allowed.

Note that the radius of a billiard ball is half the distance between p0x,p0y and p1x,p1y, so it's not given as an explicit part of the input.

Write a program or function that outputs or returns 4 numbers in any convenient Cartesian representation: the post-collision values of v0x,v0y,v1x,v1y.



A possible algorithm is:


find the normal line that passes through both centres
find the tangent line that passes through the midpoint between the two centres and is perpendicular to the normal line
change coordinate system and break down v0x,v0y and v1x,v1y into their tangential and normal components v0t,v0n and v1t,v1n
swap the normal components of v0 and v1, preserving their tangential components
change back to the original coordinate system


Tests (results rounded to 5 decimal places):

   p0x   p0y   v0x   v0y   p1x   p1y   v1x   v1y -&gt;      v0x'       v0y'       v1x'       v1y'
[-34.5,-81.8, 34.7,-76.1, 96.2,-25.2, 59.2,-93.3] [  49.05873, -69.88191,  44.84127, -99.51809]
[ 36.9, 77.7,-13.6,-80.8, -7.4, 34.4, 15.1,-71.8] [   5.57641, -62.05647,  -4.07641, -90.54353]
[-51.0, 17.6, 46.1,-80.1, 68.6, 54.0,-35.1,-73.9] [ -26.48927,-102.19239,  37.48927, -51.80761]
[-21.1,-52.6,-77.7, 91.5, 46.0, 94.1, 83.8, 93.7] [ -48.92598, 154.40834,  55.02598,  30.79166]
[ 91.3, -5.3, 72.6, 89.0, 97.8, 50.5, 36.2, 85.7] [  71.73343,  81.56080,  37.06657,  93.13920]
[-79.9, 54.9, 92.5,-40.7,-20.8,-46.9,-16.4, -0.9] [  47.76727,  36.35232,  28.33273, -77.95232]
[ 29.1, 80.7, 76.9,-85.1,-29.3,-49.5,-29.0,-13.0] [  86.08581, -64.62067, -38.18581, -33.47933]
[ 97.7,-89.0, 72.5, 12.4, 77.8,-88.2, 31.5,-34.0] [  33.42847,  13.97071,  70.57153, -35.57071]
[-22.2, 22.6,-61.3, 87.1, 67.0, 57.6,-15.3,-23.1] [ -58.90816,  88.03850, -17.69184, -24.03850]
[-95.4, 15.0,  5.3, 39.5,-54.7,-28.5, -0.7,  0.8] [  21.80656,  21.85786, -17.20656,  18.44214]
[ 84.0,-26.8,-98.6,-85.6,-90.1, 30.9,-48.1, 37.2] [ -89.76828, -88.52700, -56.93172,  40.12700]
[ 57.8, 90.4, 53.2,-74.1, 76.4,-94.4,-68.1,-69.3] [  51.50525, -57.26181, -66.40525, -86.13819]
[ 92.9, 69.8,-31.3, 72.6,-49.1,-78.8,-62.3,-81.6] [-123.11680, -23.48435,  29.51680,  14.48435]
[-10.3,-84.5,-93.5,-95.6, 35.0, 22.6, 44.8, 75.5] [ -11.12485,  99.15449, -37.57515,-119.25449]
[ -3.9, 55.8,-83.3,  9.1, -2.7,-95.6, 37.7,-47.8] [ -82.84144, -48.75541,  37.24144,  10.05541]
[-76.5,-88.4,-76.7,-49.9, 84.5, 38.0,  4.2, 18.4] [   6.52461,  15.43907, -79.02461, -46.93907]
[ 64.2,-19.3, 67.2, 45.4,-27.1,-28.7, 64.7, -4.3] [  59.66292,  44.62400,  72.23708,  -3.52400]
[  9.8, 70.7,-66.2, 63.0,-58.7, 59.5, 83.7,-10.6] [  68.07646,  84.95469, -50.57646, -32.55469]
[ 62.9, 46.4, 85.0, 87.4, 36.3,-29.0,-63.0,-56.3] [  23.53487, -86.82822,  -1.53487, 117.92822]
[ -5.5, 35.6, 17.6,-54.3, -2.2, 66.8,-15.2, 11.8] [  24.15112,   7.63786, -21.75112, -50.13786]


Shortest wins. No loopholes.



thanks @Anush for helping fix the diagram's background colour
Python 3, 67 66 bytes, 53 bytes

def f(p,v,q,w):p-=q;d=((v-w)/p).real*p;return v-d,w+d

Try it online!
-1 byte thanks to @ngn
-13 bytes thanks to @Neil
This function f takes four complex numbers as input and returns two complex numbers. The ungolfed version is shown in the following.
Ungolfed

def elastic_collision_complex(p1, v1, p2, v2):
    p12 = p1 - p2
    d = ((v1 - v2) / p12).real * p12
    return v1 - d, v2 + d

Try it online!
The computation formula is derived based on the 2D-vector formula on wiki.
Since \$m_1=m_2\$, the formula can be simplified to
$$\left\{\begin{array}{l} v_1'=v_1-dv \\ v_2'=v_2+dv\end{array}\right.$$
Let \$x_{12}=x_1-x_2\$ and \$v_{12}=v_1-v_2\$, we have
$$dv=\frac{\langle v_{12}, x_{12}\rangle}{\|x_{12}\|^2}x_{12} = \frac{Re(v_{12}\cdot\overline{x_{12}})}{x_{12}\cdot \overline{x_{12}}}x_{12}=Re\left(\frac{v_{12}\cdot\overline{x_{12}}}{x_{12}\cdot \overline{x_{12}}}\right)x_{12}=Re\left(\frac{v_{12}}{x_{12}}\right)x_{12}$$
In the ungolfed program, p12, v1 - v2, d correspond to \$x_{12}\$, \$y_{12}\$, and \$dv\$, respectively.