Billiard simulation–part 6: balls rotation

June 2, 2013 at 11:09 AMMichele Mottini

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) {;
    // 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.arc(0, 0, this.radius, 0, Math.PI*2);
    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
        // 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.scale(cosTheta, 1);
        // Draw the white circle
        ctx.fillStyle = whiteGradient;
        ctx.arc(0, 0, this.circleRadius, 0, 2 * Math.PI);
  } // 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

Posted in: Programming

Tags: , , ,

Billiard simulation – part 2: detect collisions

May 4, 2013 at 11:07 AMMichele Mottini

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]

Posted in: Programming

Tags: , , ,

Billiard simulation - part 1: model and animation

April 28, 2013 at 6:25 PMMichele Mottini

The simulation model is pretty straightforward: the table is just a rectangle, represented by the coordinates of its top left corner and its width and height; the balls are spheres, represented by the position of their center \(\vec p_k=(x_k, y_k)\), their velocity \(\vec v_k=(v_k, w_k)\) and their radius \(r_k\) – where \(k\) is the index identifying the ball.

Billiard table with balls

I kept the radius as a property of each ball instead than a global parameter, to be able to model balls of different size.

This is just the geometrical part of the model; I’ll introduce additional parameters for the physics part (mass, friction, restitution) as I develop the physics model.

The animation logic is simple as well: the screen needs to be updated with a certain frequency – let’s say 30 times a second to provide a smooth display. To do this a timer calls an update function every \(\Delta t = 1/30\) seconds, and this function computes the new positions of all the balls and then re-draws everything.

The update function has to compute the change of the balls positions in the \(\Delta t\) seconds elapsed since the last update. For a ball moving freely this change is simply its velocity times \(\Delta t\). The balls do not move freely though: they can hit other balls or the sides of the table, and every time there is a collision the velocity of the colliding balls(s) changes. The complete update algorithm then is:

  1. Compute the time \(t_{min}\) of the first collision(s)
  2. If \(t_{min}\)  is equal or greater than \(\Delta t\) then move all the balls freely for \(\Delta t\) seconds and we are done (no collisions to consider)
  3. Otherwise: move all the balls freely for \(t_{min}\) seconds
  4. Compute the new velocities of all the balls that collide
  5. Repeat from (1), reducing \(\Delta t \) to  \(\Delta t – t_{min}\)

Next time: detecting collisions.

Posted in: Programming

Tags: , ,