mimoSite Programming and some other things

I just finished Coursera’s Exploring Quantum Physics. It was pretty good: the Coursera Web sites works well, both lectures and homework exercises were good and the discussion forums were helpful.

The course starts from the basics, but it is not an introductory course: it requires a decent mathematical background (linear algebra and calculus, Fourier transforms) and some previous knowledge of quantum physics is useful, although not strictly necessary.

The course is taught by two teachers with very different styles covering quite different material, this makes it a bit disjointed. I would have preferred a single teacher thorough the whole course.

The program covers various advanced subjects like the Feynman path integral, quantum localization, superconductivity, the Dirac equation, the time-dependent Schrodinger equation and so on. The problem is that having to start more or less from the basics it is impossible to cover any of these advanced subjects in any detail: they are just introduced but not explained enough to get a real understanding. I think it would be better to limit the course to the more ‘standard’ introductory stuff, with maybe just one or two advanced topics covered in more details.

Through a post in one of the Coursera’s forum I discovered the ‘Theoretical minimum’ Stanford lectures. It is not a course but just a collection of live lectures by the same professor. I watched a couple of them on quantum field theory, the first particle physics series and now I am watching the quantum entanglement ones.

The lectures are for a general audience, so they usually start from the very basics and then build the whole theory step by step. The teacher – professor Susskind – manages to introduce very complex stuff in an approachable way.

Being live lectures, there tend to be quite a lot of question from the audience, backtracking and interruptions – these makes them much longer that the recorded Coursera lectures, but I found that this makes the material more easily absorbed (by me at least).

They are very good introduction to modern physics for anyone with some scientific background and enough time to watch an entire series: watching individual lectures would probably not work.

I did not find anything ‘pre cooked’ to convert my billiard simulator to use special relativity, so I started to work on the necessary math from scratch – that turned out to be problematic (to say the least).

For example, the collision of a ball with the side of the table become something like this:

because – in the frame of reference of the table – the ball contracts along its direction of motion. This deformation affects the computation of the collision time and changes the dynamic of the collision. The same collision seen in the frame of reference of the ball becomes instead:

because now it is the distance from the ball to the side that contracts, changing the collision angle.

A ball-ball collision is even worse, being something like this (in the table frame of reference):

I found some pretty good notes on relativistic dynamics by an Oberlin college professor. The last chapter has a section about ‘hard sphere forces’ – e.g. billiard balls – that concludes that the whole idea of handling them within special relativity is ‘ludicrous’. At the same  time I found one article about rigid bodies in special relativity. I think that some kind of solution should be possible, but I have to study more.

For the moment being I am going through the relativistic dynamics notes (including the exercises!) and then I’ll have a look at the rigid body article.

I want to modify my billiard simulator to use special relativity instead than classical mechanics. The idea is to add the speed of light as a parameter, so that it is possible to see the relativistic effects increase as the speed of light decreases.

I need to derive the relativistic formulas for ball-side and ball-ball collisions, as well as adjusting the collision detection to account for relativistic time dilation and length contraction. I did not find anything much ‘ready made’ in various books and papers found on the Internet, so I’ll start working out the simple cases and (hopefully) build from that.

The simplest collision is ball-ball head-on:

It can be treated as one-dimensional, ignoring the vertical coordinate. The energy-momentum conservation equations are:

$$ \frac{E}{c^2} = m_0\gamma_0 + m_1\gamma_1 = m’_0\gamma’_0 + m’_1\gamma’_1 = \frac{E’}{c^2} $$

$$ p = m_0 \gamma_0 v_0 + m_1 \gamma_1 v_1 = m’_0 \gamma’_0 v’_0 + m’_1 \gamma’_1 v’_1 = p’$$

where \(m_0\) and \(m_1\) are the rest masses of the two balls and

$$ \gamma_0 = \frac{1}{\sqrt{1 - v_0^2 / c^2}},  \gamma_1 = \frac{1}{\sqrt{1 - v_1^2 / c^2}}$$

with \(c\) the speed of light. The primed symbols indicate the same quantities after the collision.

The collision can be not perfectly elastic, and so the rest masses of the balls can change, i.e. some of the kinetic energy gets converted into internal energy of the balls – hence \(m’_0\) and \(m’_1\) instead than just \(m_0\), \(m_1\) in the right-hand sides of the conservation equations.

These equations are difficult to solve for \( v’_0 \) and \( v’_1 \). The standard way to simplify them is to switch to the rest frame of the two-balls system, i.e. the frame of reference where the total momentum is zero. Indicating with \(\bar v_0\) and \(\bar v_1\) the ball velocities in this frame of reference, the momentum conservation equation becomes:

$$ m_0 \bar\gamma_0 \bar v_0 + m_1 \bar \gamma_1 \bar v_1 = 0 = m’_0 \bar\gamma’_0 \bar v’_0 + m’_1 \bar\gamma’_1 \bar v’_1 $$

(note that the rest masses are Loretz-invariant, so they stay the same). This equation has solutions

$$ m’_0 \bar\gamma’_0 \bar v’_0 = -\epsilon m_0 \bar\gamma_0 \bar v_0 $$

$$ m’_1 \bar\gamma’_1 \bar v’_1 = -\epsilon m_1 \bar\gamma_1 \bar v_1 $$

where \(\epsilon\) is a positive constant, and the minus sign indicates that the velocities change direction after the collision. Even with this simplification, the single conservation of energy equation is not enough to determine both the rest masses and the velocities after the collision, so we need to make some further assumption.  In the classical, perfectly elastic case the energy of each ball in the zero momentum frame of reference is conserved after the collision. Assuming that this is still the case in a relativistic collision gives the equations:

$$ m_0 \bar\gamma_0 = m’_0 \bar\gamma’_0 $$

$$ m_1 \bar\gamma_1 = m’_1 \bar\gamma’_1 $$

Solving for \( \bar v’_0 \), \(m’_0\) and \( \bar v’_1 \), \(m’_1\) gives:

$$ \bar v’_0 = –\epsilon \bar v_0 $$

$$ \bar v’_1 = -\epsilon \bar v_1 $$

$$ m’_0 = m_0 \frac{\bar\gamma_0}{\bar\gamma’_0} = m_0 \sqrt\frac{1-\epsilon^2 \bar v_0^2 / c^2}{1- \bar v_0^2 / c^2} $$

$$ m’_1 = m_1 \frac{\bar\gamma_1}{\bar\gamma’_1} = m_1 \sqrt\frac{1-\epsilon^2 \bar v_1^2 / c^2}{1- \bar v_1^2 / c^2} $$

The first two equations are exactly as in the classical case, with \(\epsilon\) being the coefficient of restitution.

To compute the final velocities in the frame of reference of the billiard table it is necessary to convert the initial velocities to the zero-momentum frame of reference, apply the formulas above and then convert back the results to the original frame of reference.

The Lorentz transformation of the momentum \(p\) of a body with energy \(E\) to a frame of reference with velocity \(u\) is:

$$ p’ = \gamma_u (p – \frac{u}{c^2}E) $$

so the velocity of the zero-momentum frame of reference is

$$ \begin{equation}\tag{U} u = \frac{p}{E/c^2} = \frac{m_0\gamma_0 v_0 + m_1 \gamma_1 v_1}{m_0 \gamma_0 + m_1 \gamma_1} \end{equation}$$

and the velocity transformations are:

$$ \bar v = \frac{v-u}{1-vu / c^2} $$

$$ v’ = \frac{\bar v’+u}{1+\bar v’u / c^2} $$

hence:

$$ \begin{equation}\tag{V} v’ = \frac{u(1-u v / c^2) + \epsilon (u – v)}{1- u v / c^2+\epsilon u / c^2 (u-v)} \end{equation} $$

$$ \begin{equation}\tag{M} m’ = m \sqrt{\frac{1 – 2(1-\epsilon^2)vu/c^2 – \epsilon^2 (v^2 + u^2) / c^2 + v^2u^2/c^4}{1- (v^2+u^2)/c^2 + v^2u^2/c^4}} \end{equation} $$

Classical limit

In the classical limit (U) becomes:

$$ u = \frac{m_0 v_0 + m_1 v_1}{m_0 + m_1} $$

and (V) becomes:

$$ v’ = (1 + \epsilon) u – \epsilon v$$

combining the two:

$$ v’ = \frac {1+\epsilon} {m_0+m_1} [m_0 (v_0 – v) + m_1 (v_1 – v)] + v$$

that gives:

$$ v’_0 = –m_1 \frac {1+\epsilon} {m_0 + m_1} (v_0 – v_1) + v_0 $$

$$ v’_0 = m_0 \frac {1+\epsilon} {m_0 + m_1} (v_0 – v_1) + v_1 $$

that is the same solution as in the two dimensional classical case when \(\alpha = 0\).

Inelastic collision

For a perfectly inelastic collision \(\epsilon = 0\) and so (V) becomes simply \( v’ = u \): the two balls have always the same velocity after the collision, i.e. they ‘stick together’ as expected.

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.

The velocity of the balls progressively decreases due to three different effects: loss of energy during collisions, friction with the table (rolling friction or rolling resistance) and friction with the air (drag).

The loss of energy during collisions is already accounted for by the coefficient of restitution used in the collisions formulas.

The general formula for the rolling friction is:

$$ F_{rr} = c_{rr}N $$

where \(F_{rr}\) is the resulting rolling friction force, \(c_{rr}\) is the coefficient of rolling friction – that depends on the materials of the ball and the table – and \(N\) is the downward force exerted by the ball on the table. \(N\) is equal to the mass \(m\) of the ball times the gravitational acceleration \(g\): \(N = mg\). The deceleration \(a_{rr}\) of the ball due to the rolling friction is the friction force divided by the mass of the ball: \(a_{rr}=F_{rr}/m\). Combining these two equation the rolling friction deceleration is:

$$ a_{rr} = c_{rr}g $$

The formula for the air drag is:

$$ F_d = \frac{1}{2}\rho_a v^2 c_d A $$

where \(F_d\) is the resulting drag force, \(\rho_a\) is the air density, \(v\) is the ball speed, \(c_d\) is the drag coefficient and \(A\) is the cross-sectional area of the ball. The deceleration \(a_{d}\) of the ball due to the air drag is the air drag force divided by the mass of the ball:

$$ a_d = \frac{F_d}{m} = \frac{F_d}{\frac{4}{3}\pi\rho_b r^3} $$

where \(\rho_b\) is the density of the ball and \(r\) is its radius. The cross-sectional area of the ball is \(A=\pi r^2\) – putting all this together the air drag deceleration is:

$$ a_d = f_d \frac{v^2}{r}$$

where the air drag factor \(f_d\) is defined as

$$ f_d = \frac{3\rho_a c_d}{8\rho_b} $$

The drag coefficient \(c_d\) for a sphere varies quite a bit with the velocity of the object – see here for example, but overall the air drag is quite small, so even using a constant would not change the result much; hence given the air density and the ball density \(f_d\) can be considered a constant as well.

Assuming that the time interval used to update the simulation is small compared with the rate of change of the ball velocities, the speed after the interval is simply the original speed decremented by the deceleration times the interval:

$$ v_t = v_0 – (a_{rr}+a_d)t $$

The resulting code is:

  /**
   * Updates the ball velocity, applying rolling resistance and air drag for a specified time interval
   * @param t the time interval to use. Only first-order effects are considered, so the time interval 
   * must be small compared to the rate of change of the velocity
   * @param airDragFactor total air drag factor computed from the ball and air densities and the 
   * air drag coefficient (assumed to be constant)
   * @param rollingResistanceDeceleration total rolling resistance deceleration computed from the 
   * rolling resistance coefficient and the gravitational acceleration 
   */
  updateVelocity(t: number, airDragFactor: number, rollingResistanceDeceleration: number) {
    var speed2 = this.v * this.v + this.w * this.w;
    if (speed2 > 0) {
      var airResistanceDeceleration = airDragFactor * speed2 / this.radius;
      var totalDeceleration = airResistanceDeceleration + rollingResistanceDeceleration;
      var speed = Math.sqrt(speed2);
      var newSpeed = speed - totalDeceleration * t;
      if (newSpeed <= 0) {
        // The ball stopped
        this.v = 0;
        this.w = 0;
      } else {
        // Update the speed, keeping the velocity direction the same (the air drag and rolling resistance 
        // forces are in the exact opposite direction of the velocity)
        this.v = this.v * newSpeed / speed;
        this.w = this.w * newSpeed / speed;
      }
    }
  } // updateVelocity

The last part covered the collisions between balls, what about collisions with the sides of the table? If the balls are not spinning and there is no ball-table friction the component of the ball velocity that is parallel to the table side remains the same, and the component perpendicular to the side changes direction and intensity based on the coefficient of restitution. Using primed letters for the velocity after the collisions, the formula for a collision with the left or right side is:

$$ v’ = –cv $$

$$ w’ = w $$

where \(c\) is the ball-table coefficient of restitution (that in general is different from the ball-ball one). The formula for a collision with the ‘top’ (as seen on the screen) or ‘bottom’ side is:

$$ v’ = v $$

$$ w’ = -cw $$

With this last piece, the complete update code is:

  update(dt: number) {
    while (dt > 0) {
      var firstCollisions = this.detectCollisions(dt, 0, 400, 0, 300);
      if (firstCollisions.t > 0) {
        // The balls move freely up to the time of the first collision: update their positions
        for (var i = 0; i < this.balls.length; i++) {
          var ball = this.balls[i];
          ball.updatePosition(firstCollisions.t);
        }
      }
      // Compute the new velocities after the collisions
      for (var i = 0; i < firstCollisions.collisions.length; i++) {
        var collision = firstCollisions.collisions[i];
        switch (collision.type) {
          case "x":
            collision.b1.v = -collision.b1.v * this.parameters.sideRestitution;
            break;
          case "y":
            collision.b1.w = -collision.b1.w * this.parameters.sideRestitution;
            break;
          case "b":
            collision.b1.collide(collision.b2, this.parameters.ballRestitution);
            break;
        }
      }
      // Continue with the remaining time
      dt -= firstCollisions.t;
    }
    this.draw();
  } // update

where ‘detectCollisions’ is a function that uses the code described previously to detect the first collision – or collisions – that occur within the specified interval. It returns an object containing the time of the collision and an array of the balls that collide. Each element of this array has two members ‘b1’ and ‘b2’ containing the colliding balls (‘b2’ is null for collisions with the sides) plus a ‘type’ member that indicates if the collision is with the left or right sides (“x”), with the top or bottom sides (“y”) or between two balls (“b”).

Here is the ‘detectCollisions’ code:

  private detectCollisions(dt: number, minx: number, maxx: number, miny: number, maxy: number) {
    var result = {
      t: dt,
      collisions: <{ type: string; b1: Ball; b2: Ball; }[]>[]
    }
    var addCollision = (t, collision) => {
      if (t === result.t) {
        // The new collision happens at the exact same time of the current one, it has to be added to the list
        result.collisions.push(collision);
      } else {
        // The new collision happens before the current one, so it replaces the entire list
        result.collisions = [collision];
        result.t = t;
      }
    }
    for (var i = 0; i < this.balls.length; i++) {
      // Collisions with the sides
      var ball = this.balls[i];
      var t = ball.sideXCollisionTime(minx, maxx);
      if (t && t <= result.t) {
        addCollision(t, { type: "x", b1: ball, b2: <Ball>null });
      }
      t = ball.sideYCollisionTime(miny, maxy);
      if (t && t <= result.t) {
        addCollision(t, { type: "y", b1: ball, b2: <Ball>null });
      }
      // Ball-ball collisions
      for (var j = i + 1; j < this.balls.length; j++) {
        var otherBall = this.balls[j];
        t = ball.collisionTime(otherBall);
        if (t && t <= result.t) {
          addCollision(t, { type: "b", b1: ball, b2: otherBall });
        }
      }
    }
    return result;
  } // detectCollisions

As promised in the title this code has a problem though: it works fine in most cases but fails when there are multiple ball-ball collisions at the same time.

Have a look at this: the ball coming from the left hits the two balls at rest at the same time, but the code evaluates the collisions one at a time: it first consider the collisions with the ball at the top, updating the velocity of the moving ball to move down, and then uses this velocity to compute the second collision. The result is clearly wrong: the moving ball should bounce straight back, and the result should not depend on the order of the balls.

I have no solution for this at the moment being.

..or more formally: elastic collision of spheres in two dimensions.

Let’s start with a simple case: ball 1 moving with velocity \((\bar v_1, \bar w_2)\) hitting ball 2 at rest. To make things even easier ball 1 and ball 2 are aligned horizontally at the moment of impact – i.e.  they have the same \(y\) coordinate:

In our model (at the moment being at least) balls are not spinning, and there is no ball-ball friction, so the force between balls in such a collision is only along the horizontal line connecting their centers. This means that the vertical components of the velocities before and after the collision are the same:

$$ \bar{w}’_1 = \bar{w} $$

$$ \bar{w}’_2 = 0 $$

- where the primed values are the ones after the collision. The horizontal components of the velocities on the other hand are related by the ball-ball coefficient of restitution \(c\):

$$ \bar{v}’_2 - \bar{v}’_1 = c(\bar{v}_1 - \bar{v}_2) = c\bar{v}_1  $$

(\(\bar{v}_2\) is zero). The horizontal momentum must be the same before and after the collision:

$$ m_1\bar{v}’_1 + m_2\bar{v}’_2 = m_1\bar{v}_1 $$

- where \(m_1\) and \(m_2\) are the masses of the two balls. Combining these two equation gives:

$$ \bar{v}’_1 = \frac{m_1 – cm_2}{m_1+m_2}\bar{v}_1 $$

$$ \bar{v}’_2 = \frac{(1+c)m_1}{m_1+m_2}\bar{v}_1 $$

that is the result for this simple case: the velocities after the collision in terms of the (known) velocities before the collision.

What abut the general case in which both balls move and are not aligned?

This second case can be converted in the simple one with a coordinate transformation: subtract the velocity of the second ball from all velocities – i.e. use the frame of reference of the second ball - and rotate the axis by the angle \(\alpha\):

$$ \bar{v} = (v-v_2)\cos\alpha + (w-w_2)\sin\alpha $$

$$ \bar{w} = -(v-v_2)\sin\alpha + (w-w_2)\cos\alpha $$

where \((v, w)\) is an arbitrary velocity and \((\bar{v}, \bar{w})\) is the corresponding transformed velocity. The angle \(\alpha\) is*:

$$ \alpha = \arctan \frac{y_1– y_2}{x_1 – x_2} $$

where \((x_1, y_1)\) and \((x_2, y_2)\) are the positions of the two balls.

The inverse transformation is:

$$ v = \bar{v}\cos\alpha - \bar{w}\sin\alpha + v_2 $$

$$ w = \bar{v}\sin\alpha + \bar{w}\cos\alpha + w_2 $$

Now the procedure to solve the general case is: transform the original velocities in their ‘bar’ equivalent; apply the formula of the simple case; apply the inverse transformation to go back to the velocities in the original coordinate system (without ‘bar’). These steps produce the following results:

$$ v’_1 = \frac{1}{m_1+m_2}\Big[\big(m_1-m_2(c\cos^2\alpha-\sin^2\alpha)\big)(v_1-v_2)-(1+c)m_2\sin\alpha\cos\alpha(w_1-w_2)\Big]+v_2$$

$$ w’_1 = \frac{1}{m_1+m_2}\Big[-(1+c)m_2\sin\alpha\cos\alpha(v_1-v_2) + \big(m_1-m_2(c\sin^2\alpha-\cos^2\alpha)\big)(w_1-w_2)\Big]+w_2$$

$$ v’_2 = \frac{(1+c)m_1}{m_1+m_2}\cos\alpha\Big[(v_1-v_2)\cos\alpha + (w_1-w_2)\sin\alpha\Big] + v_2 $$

$$ w’_2 = \frac{(1+c)m_1}{m_1+m_2}\sin\alpha\Big[(v_1-v_2)\cos\alpha + (w_1-w_2)\sin\alpha\Big] + w_2 $$

The first two equations are quite unwieldy, but adding and subtracting \(v_1\) from the first one and \(w_1\) from the second one, and then doing some algebra produces:

$$ v’_1 = –m_2a\cos\alpha+ v_1 $$

$$ w’_1 = –m_2a\sin\alpha+ w_1 $$

$$ v’_2 = m_1a\cos\alpha+ v_2 $$

$$ w’_2 = m_1a\sin\alpha+ w_2 $$

where

$$ a = \frac{1+c}{m_1+m_2}\Big[(v_1-v_2)\cos\alpha + (w_1-w_2)\sin\alpha\Big] $$

A nice symmetric solution – and correct to boot.

Assuming that all the balls have the same density, their mass is a fixed constant times their radius cubed. In the formulas above the masses appear always as a quotient, so they can be replaced with the balls’ cubed radiuses.

The final code is this:

  /**
   * Updates the velocities of this ball and another one after a collision
   * The coordinate of the balls must be at the collision point.
   * @param otherBall second colliding ball
   * @param restitution coefficient of restitution for a ball-ball collision
    */
  collide(otherBall: Ball, restitution: number) {
    var dx = this.x - otherBall.x;
    var dy = this.y - otherBall.y;
    var dv = this.v - otherBall.v;
    var dw = this.w - otherBall.w;
    var alpha = Math.atan2(dy, dx);
    var sinAlpha = Math.sin(alpha);
    var cosAlpha = Math.cos(alpha);
    var m1 = Math.pow(this.radius, 3);
    var m2 = Math.pow(otherBall.radius, 3);
    var a = (1 + restitution) / (m1 + m2) * (cosAlpha * dv + sinAlpha * dw);
    this.v = -m2 * a * cosAlpha + this.v;
    this.w = -m2 * a * sinAlpha + this.w;
    otherBall.v = m1 * a * cosAlpha + otherBall.v;
    otherBall.w = m1 * a * sinAlpha + otherBall.w;
  } // collide

(*) That ‘naïve’ formula for \(\alpha\) does not work when \(x_1-x_2\) is zero, nor it handles the various possible combination of signs of the numerator and denominator, but in the code all that is handled by Math.atan2().

I failed with the Udacity course (see previous episode), but I still want to improve my JavaScript / HTML5 / canvas skills. 

I am a software developer by trade, but a physicist by training: physics was my passion before computers sort-of-displaced it. I always had an itch to get back into it.

Hence the idea of writing physics simulation code running in a browser, and billiard seems a reasonably simple but interesting case.

The plan is to write a series of posts here describing the development of the simulation, writing about both the physics and the software parts. It will be a work-in-progress kind of thing: I don’t have a finished program that can be described from A to Z; I am going to describe the various pieces as I implement (and possibly ditch) them.

The ambition is to start with classical mechanics and then ‘graduate’ to relativity and/or quantum mechanics. It would be interesting to have a relativistic simulation and be able to play with \(c\) and see how the classical behavior changes into the relativistic one.

The simulation will run entirely in the browser, without a server part. The code will be in TypeScript, compiled into JavaScript for deployment.

Why TypeScript instead than plain JavaScript? I always used strongly typed languages (C#, C++, Pascal, C), and JavaScript seems a bit too loose to me, especially for bigger projects. TypeScript seems a very good solution: offers strong typing, classes and modules, but it is very similar to JavaScript and can use native JavaScript libraries.

The snippets of code I am going to place here will most likely run unchanged – or with very minor changes - as pure JavaScript, so whoever is interested in the simulation and physics part won’t have to learn TypeScript just to follow the implementation.

The simulation I developed so far can be seen here, the source code is on GitHub here. Everything is still very very rough – hopefully it will improve.