Phrase of the week

June 30, 2013 at 10:23 AMMichele Mottini

He has a great work ethic when he is in the mood

My wife

Posted in: Opinion

Tags:

Billiard simulation–part 8: double collision

June 30, 2013 at 4:47 AMMichele Mottini

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:

doublecollision

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.

Posted in: Physics | Programming

Tags: , , , ,

Phrase of the week

June 23, 2013 at 12:22 PMMichele Mottini

A Book on Those Geometric Constructions Which Are Necessary for a Craftsman

A Book on What Is Necessary from the Science of Arithmetic for Scribes and Businessmen

Titles of X century mathematical texts by Abū al-Wafā' Būzjānī 

Posted in: Opinion

Tags:

Phrase of the week

June 16, 2013 at 2:31 PMMichele Mottini

And so I dream of a day when the only people who suffer from money-losing investments are the money-losing investors, whose only penance is lost money.

Ryan Avent of The Economist in 'The wages of sin'

Posted in: Opinion

Tags:

Tolkien on the NSA wiretapping

June 16, 2013 at 2:22 PMMichele Mottini

[Sam]  ‘I wish you’d take his Ring. You’d put things to right. You’d stop them digging up the gaffer and turning him adrift. You’d make some folk pay for their dirty work.’

‘I would,’ she [Galadriel] said. ‘That is how it would begin. But it would not stop with that’

‘The Lord of the Rings’, book two, chapter VII.

Posted in: Opinion

Tags:

Billiard simulation–part 7: sound

June 15, 2013 at 11:18 AMMichele Mottini

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.

Posted in: Programming

Tags: , ,

Phrase of the week

June 9, 2013 at 7:21 AMMichele Mottini

...activist investors seeking to crack Japan have fared only slightly better than the Christian missionaries of the 16th century (who were crucified).

The Economist in 'Goodbye, Mr Bond?'

Posted in: Opinion

Tags:

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)\):

thetaphi

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’:

fromabove

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:

fromside

(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\):

circle

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

Posted in: Programming

Tags: , , ,

Phrase of the week

June 2, 2013 at 9:50 AMMichele Mottini

The resource that is in shortest supply is usually time, since there is no way to create more of it.

Raymon Chen in Microspeak: booked

Posted in: Opinion

Tags: