mimoSite Programming and some other things

For the billiard simulation user interface I wanted to select a value with a slider, displaying it at the same time – so I created a new jQuery UI custom widget based on the standard slider one. It is used to edit the coefficients of restitution and the rolling resistance coefficient here.  The Javascript source code is here, and the TypeScript definition file is here.

The widget is pretty basic: it supports only one value, the slider must be horizontal and the value is always displayed on the right of the slider.

Some problem I had writing the widget:

Case of the widget name

I originally called the widget ‘sliderWithValue’ using camel case. The name of the widget is used by default as the prefix for the complete name of the widget events – but converted to lowercase – so the name of the change event is ‘sliderwithvaluechange’ and NOT ‘sliderWithValuechange’. The documentation does not mention this conversion to lowercase.

I renamed the widget ‘valueslider’ to avoid this problem.

Setting the width of a ‘span’ element.

To space correctly the label containing the value I needed to convert a width from ems to pixels, to do this the only way I found is to create a dummy HTML element, setting its widths in ems and then getting its pixel width. The problem is that using Internet Explorer 10 this trick works fine using ‘span’ element, whereas in Chrome and Safari attempts to set the width of  this empty ‘span’ elements are ignored – it is necessary to use a ‘div’.

Events and widget ‘destroy’ method

The jQuery UI documentation instructs to wire events either directly in the widget options:

$( ".selector" ).slider({
  change: function( event, ui ) {}
});

or using the ‘on’ function:

$( ".selector" ).on( "slidechange", function( event, ui ) {} );

but they are not exactly equivalent: event handlers defined using ‘on’ are bot removed by the ‘destroy’ method of the widget, so if the code creates and destroys widgets dynamically it has to take care of removing the handlers defined using ‘on’.

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.

Another way to make the simulation a bit more realistic is to produce some sound when balls hit each other and when they hit the side of the table.

Playing a sound in the browser turns out to be super-easy – create ‘Audio’ objects loading sound files:

  private audioBallBall = new Audio("sounds/ball-ball.mp3");
  private audioBallSide = new Audio("sounds/ball-side.mp3");

and play them in the loop that process the collisions:

      // 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":
            this.audioBallSide.volume = Math.min(this.maxSpeed, Math.abs(collision.b1.v)) / this.maxSpeed;
            this.audioBallSide.play();
            collision.b1.v = -collision.b1.v * this.parameters.sideRestitution;
            break;
          case "y":
            this.audioBallSide.volume = Math.min(this.maxSpeed, Math.abs(collision.b1.w)) / this.maxSpeed;
            this.audioBallSide.play();
            collision.b1.w = -collision.b1.w * this.parameters.sideRestitution;
            break;
          case "b":
            var dx = (collision.b1.x - collision.b2.x);
            var dy = (collision.b1.y - collision.b2.y);
            var relativeSpeed = Math.abs((collision.b1.v - collision.b2.v) * dx + (collision.b1.w - collision.b2.w) * dy) / Math.sqrt(dx*dx + dy*dy);
            this.audioBallBall.volume = Math.min(this.maxSpeed, relativeSpeed) / this.maxSpeed;
            this.audioBallBall.play();
            collision.b1.collide(collision.b2, this.parameters.ballRestitution);
            break;
        }
      }

The volume is adjusted to decrease as the collision speed decreases, with 100% volume corresponding to some (rather arbitrary) maximum speed. I did not investigate the actual physics – I suspect that the sound volume does not simply decrease linearly with the speed, but it sounds good enough.

The audio files are loaded asynchronously by the browser, so in theory they should be played only after receiving a notification that loading completed (see here for example), but being very small file this is not really necessary.

I recorded myself the two different sounds used for ball-ball and ball-side collisions – hitting two different surfaces with a pen, not very realistic but good enough for now (probably the title of this post should be ‘good enough’). I used Audacity to record the sound and cut them – I didn’t do any other editing.

There are libraries of sound effects that can be bought quite cheaply, including billiard sound – see this for example, but they appear to be more complex sounds, not suitable for just a single impact.

UPDATE

The code above contained a bug in the computation of the relative speed for ball-ball collisions, now fixed.

The simulation is two-dimensional: the table and balls are seen from straight above. Drawing the balls as simple circles produces a very un-realistic effect though; to improve things it is necessary to render in some way the rotation of the balls. The easiest way to do this is to add some feature to a ball surface – like a white circle – and change its position and size as the ball rotates.

A point on the surface of a ball is identified by two angles \((\varphi, \vartheta)\):

As the balls moves, it rotates by an angle \(\delta = \Delta s / r\), where \(r\) is the ball radius and \(\Delta s = \sqrt{\Delta x^2 + \Delta y^2}\) is the distance travelled. The direction of the movement is the direction of the ball velocity, so this is the situation as seen ‘from above’:

where \(\alpha = \arctan(w/v)\). Let’s rotate the \(x, y\) axes by the angle \(\alpha\), so that the movement is along the new \(\bar x\) axis: to do it simply subtract \(\alpha\) from \(\varphi\). Using these new coordinates the situation as seen ‘from the side’ is:

(the \(\bar y\) axis goes away from you ‘through’ the screen). The initial coordinate of the point are:

$$ \bar x = r \sin \vartheta cos (\varphi-\alpha) $$

$$ \bar y = r \sin \vartheta sin (\varphi-\alpha) $$

$$ z = r \cos \vartheta $$

after the ball moves the \(\bar y\) coordinate is the same, whereas \(\bar x\) and \(z\) are rotated by \(\delta\) clockwise:

$$ \bar{x}’ = \bar x \cos\delta + z \sin\delta = r \sin \vartheta cos (\varphi-\alpha) \cos\delta + r \cos \vartheta \sin\delta $$

$$ \bar y’ = r \sin \vartheta sin (\varphi-\alpha) $$

$$ z’ = –\bar x \sin\delta + z \cos\delta = –r \sin \vartheta cos (\varphi-\alpha) \sin\delta + r \cos \vartheta \cos\delta $$

The angles after the ball moves are then:

$$ \varphi’ = \alpha + \arctan \frac{\bar y’}{\bar x’} = \alpha + \arctan \frac {\sin \vartheta sin (\varphi-\alpha)} {\sin \vartheta cos (\varphi-\alpha) \cos\delta + \cos \vartheta \sin\delta} $$

$$ \vartheta’ = \arccos \frac{z’}{r} = \arccos(-\sin \vartheta cos (\varphi-\alpha) \sin\delta + \cos \vartheta \cos\delta) $$

and here is the corresponding code:

  /**
   * Updates the ball position and orientation, applying the current velocity for a specified time interval
   * @param t the time interval to use. The velocity is considered constant, so the time interval must be 
   * small compared to the rate of change of the velocity
   */
  updatePosition(t: number) {
    var dx = this.v * t;
    var dy = this.w * t
    // Update the ball position
    this.x += dx;
    this.y += dy;
    var ds2 = dx * dx + dy * dy;
    if (ds2 > 0) {
      // Update the ball orientation
      var delta = Math.sqrt(ds2) / this.radius;
      var sinDelta = Math.sin(delta);
      var cosDelta = Math.cos(delta);
      var alpha = Math.atan2(this.w, this.v);
      var sinTheta = Math.sin(this.theta);
      var cosTheta = Math.cos(this.theta);
      var phiAlpha = this.phi - alpha;
      var sinPhiAlpha = Math.sin(phiAlpha);
      var cosPhiAlpha = Math.cos(phiAlpha);
      this.phi = alpha + Math.atan2(sinTheta * sinPhiAlpha, sinTheta * cosPhiAlpha * cosDelta + cosTheta * sinDelta);
      this.theta = Math.acos(-sinTheta * cosPhiAlpha * sinDelta + cosTheta * cosDelta);
    }
  } // updatePosition

A circle drawn on the surface of the ball that is centered at the point \((\varphi, \vartheta)\) looks like a circle when \(\vartheta=0\), and then get progressively ‘squashed’ as \(\vartheta\) increases, until is disappears when \(\vartheta>\pi/2\):

The circle is actually a spherical one, that should be projected on a horizontal plane to compute the exact shape to be drawn; approximately though the circle becomes an ellipse as \(\vartheta\) increases – with the axis in the direction from its center to the center of the ball decreasing with \(\cos \vartheta\). Using this approximation the code to draw the ball is this:

  /**
   * draw: draws the ball
   * @param ctx the canvas rendering context to use to draw the ball
   */
  draw(ctx: CanvasRenderingContext2D) {
    ctx.save();
    // Move the coordinates to the center of the ball - simplifies everything else
    ctx.translate(this.x, this.y);
    // Gradient from the ball color to black, used to shade the ball color to give an illusion of depth
    var ballColorGradient = ctx.createRadialGradient(0, 0, 0, 0, 0, this.radius*3);
    ballColorGradient.addColorStop(0, this.color);
    ballColorGradient.addColorStop(1, "black");
    // Gradient from white to black, used to shade the white circle to give an illusion of depth
    var whiteGradient = ctx.createRadialGradient(0, 0, 0, 0, 0, this.radius * 3);
    whiteGradient.addColorStop(0, "white");
    whiteGradient.addColorStop(1, "black");
    // Draw the ball: a circle filled with the ball color shading to black
    ctx.fillStyle = ballColorGradient
    ctx.beginPath();
    ctx.arc(0, 0, this.radius, 0, Math.PI*2);
    ctx.fill();
    if (this.theta < Math.PI / 2) {
      // Draw the white circle if it is visible
      var d = this.radius * Math.sin(this.theta);
      var cosTheta = Math.cos(this.theta);
      var s = this.circleRadius * cosTheta;
      if (d - s < this.radius) {
        var cosPhi = Math.cos(this.phi);
        var sinPhi = Math.sin(this.phi);
        // Clip to the ball's circle - do not want to draw parts of the white circle that fall outside the ball borders
        ctx.clip();
        // Move the coordinates to the center of the white circle
        ctx.translate(d * cosPhi, d * sinPhi);
        // Compress the coordinates by cosTheta in the direction between the center of the white circle and the center of the ball
        ctx.rotate(this.phi);
        ctx.scale(cosTheta, 1);
        // Draw the white circle
        ctx.fillStyle = whiteGradient;
        ctx.beginPath();
        ctx.arc(0, 0, this.circleRadius, 0, 2 * Math.PI);
        ctx.fill();
      }
    }
    ctx.restore();
  } // draw

Note how the code uses gradients to fill the ball and the white circle to give a more ‘tridimensional’ look.

Aside: these first experiences of using the HTML canvas are very positive: the API is complete, easy to use and works well in multiple ‘modern’ browsers

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().

First of all, what exactly ‘detect collisions’ means? As explained in the previous post the update algorithm needs the time of the next collisions between two balls and between a ball and the sides of the table – so ‘detecting collisions’ means computing these times for all the balls.

The initial position of the center of the \(k\)-th ball is \((x_k, y_k)\) and its velocity is \((v_k, w_k)\). The position of the center of each ball varies with time, at time \(t\) it is \((x_k+v_kt, y_k+w_kt)\).

A ball collides with one side of the table when the distance between its center and the side equals the ball radius. A ball moves toward the left or right sides with its horizontal speed \(v_k\) – if this speed is \(0\) the ball won’t collide, otherwise the collision times are \((x_{min}+r_k-x_k)/v_k\) and \((x_{max}-r_k-x_k)/v_k\), where \(x_{min}\) and \(x_{max}\) are the horizontal coordinate of the left and right sides (and hence the minimum and maximum values for the \(x\) coordinates). Only one of these values will be positive, and that is the result (if the ball moves towards the right \(v_k\) is positive and it will collide with the right side in the future – collision time positive, and had collided with the left side in the past – collision time negative; vice-versa for a ball moving to the left).

The collisions time with the top and bottom sides are computed in the same way – replacing \(x_k\) with \(y_k\), \(v_k\) with \(w_k\), \(x_{min}\) with \(y_{min}\) and \(x_{max}\) with \(y_{max}\). 

The complete code for collisions between a ball and the table sides is:

  /**
   * Returns the time of the first future collision with a side of the table. 
   * Returns null if there is no collision.
   * @param x ball horizontal or vertical position
   * @param v ball horizontal or vertical speed
   * @param min minumum horizontal or vertical position (i.e. co-ordinate of one side of the table)
   * @param max maximum horizontal or vertical position (i.e. co-ordinate of the other side of the table)
   */
  private sideCollisionTime(x: number, v: number, min: number, max: number) {
    if (v === 0) {
      // Not moving towards the sides: no collision
      return null;
    }
    var result;
    if (v < 0) {
      // Moving up/left - check agains the minumum 
      result = (min - x + this.radius) / v;
    } else {
      // Moving down/right - check agains the maximum
      result = (max - x - this.radius) / v;
    }
    if (result >= 0) {
      return result;
    }
    return null;
  } // sideCollisionTime

  /**
   * Returns the time of the first future collision with the left or right sides of the table. 
   * Returns null if there is no collision.
   * @param min minumum horizontal position (i.e. x co-ordinate of the left side of the table)
   * @param max maximum horizontal position (i.e. x co-ordinate of the right side of the table)
   */
  sideXCollisionTime(min: number, max: number) {
    return this.sideCollisionTime(this.x, this.v, min, max);
  } // sideXCollisionTime

  /**
   * Returns the time of the first future collision with the top or bottom sides of the table. 
   * Returns null if there is no collision.
   * @param min minumum vertical position (i.e. y co-ordinate of the top side of the table)
   * @param max maximum vertical position (i.e. x co-ordinate of the bottom side of the table)
   */
  sideYCollisionTime(min: number, max: number) {
    return this.sideCollisionTime(this.y, this.w, min, max);
  } // sideYCollisionTime

What about collisions between balls? This is a bit more complex – two balls collide when their distance is equal to the sum of their radiuses:

The square of the distance between balls \(k\) and \(l\) at time \(t\) is:

$$s^2(t)=(x_k + v_kt - x_l - v_lt)^2 + (y_k + w_kt - y_l - w_lt)^2 = (\Delta x + \Delta vt)^2 + (\Delta y + \Delta wt)^2$$

where \(\Delta v = v_k – v_l\), \(\Delta w = w_k – w_l\), \(\Delta x = x_k – x_l\) and \(\Delta y = y_k – y_l\). Equating this to the squared sum of their radiuses produces:

$$(r_k + r_l)^2 = (\Delta x + \Delta vt)^2 + (\Delta y + \Delta wt)^2$$

that is a quadratic equation in \(t\). Solving it gives the collision time:

$$t_{kl} = \frac{-p \pm \sqrt{q}}{\Delta v^2 + \Delta w^2}$$

where \(p=\Delta v\Delta x + \Delta w\Delta y\) and  \(q = (\Delta v^2+\Delta w^2)(r_k+r_l)^2-(\Delta v\Delta y - \Delta w\Delta x)^2\).

If \(q\) is negative the equation has no (real) solutions, that means that the balls do not collide. If \(q\) is zero the balls have only one collision: they ‘glance’ each other. If \(q\) is greater than zero there are two ‘collisions’: the balls touch, then go through each other and touch again on the other side. In the complete simulation after the first collision the balls change course, so we are interested only in the first – i.e. minimum – value. Finally, one or both solutions can be negative – corresponding to collisions that happened in the past; we want only the collisions in the future, so the result is the minimum non-negative solution.

With this the algorithm is complete, but there is still a problem: it all works fine if the balls never overlap, but if they do we will detect spurious collisions.

If the balls do not overlap to begin with, and the simulation detect the collisions and make the balls bounce away when they touch, how is it possible to even have overlaps? The answer is floating-point arithmetic: the actual code uses floating-point values, that introduces rounding errors in the computation.

What can happen (and actually DO happen) is:

1. two balls move toward each other,
2. they reach the collision point,
3. the simulation updates their velocity so that they ‘bounce’,
4. due to rounding errors there is a slight overlap in their position,
5. this cause a new immediate (spurious) collision to be detected, with a new update in the balls velocities

…and so on. The final effect is that sometimes balls touch and then stick together, or start to bounce inside each other – sort of fun, but definitely not billiard.

The solution for this is to check the value of \(p\): only if it is negative the balls are moving toward each other, otherwise are moving apart (or parallel); so when \(p\) is not negative we always return ‘no collision’, this excludes the spurious collisions at step (5) above and fixes the simulation.

The complete code for collisions between two balls is:

  /** 
   * Returns the time of the first future collision between this ball and another ball. 
  *  Returns null if the balls do not collide
   * @param other the other ball
  */
  collisionTime(other: Ball) {
    var dv = this.v - other.v;
    var dw = this.w - other.w;
    var dx = this.x - other.x;
    var dy = this.y - other.y;
    var p = (dv * dx + dw * dy);
    if (p >= 0) {
      // The balls are moving apart or parallel
      return null;
    }
    var dv2 = dv * dv + dw * dw;
    var r = this.radius + other.radius;
    var q = dv2 * r * r - Math.pow(dv * dy - dw * dx, 2);
    if (q < 0) {
      // No real solution: the balls do not collide
      return null;
    }
    q = Math.sqrt(q);
    if (p <= q) {
      return (-p - q) / dv2;
    }
    return (-p + q) / dv2;
  } // collisionTime

How comes that the sign of \(p\) indicates if the balls are moving apart or not? The time derivative of the distance between the balls is

$$\frac{ds}{dt}\bigg|_{t=0} = \frac{\Delta x\Delta v + \Delta y\Delta w}{s} = \frac{p}{s}$$

\(s\) is positive, so if \(p\) is negative the derivative is negative and the distance is decreasing – i.e. the balls are moving towards each other. 

Another way to prove it is noticing that \(s\) is the scalar product of the relative position and speed of the balls, hence the sign of \(p\) is the sign of the cosine of the angle \(\varphi\) between those two vectors. If \(p\) is negative the angle is greater than \(90^\circ\), that means that the direction of the relative velocity is opposite of the direction of the relative position – i.e. the balls are moving towards each other.

[Calculus agrees with geometry – that is good]

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.

In the previous post I dealt with overloading constructors. The exact same system can be used to overload class methods - this code:

class DateHour {

  private date: Date;
  private relativeHour: number;

  . . .

  init(year: number, month: number, day: number, relativeHour: number);
  init(date: Date, relativeHour: number);
  init(dateOrYear: any, monthOrRelativeHour: number, day?: number, relativeHour?: number) {
    if (typeof dateOrYear === "number" && monthOrRelativeHour && day) {
      this.date = new Date(dateOrYear, monthOrRelativeHour, day);
      this.relativeHour = relativeHour;
    } else {
      var date = <Date> dateOrYear;
      this.date = new Date(date.getFullYear(), date.getMonth(), date.getDate());
      this.relativeHour = monthOrRelativeHour;
    }
  }
}

defines an 'init' method with two overloads, one taking a year/moth/day triplet and the other a date object. The one and only implementation is hidden by the overloads, and checks the type of the parameters to determine which of the two overloads to implement.

Finally, the same system can be used to overload stand-alone functions:

function initDateHour(year: number, month: number, day: number, relativeHour: number);
function initDateHour(date: Date, relativeHour: number);
function initDateHour(dateOrYear: any, monthOrRelativeHour: number, day ? : number, relativeHour?: number) {
  if (typeof dateOrYear === "number" && monthOrRelativeHour && day) {
    // Handle the year/moth/day case    
  } else {
    // Handle the date object case
  }
}

Part 1