(December 2013)

Click/tap the button above to pull the string at its middle. I have deliberately set the wave's rate of energy loss quite low, so by clicking the button repeatedly, you will be able to see the individual waves you are creating interfere with each other.

Can you can hit the button with the proper frequency to maximize the wave amplitude? That is, can you find the string's resonance frequency?

Try to make the waves so big they reach the "ceiling" :‑)

Note that due to the algorithm I am using, the number of waves doesn't impact the simulation speed. You may find the theory behind this implementation interesting - I certainly enjoyed figuring it out.

I based this Javascript port on my original, libSDL-based GPL implementation, which compiles and runs under all major OSes. As for the JS code, I tested it under Chrome/Firefox/Opera on a puny Atom 330 running ArchLinux, and also on an iPhone 3GS's Safari. It should therefore work fine on all modern browsers (that support the HTML5 canvas API) - if you are on an older browser, it's time to upgrade...

Click with the mouse (or tap, if you use a phone / tablet) anywhere you want in the blue pool above. You will see waves generated and reflected around the borders, interfering with themselves and any other waves you generate. Don't click only once - click many times, and marvel at the waves' interference patterns :‑)

Note that due to the algorithm I am using, the number of waves doesn't impact the simulation speed. You may find the theory behind this implementation interesting - I certainly enjoyed figuring it out.

I based this Javascript port on my original, libSDL-based GPL implementation, which compiles and runs under all major OSes. As for the JS code, I tested it under Chrome/Firefox/Opera on a puny Atom 330 running ArchLinux, and also on an iPhone 3GS's Safari. It should therefore work fine on all modern browsers (that support the HTML5 canvas API) - if you are on an older browser, it's time to upgrade...

## Theory

Assume we have a wave, $f_n(i)$, where $i$ refers to the horizontal screen axis (i.e. the X-pixel coordinate), and $f_n(i)$ is the Y-coordinate of that pixel, in frame $n$. We want to find what that Y-coordinate will change to, in frame $n+1$.

Consider the individual pixels to represent the water molecules, and assume that each molecule is influenced only by its two neighbours (left and right), as if they are connected to it with springs. That is, assume that a given pixel's height $f_n(i)$, is only influenced by the attractive force of its two neighbours:

• the pixel to the left:  $f_{n,L} = f_n(i-1)$
• the pixel to the right: $f_{n,R} = f_n(i+1)$.

With $f$ being the position, $v$ the velocity and $a$ the acceleration, we write down the laws governing motion, for $dt=1$ (yes, I know, dt is supposed to be infinitesimal ; it doesn't matter, we are engineers playing with computers, we make our own worlds ; keep reading :‑) Basically, we will be applying simple Euler integration with a huge $dt$... let's see what happens:

First, the pixel's position in the new frame is equal to the position in this frame, plus the vertical velocity times dt:

$$f_{n+1} = f_n + v_n \label{eq:f}$$

The velocity is equal to the old velocity plus the old acceleration times dt:

$$v_n = v_{n-1} + a_{n-1} \label{eq:v}$$

And since springs cause acceleration that is linearly related to the "pulled" distance:

$$a_{n-1} = k (f_{n-2,L} - f_{n-2}) + k (f_{n-2,R} - f_{n-2}) \label{eq:a}$$

...where $k$ is the spring's coefficient - the larger it is, the more acceleration is caused by the spring.

If you prefer the full expanded notation for pixel $i$:

$$a_{n-1}(i) = k (f_{n-2}(i-1) - f_{n-2}(i)) + k (f_{n-2}(i+1) - f_{n-2}(i)) \nonumber$$

These equations can be used to perform the simulation, and they would work fine. To achieve even greater simulation speeds, however - especially for the 2D simulation of water waves - we'll violate the universe's laws some more: we will use the acceleration at time $n+2$ (instead of $n-1$), which basically means that we will approximate the current acceleration rate with the one we would have in three frames' time. We will also set $k$ to 1, to further simplify the calculation.

Replacing \eqref{eq:v} in \eqref{eq:f}, using $a_{n+2}$ instead of $a_{n-1}$, and then applying \eqref{eq:a} leads to the following:

$$\begin{eqnarray} f_{n+1} & = & f_n + v_n\nonumber \\ & = & f_n + v_{n-1} + a_{n-1}\nonumber \\ & = & f_n + (f_n - f_{n-1}) + a_{n-1}\nonumber \\ & = & f_n + (f_n - f_{n-1}) + a_{n+2}\nonumber \\ & = & 2 f_n - f_{n-1} + a_{n+2}\nonumber\\ & = & 2 f_n - f_{n-1} + f_{n,L} - f_n + f_{n,R} - f_n\nonumber\\ & = & f_{n,L} + f_{n,R} - f_{n-1} \end{eqnarray}$$

In our fictitious universe, simulating a 1D-wave is amazingly simple: we add the $y$-coordinates of the two neighbouring pixels, and subtract the old value we used to have in the previous frame! Since the springs are losing power with each iteration, we will also scale the output of the last equation so that it diminishes over time (e.g. multiply by a factor of 0.99): $$$$\boxed{ f_{n+1}(i) = 0.99 (f_n(i-1) + f_n(i+1) - f_{n-1}(i)) \nonumber }$$$$

For the corresponding two-dimensional problem (the 2D simulation of water waves), we just average the effects of the X- and Y- coordinate waves:

$$\begin{eqnarray} f_{n+1}(i,j) & = & \frac{1}{2} & (f_{n+1,HorizontalWave} + f_{n+1,VerticalWave}) \nonumber \\ & = & \frac{1}{2} (& 0.99 (f_n(i-1,j) + f_n(i+1,j) - f_{n-1}(i,j)) + \nonumber \\ & & & 0.99 (f_n(i,j-1) + f_n(i,j+1) - f_{n-1}(i,j)))\nonumber \end{eqnarray}$$

So, for real-time simulation of two-dimensional waves, we only need this equation: $$$$\boxed{ f_{n+1}(i,j) = 0.99 ( \frac{f_n(i-1,j) + f_n(i+1,j) + f_n(i,j-1) + f_n(i,j+1)}{2} - f_{n-1}(i,j) ) }$$$$ Which translates to: we add the 4 neighbouring pixels, divide by 2, subtract the previous frame value, and scale by 0.99. Simple!

In the past, almost all of my graphics hacks were coded in an OS-agnostic manner, via libSDL and autoconf/automake. I could therefore execute them under Linux, Windows and OS/X, and rightfully claim that my code works everywhere.

There was a catch, though - to compile and run them, you... well, you had to be a person that loves coding, like me. I could provide you with pre-compiled binaries (and I do) ; but then you'd have to trust me that I am not some malevolent Greek, inventing Trojan Horses for fun and profit Timeo Danaos et dona ferentes, or, in my native language: "Φοβού τους Δαναούς και δώρα φέροντας"  :‑)

But now all that is in the past. My libSDL days are over ; when I feel the urge to hack stuff for my blog, I'll do it in the planet's portable assembly language from now on... that is, in Javascript :‑)

If my usage of the canvas API is incorrent or I am missing something else, please feel free to comment - below you will find the complete simulation code, for your easy perusal.

Enjoy!

//////////////////////////////////////
// Global variables
//////////////////////////////////////
//
// This is coded just for fun ; no, I don't use globals in production code.
// Honest! :-)

// The working image buffer dimensions.
// I use smallish resolution for the HTML5 canvas native resolution,
// so that the code will work even on puny mobile phones.
// (i.e. I don't want the code to require beastly CPUs)
g_width = 320;
g_height = 240;

// The global setInterval timer that calls the 'simulate' function.
g_timer = -1;

// Which tab is currently visible?
g_activeTab = 1;

// Has MathJax been configured in the past?
g_mathJax_configured = false;

// HTML5 canvas variables, shared between the two simulations.
g_dc = null;
g_canvas = null;
g_imageData = null;

// Index of the currently active buffer in our double-buffer scheme
g_1dwave_idx = 0;
g_2dwave_idx = 0;

// The two double-buffers of data - 1d and 2d
g_1dwave_y = [ new Float32Array(g_width), new Float32Array(g_width)];
g_2dwave_y = [ new Float32Array(g_width*g_height), new Float32Array(g_width*g_height)];

// Lookup table for the 2d simulation.
g_2d_offsets = new Int32Array(g_height);
for(var ii=0; ii<g_height; ii++) {
g_2d_offsets[ii] = g_width*ii;
}

// The backlog of pixel(s) to clear when we pull the string
g_pullBacklog = [];

// Function that checks whether we have a working HTML5 canvas
// and then sets up the canvas related globals and the timer.
function setupCanvasAndStartSimulation() {
var is1Dsim = g_activeTab === 0;
g_canvas = document.getElementById(is1Dsim?"1dwaveDC":"2dwaveDC");
if (!g_canvas.getContext) {
alert("You are using an old browser - please use a recent Firefox / Chrome / Opera / Safari");
return false;
}
g_dc = g_canvas.getContext("2d");
if (!is1Dsim) {
g_canvas.removeEventListener('mousedown', rainDrop);
}
// ...and start the simulation!
if (g_timer === -1) {
g_timer = setInterval(function() { simulate(); }, 20);
}
return true;
}

////////////////////////////////////////////////////////////////////
// 1D wave simulation
////////////////////////////////////////////////////////////////////

// Quickly plot two green or black pixels in the 1D canvas image data
function qplot(x, y, color) {
if (y<0 || x<0 || y>=g_height-1 || x>=g_width)
return;
var ofs = Math.floor(y)*g_width*4 + Math.floor(x)*4;
g_imageData[ofs + 1] = color;
g_imageData[ofs + 3] = 255;
g_imageData[ofs + 4*g_width + 1] = color;
g_imageData[ofs + 4*g_width + 3] = 255;
}

// Pull the string at the middle
function pull() {
// If there's no simulation running, set it up.
if (g_timer === -1) {
if (!setupCanvasAndStartSimulation())
return;
}
// Now, pull the middle string sample to minimum height.
// Also, mark its old value in the backlog, so we'll clear
// that old pixel in the main loop.
var yOldest = g_1dwave_y[(g_1dwave_idx+1) % 2];
g_pullBacklog.push(yOldest[Math.floor(g_width/2)]);
// And now, pull!
yOldest[Math.floor(g_width/2)] = g_height/2;
}

// The theory behind the 1d wave calculation is the following:
// Assume we have a wave function, f(x), (x=horizontal pixel) as it is at time (n).
// We want to find how the values will change at time (n+1).
//
// Consider the individual pixels to represent the water molecules.
// Assume that each molecule is influenced only by its two neighbours,
// as if they are connected to it with "springs":
//
//      xxx                   xxx
//     x    x                R
//    x      x              O._
//            x            L |\
//  ---------------------------\------
//              x        x      \
//               x      x        ---- This point (O) is only influenced
//                x    x               by the forces of two "springs":
//                 xxxx                the one attached to its left
//                                     neighbour(L) and the one attached
//                                     to its right neighbour(R)
// Using Euler integration of Neutonian laws, we will have for the point O:
// (p=position, v=velocity, a=acceleration)
//
//
// p    = p  + v        (new position = old position + velocity)
//  n+1    n    n
//
//
// v  = v    + a        (velocity = old velocity + old acceleration)
//  n    n-1    n-1
//
//                      (springs cause acceleration linear to distance)
//
// a    = k * (pL    - p    ) + k * (pR    - p    )  (k is the spring coeff)
//  n-1          n-2    n-2            n-2    n-2
//
//
// This would work, but would require quite a lot of calculations...
// Instead, we will use the acceleration at time n+2 (instead of n-1),
// which basically means that we will delay the effects of the springs
// for three time slots - as if the springs have a "delayed" effect.
//
// We will also set k=1. Look at what happens (we replace v  first):
//                                                         n
//
//
// p    = p  + v    + a    = (aprox) p  + (p  - p   )  + a
//  n+1    n    n-1    n-1            n     n    n-1      n+2
//
//      = 2*p  - p    + pL    - p   + pR    - p
//           n    n-1     n      n      n      n
//
//      = pL  + pR  - p
//          n     n    n-1
//
// So the end result is amazingly simple: we add the values of the two
// neighbouring pixels, and subtract the old value we used to have!
//
// Since the springs are losing power with each iteration, we will also
// scale the output of the last equation so that it diminishes over time
// (e.g.  multiply by a factor of 0.995):

function simulate1d() {
var image = g_dc.getImageData(0, 0, g_width, g_height);
g_imageData = image.data;
// Double buffering two string 'instances'
var yOldest = g_1dwave_y[g_1dwave_idx];
g_1dwave_idx = (g_1dwave_idx+1) % 2;
var yOld=g_1dwave_y[g_1dwave_idx], i=0;

// Erasing the pixels of the 'pull' backlog
for(; i<g_pullBacklog.length; i++) {
qplot(g_width/2, g_height/2 + g_pullBacklog[i], 0);
}
g_pullBacklog = [];

for(i=1;i<g_width-1; i++) {
// Erasing the old string instance pixels
qplot(i, g_height/2 + yOld[i], 0);
// Computing our awesome equation
var val = yOld[i-1] + yOld[i+1] - yOldest[i];
yOldest[i] = 0.995*val;
// Plotting the new string instance pixels
qplot(i, g_height/2 + yOldest[i], 255);
}
// Display the freshly computed HTML5 canvas image
g_dc.putImageData(image, 0, 0);
}

////////////////////////////////////////////////////////////////////
// 2D wave simulation
////////////////////////////////////////////////////////////////////

// When the user clicks anywhere in our blue pool...
function rainDrop(evt) {
// ...first figure out where he clicked, in terms of our canvas' native coordinates:
evt = evt || window.event;
var x = g_width*(evt.pageX - g_canvas.offsetLeft - $('#tabs')[0].offsetLeft)/g_canvas.clientWidth; var y = g_height*(evt.pageY - g_canvas.offsetTop -$('#tabs')[0].offsetTop)/g_canvas.clientHeight;

// And then, pull that water molecule to maximum height...
var yOldest = g_2dwave_y[(g_2dwave_idx+1) % 2];
yOldest[g_2d_offsets[Math.floor(y)]+Math.floor(x)] = 8.0;

// In fact, the maximum height is 1.0 - but by pulling it to 8,
// we will in fact be displaying this wave as 8 closely
// diminishing ones (see the scaling by 127 below and what
// it means for values larger than 255...)
}

// (Continuing from the comments of 'simulate1d' above...)
//
// ...for the corresponding two-dimensional problem,
// we just average the effects of the X- and Y- coordinate waves:
//
// p   (i,j) = 0.5 (p                  (i,j) + p                (i,j)
//  n+1              n+1,HorizontalWave         n+1,VerticalWave
//
//           = 0.5 (0.999(p (i-1,j) + p (i+1,j) - p   (i,j)) +
//                         n           n           n-1
//
//                + 0.999(p (i,j-1) + p (i,j+1) - p   (i,j)))
//                         n           n           n-1
//
// p   (i,j) = 0.999(0.5(p (i-1,j) + p (i+1,j) + p (i,j-1) + p (i,j+1))  - p   (i,j) )
//  n+1                   n           n           n           n             n-1
//
// Which translates to: we add the 4 neighbouring pixels, divide by 2,
// subtract the previous frame value, and scale by our energy loss coefficient.
//
// Simple! :-)

function simulate2d() {
var image = g_dc.getImageData(0, 0, g_width, g_height);
// Double buffering of the two pool 'instances'
var yOldest = g_2dwave_y[g_2dwave_idx];
g_2dwave_idx = (g_2dwave_idx+1) % 2;
var yOld=g_2dwave_y[g_2dwave_idx];
// Now apply the formula...
for(var i=1; i<g_height-1; i++) {
// ...and pre-compute the line offsets, for speed.
// (probably unnecessary in the case of V8,
//  but I did do this in the libSDL code, so I ported it over)
var     lineOffset = g_2d_offsets[i];
var prevLineOffset = g_2d_offsets[i-1];
var nextLineOffset = g_2d_offsets[i+1];
for(var j=1; j<g_width-1; j++) {
var val = (
yOld[prevLineOffset + j-1] + yOld[prevLineOffset + j+1] +
yOld[nextLineOffset + j-1] + yOld[nextLineOffset + j+1])/2 - yOldest[lineOffset + j];
yOldest[lineOffset + j] = (0.99999*val);
// Scale between 0 and 255...
var color = (127.0+127.0*0.99999*val);
// R G [B] [A] - water is blue, no?
image.data[4*(lineOffset+j) + 2] = color;
image.data[4*(lineOffset+j) + 3] = 255;
}
}
g_dc.putImageData(image, 0, 0);
}

// The function called by setInterval every 20ms (maximum refresh rate: 50Hz).
function simulate() {
// Depending on the active TAB, call the proper simulator.
if (g_activeTab === 0) {
simulate1d();
} else if (g_activeTab === 1) {
simulate2d();
} else {
// In the rest of the TABs (the 'Theory' and 'Javascript' sections)
// disable the simulation (don't waste tablet/phone batteries!)
if (g_timer !== -1) {
clearInterval(g_timer);
g_timer = -1;
}
}
}

// Called from some HTML side links to switch the active jQueryUI tab.
function switchToTab(tab) {
jQuery('#tabs').tabs("option", "active", tab);
return false;
}

// The main entrypoint:
// Create a jQueryUI tab, and show the first tab ('1D wave'):
jQuery( "#tabs" ).tabs({
activate: function(event, ui) {
g_activeTab = {
'1-dimension wave (string)': 0,
'2-dimension waves (water)': 1,
'The theory':                2,
'The code':                  3
}[ui.newTab[0].children[0].text];
if (g_activeTab === 0 || g_activeTab === 1)
setupCanvasAndStartSimulation();
else if (g_activeTab === 2) {
if (g_mathJax_configured === false) {
g_mathJax_configured = true;
MathJax.Hub.Configured();
}
}
},
active: 1
});
setupCanvasAndStartSimulation();
});
`