Issue #7 - December 2000

by Zip <zippy_@hotmail.com>

In this tutorial, I will explain the math behind making a simple ripple effect, and hopefully provide a working demonstration. I can't think of any practical application of such an effect beyond just a nice graphics demo, but it's quite nice, nonetheless. We will start with the one dimensional case, then build upon that for a 2D implementation. It is assumed that your are no less than an intermediate programmer and that you have a solid background in math (specifically trigonometry and algebra), though the topics discussed here should be easily understood. All I know of this effect is what I found on my own, so it's not guaranteed to be physically accurate in any way. But it looks neato. ;)

The ripple follows a sine wave, and the amplitude of the wave decreases as the distance from the origin increases. The viewer looks directly down at the wave, and the light ray is refracted by .5 pi radians from the angle of the instantaneous slope of that point. So we just take the inverse of the instantaneous slope, and project that to the X axis. Perhaps a few diagrams will help.

y^ pi/2 y = sin x | | | . . | . . | . . | . . | . . | . . pi radians 2 pi radians |. . / \ ..---------------------.----------------------.-> x |\ . . | 0 . . | . . | . . | . 3pi/2 . | . | . | . . | v

As you can see, the slope of the curve is steepest at x = 0 and x = pi, and least steep at x = pi/2 and x = 3pi/2. Also, notice that the instantaneous slopes of the curve range from 1 to -1, as do the values of the sine function, and that the instantaneous slope when x = 0 is 1, 0 at x = pi/2, etc.

And here's the cosine function:

y^ y = cos x | |. . |\ . . \ | 0 . . 2pi | . . | . . | . pi/2 3pi/2 . | . / \ . +----------.----------------------.-----------------> x | . . | . . | . . | . . | . pi . | . | . | . . | v

As you can probably see, the cosine function is the instantaneous slope (or first derivative, I believe) of the sine function. Now let's look at how light is refracted in a wave. The method I will describe here makes a big assumption, and neglects the fact that different liquids may refract light differently than others. However, this is generally good enough for our purposes. Let's say that a light ray going straight down to the wave is refracted at 90 degrees, or pi/2 radians, to the instantaneous slope at that point. The ray is then projected to the X axis. For example:

y^ <O> / | | / -instantaneous slope | | / . . | | /. . | | / . | |. inverse of . | .\ / instantan- . | / \ eous slope . |/ \ . ..-----\---------------.----------------------.-> x | . . | . . | . . | . . | . . | . . | . . | v

Now let's move the sine wave up a bit, since the X axis will be our background picture, and wave usually don't go below the solid object under the liquid. The distance we move the graph of the sine wave upward will be the depth of the liquid, which I will call d.

y^ y = d + sin x | | . . | . . | . . | . . | . . | . . |. . .. . . -+ | . . | | . . | | . . | | . . | | . . +- d | . . | | . . | | | | | |-----------------------------------------------> x -+ v

Now we need a formula to project the vector to the X axis. We shall derive this formula from the equation for a line, y=mx+b. First, let's consider what's given, what we've already figured out, and what remains to be solved. x1 is given, and y1 is the sine of X plus the depth. y2 is always 0, since it's being projected to the X axis. m1 is the inverse of the instantaneous slope of the line, which is one divided by the cosine of X, and m2 is the same value. Now you may be wondering how we would calculate b, the Y intercept. All we need to do with that is express it in terms of x, y, and m. The Y intercept can be expressed as y1-mx. Now I shall solve the equations for x2.

First, a few variable definitions:

x = x1 m = m1 = m2 = 1 / cos(x) y = y1 = d + sin(x) y2 = 0 u = x2 b = b1 = b2 = y - mx

Now the derivation:

y1 = m1*x1 + b1 y = mx + b y - b = mx (y - b) / x = m y2 = m2*x2 + b2 0 = mu + b -b = mu -b/u = m ..·. (y - b) / x = -b/u (u(y - b)) / x = -b u(y - b) = -bx u = -bx / (y - b) u = ((-y-mx)x) / (y - (y-mx)) u = -((y-(x/cos(x)))x) / (y - (y-(x/cos(x)))) u = -((y-(x/cos(x)))x) / (x/cos(x)) u = -(cos(x)(y-(x/cos(x)))x) / x u = -cos(x)(y - (x / cos(x))) u = -cos(x)y - (x*cos(x))/cos(x) u = -(cos(x)y - x) u = x - cos(x)y u = x - cos(x)(d + sin(x))

I hope that wasn't too hard to follow. You may be asking what good this is, as it obviously can't be used for a ripple effect with a background. Therefore, we shall use what we derived here to develop routines for the 2D case, which will appear 1D as it will be displayed here, but could undoubtedly be utilized in a 3D engine without too much trouble. If it still doesn't make sense, however, then run this little program. It is a graphic illustration of what we're doing here.

SCREEN 12 PI = 3.14159 DEG = PI / 180 DIM BG(639) AS INTEGER FOR A = 0 TO 639 BG(A) = 1 + RND * 14 PSET (A, 440), BG(A) NEXT A FOR X = 0 TO PI * 200 XX = X / 100 YY = SIN(XX) Y = 240 - YY * 200 XX2 = XX + (1 + SIN(XX)) * COS(XX) X2 = XX2 * 100 LINE (X, Y)-(X2, 440), 8, , &HAAAA LINE (X, 0)-(X, Y), 8, , &HAAAA IF X2 >= 0 AND X2 <= 639 THEN PSET (X, Y), BG(X2) WAIT &H3DA, 8 WAIT &H3DA, 8, 8 LINE (X, Y)-(X2, 440), 0 LINE (X, 0)-(X, Y), 0 IF X2 >= 0 AND X2 <= 639 THEN PSET (X, Y), BG(X2) PSET (X2, 440), BG(X2) END IF NEXT X

What?! All that work for one little formula?? That's right. Also, you may have noticed that the above program uses the formula u=x+(d+sin(x))cos(x) instead of u=x-(d+sin(x))cos(x). This is because the positive direction is down instead of up on the computer screen. Now maybe it makes a little more sense, but if you're still confused, then just read all that again until you understand it.

Although I'm calling this the 2D case, it's actually 3D, because the wave moves up and down in the 3rd dimension. For simplicity, I'm going to call the horizontal axis the X axis, the vertical the Y axis, and the axis going directly into the screen the Z axis, though this is not how they are normally named in math classes. To try to help you visualize how the wave would appear in 3D space, look at the picture below.

y ^ | | | | | | ..... ..ooooo.. .ooOOOOOoo. .oOOoooooOOo. .oOoo.....ooOo. .oOo.. | ..oOo. <------------------.oOo.---+---.oOo.------------------> x .oOo.. | ..oOo. .oOoo.....ooOo. .oOOoooooOOo. .ooOOOOOoo. ..ooooo.. ..... | | | | | | v

Notice that the ripple will be along a single plane, but won't be right on it. That's why I'm calling it 2D and 3D at the same time. In the 1D case, the height of the wave was the sine of the distance from the origin, which was X. However, in 2D, the distance of a point from the origin is the square root of the sum of the squares of the x and y values because of the Pythagorean theorem (z = sqr(x²+y²)). The next part may get a tiny bit confusing, but not terribly so. You may remember that the angle of a vector is equal to the arctangent of its slope, so we have theta=arctan(y/x). What we want is a 2D vector that will point to the correct point on our background picture after being refracted through the sine wave. I'll call the x & y components of this vector u & v, respectively. From the 1D case, we know that:

s = d - y/m

where s is displacement, d is distance along the X or Y axis, m is the inverse of the instantaneous slope of that point, and y is the height of the wave at that point. Now we shall solve this equation for the X and Y components. But before we do that, let's think about what notation to use. We have a vector to the point before refraction in polar notation (angle & magnitude), but we need rectangular coordinates for the final vector. Therefore, we can express the vector from the origin to the point on the wave as:

V = (sqr(x²+y²)cos(arctan(y/x)))i + (sqr(x²+y²)sin(arctan(y/x)))j

So what did that accomplish? Good question. Er, oh yeah. Think about it. When you refract a light ray, its angle in the plane of the ripple does not change; only its vertical angle will change. Therefore, all we have to do is calculate how much this vertical angle changes and project it to the x-y plane. Now, solving the equation above for s:

s = d - y/m s = sqr(x²+y²) - (d + sin(sqr(x²+y²)))cos(sqr(x²+y²)) ^ (depth)

What is this exactly, you ask? It is the magnitude of the vector after refraction. But what's the angle? It's the same as it was before, because the angle in the x-y plane does not change. Using this new magnitude, we get the vector:

V = ((sqr(x²+y²) - (d + sin(x²+y²))cos(x²+y²))cos(arctan(y/x)))i + ((sqr(x²+y²) - (d + sin(x²+y²))cos(x²+y²))sin(arctan(y/x)))j

Ok, I know that looks really nasty and complicated at first, but it's really not that bad. All I did is take apply the concepts from the 1D case to the 2D case.

This part is quite simple. All you do is scan an area of pixels around the origin and run the coordinates through those formulas. Keep in mind that the values can be negative, so you will want to add some value (i.e. half the size of the screen) to the values calculated. Pseudo-code would be:

for y = -100 to 100 for x = -100 to 100 u = (sqr(x²+y²) - (d + sin(sqr(x²+y²)))*cos(sqr(x²+y²)))*cos(atn(y/x)) v = (sqr(x²+y²) - (d + sin(sqr(x²+y²)))*cos(sqr(x²+y²)))*sin(atn(y/x)) if (point(u+100,v+100) is on the screen) { plot x+100, y+100, BGpic(u+100,v+100) } next x next y

Ah, yes. I can hear you screaming now.. "You said you were gonna teach how to make a ripple, but this doesn't even move!". Ok. All we need to do for that is have a variable (w) that is decremented each frame and added to the magnitude used to calculate the sine & cosine in the magnitude of the vector:

u = (sqr(x²+y²) - (d + sin(sqr(x²+y²)+w))*cos(sqr(x²+y²)+w))*cos(atn(y/x)) v = (sqr(x²+y²) - (d + sin(sqr(x²+y²)+w))*cos(sqr(x²+y²)+w))*sin(atn(y/x))

Also, remember that the amplitude of the wave decreases as the magnitude increases. I will say that variable n is equal to the magnitude of the unrefracted ray (sqr(x²+y²)), or one if the magnitude is zero. You may also wish to have a maximum amplitude (a) other than one. Changing the above formulas to reflect these changes give us:

u = (sqr(x²+y²) - (d + a * (sin(sqr(x²+y²)+w) / n)) * a * cos(sqr(x²+y²)+w)) * cos(atn(y/x)) v = (sqr(x²+y²) - (d + a * (sin(sqr(x²+y²)+w) / n)) * a * cos(sqr(x²+y²)+w)) * sin(atn(y/x))

That's right, but there are a few things that can be optimized. Notice that the only difference between the formulas for u & v are that u is the magnitude multiplied by the cosine of the angle while v is the magnitude multiplied by the sine. Therefore, you only have to calculate the magnitude once. You can also use look-up tables for the sine and cosine. The magnitude of the unrefracted vector (sqr(x²+y²)) needs to be calculated only once, as well. Finally, you will definitely want to use an assembly language library to gain as much speed as possible. Double-buffering will help you to gain some speed, too. Even with all this, it will not run too quickly unless you have a fast computer or find more optimizations. It may help if you make the ripple be only a small part of the screen, since redrawing the entire screen one pixel at a time can be very slow.

I changed my mind. I'm only going to give you pseudo-code, as that increases the chances
of you understanding the concepts presented before making a working example. If you want
to see my ripple demo, go to http://www.angelfire.com/co/zippy15/wtrz.zip. [*Editor's Note:
The ripple demo is included in the downloadable version of this issue, in the file named
ripple.zip*]

set graphics mode dim buffer1(sizeOfScreen), buffer2(sizeOfScreen) loadPicTo "background.bmp", addressOf(buffer) w = 0 do { w = w - 1 degree (pi/180 radians) if w < 0 then w = 360 degrees (2pi radians) for y = top to bottom { deltaY = middleY - y for x = left to right { deltaX - middleX - x mag = sqr(deltaX² + deltaY²) n = mag if (!n) n = 1 i = mag + w m = mag - (depth + maxAmplitude * (sin(i) / n)) * maxAmplitude * cos(i) if (deltaX != 0) { theta = atn(deltaY/deltaX) } else { theta = atn(deltaY) } if (deltaX < 0) theta = theta - pi radians (180 degrees) u = m * cos(theta) v = m * sin(theta) if (u,v) is in bgPic boundaries { plotPixel x, y, colorAt(u,v) } } } copy buffer to screen } loop while (!keypressed)

Do something like that, and it should work. I hope you can make sense of my weird pseudo-code there. Anyway, that's pretty much all there is to it. Notice that this method does not consider reflection, but refraction only. And the angle of refraction may not be (and probably isn't) physically accurate. But, as you shall see, it's good enough.

Thanks for taking the time to read this tutorial. I hope you learned something new from it. If you have questions, comments, suggestions (i.e. something that should be better explained), etc., then email me.