How To: Make a 2D wave simulation Part 2

The discrete equation

What is it exactly that we're going to make? We're going to make a collection of snapshots of a surface over time, frames if you will. The surface will be represented as integer, yes, integer values inside an array. A 1D array. I'll get to that later.

This project and the code written here is for the C programming language. You can use any language you want and translate each line of code to any other language if you wish!

Recall the equation from the previous page.
Each point of the surface is going to have one integer index in the array. Therefore, the difference between, say, index i and the next index i+1 is 1. The difference with the previous element is also 1. That's dx and dy in the equation! dx2 and dy2 as well! They're aaaalll 1.
In the time domain, since we're talking about snapshots or frames, the difference between the current frame with the next is 1 frame, so dt and dt2 subsequently is 1.

Let's take the first derivative of the equation, with respect to time.
The first derivative, df/dt, is f(t+dt)-f(t). The second derivative therefore is (f(t+2dt)-f(t+dt))-(f(t+dt)-f(t)), simplified it gives f(t+2dt) - 2*f(t+dt) + f(t).
A few things to note here. f(t) represents a frame, a whole array. f(t+dt) represents the next frame and f(t+2dt) represents the frame after that. Therefore, for the simulation we're going to need the previous frame and the frame before the previous one. f(t+2dt) represents the frame to be computed!

From now on, the frame, or array of values for the height of each point of the surface, to be computed will be referred to as next[x+y]. The most recent one is curr[x+y] for current, and the one before that is prev[x+y].
Solving for next[x+y] we have:

next[x+y] = 2*curr[x+y] - prev[x+y] + (rest of equation);

The part before the "rest of equation" will be more useful if we define it as the time difference separately.
int dt = 2*curr[x+y] - prev[x+y];

The space part of the equation

Simlarly, the second derivative for each of the space variables x and y is (centered):

//for x:
int dx = 2*curr[x+y] - curr[(x+1)+y] - curr[(x-1)+y];

//for y:
int dy = 2*curr[x+y] - curr[x+(y+ROW)] - curr[x+(y-ROW)];

//Together:
//int ds = dx + dy;
//Or equivalently
int ds = 4*curr[x+y] - curr[(x+1)+y] - curr[(x-1)+y] - curr[x+(y+ROW)] - curr[x+(y-ROW)];

It's time to get a bit technical here and explain why use a 1D array. Simply because it's a bit more efficient!.
A 2D array, say a[ROW][COL], of row width ROW and column height COL, is laid out linearly in memory.
Think of a 2D matrix. To access the element at a[1][2] starting from the upper left element and going to the right, row by row, you'll have to get past 1 row worth of elements, plus 2, to get to the second element of the second row. Therefore, a[1][2] is nothing more than a[1*ROW+2] should you map the array linearly!!!
Therefore, the ROW in the above code is equal to just the width of the 2D surface! Let's define it as a constant:

#define ROW 128

...for a frame of width 128 elements, or a surface of 128 discrete points! You can use any size you want!
Of course, you can define the vertical width of the frame/surface, along the Y axis, with a constant like so:
#define COL 128

However, you'll notice later that we don't really need that when transversing the whole array, but we do need the size ROW*COL. Let's define it as:
#define SIZE ROW*COL

Now, we've almost written the full function that computes each frame!
With two for loops to transverse each array, we have:

#define ROW 128
#define COL 128
#define SIZE ROW*COL

const int c2 = 0x110; // will be explained later
void compute_next(int *next, int *curr, int *prev) {
  for (int x=0; x<SIZE; x+=ROW) {
    for (int y=0; y<ROW; y++) {
      int dt = 2*curr[x+y] - prev[x+y];
      int ds = 4*curr[x+y] - curr[(x+1)+y] - curr[(x-1)+y] - curr[x+(y+ROW)] - curr[x+(y-ROW)];
      next[x+y] = dt - c2*ds;
    }
  }
}

This is the core of the simulation!!!
Now, there's nothing left to do except define three buffers, one for the next, the curr, and prev arrays, and we're set! Right?

Not so fast!!!

You might have noticed already that the current code goes out of bounds for the value of x and y that describe the bounds of the array. We need what we call in the math world boundary conditions!.

For the sake of simplicity, we're gonna define the boundary conditions for the next array to be 0, but remember, this will invert the reflected wave! This is fine, unless you want to simulate a liquid's surface, in which case, we'll have to resort to another boundary condition where the velocity along the edges is 0.

The full code so far:

#include <stdio.h>

#define ROW 128
#define COL 128
#define SIZE ROW*COL

const int c2 = 0x110; // will be explained later
void compute_next(int *next, int *curr, int *prev) {
  for (int x=1; x<SIZE-ROW; x+=ROW) {
    for (int y=1; y<ROW-1; y++) {
      int dt = 2*curr[x+y] - prev[x+y];
      int ds = 4*curr[x+y] - curr[(x+1)+y] - curr[(x-1)+y] - curr[x+(y+ROW)] - curr[x+(y-ROW)];
      next[x+y] = dt - c2*ds;
    }
  }
  //edge along y=0
  for (int x=0; x<SIZE; x+=ROW) next[x] = 0;
  //edge along y=ROW-1
  for (int x=ROW-1; x<SIZE; x+=ROW) next[x] = 0;
  //edge along x=0
  for (int y=0; y<ROW; y++) next[y] = 0;
  //edge along x=SIZE-ROW
  for (int y=SIZE-ROW; y<SIZE; y++) next[y] = 0;
}

In the next pages, I'm going to show you how to use fixed point arithmetic and how to generate actual images using the data from the generated array!

Previous page Back Next page