# Wave Physics

## Recommended Posts

Hi there. I have a small problem regarding wave physics. I've taken this algorithm :
void CWave::update()
{
float n;
float temp;
float one_over_dampf = 1.0f / DampingFactor;

int h;
register int v;

//Skip the edges to allow area sampling
for(h = 1; h < mHorizontalSize; ++h)
{
for(v = 1; v < mVerticalSize; ++v)
{
n = (	mWave2[v - 1][h    ] +
mWave2[v + 1][h    ] +
mWave2[v    ][h - 1] +
mWave2[v    ][h + 1]) * 0.5f - mWave[v][h];

n -= n * one_over_dampf;

mWave[v][h] = n;
}
}

//wave=wave2

for(h = 0; h <= mHorizontalSize; ++h)
{
for(v = 0; v <= mVerticalSize; ++v)
{
temp = mWave2[v][h];
mWave2[v][h] = mWave[v][h];
mWave[v][h] = temp;
}
}
}


which is basically an update algo for a water wave. Now I've been told that water waves and explosion waves are relatively the same. The thing is, does anyone know how to prevent the waves from bouncing off the edges? I found this algo from a site that I can't remember and that guy who used it got it from the Programming Gems book. So, anyone have any ideas? Thanks a lot!

##### Share on other sites
To stop it bouncing set the values at the edges equal to the values one index in from the edges after each update (and only bother updating the internal points). Look up "boundary conditions".

Incidently - you can avoid that slow memory copy by using pointers (not 2D arrays) and just swapping pointers between updates.

##### Share on other sites
Thanks a lot!

So you will only update, lets say if h is the horizontal width and v the vertical height, x=1 to h-1 and y=1 to v-1 (if x and y are variables of a for statement) ?

And what about the corners?

Sorry about all the questions, I'm really desperate. =|

Thanks again.

##### Share on other sites
Yes, that's right. Corners don't matter since they never actually get referenced - you can set them to anything you like.

##### Share on other sites
Hi again.
I've done the changes and the waves don't bounce but I get a weird disturbance in the middle of the process.

Here's what I did:
void CWave::update(){    float n;    //float temp;    float one_over_dampf = 1.0f / DampingFactor;        int h;    register int v;        //Skip the edges to allow area sampling    for(h = 1; h < mHorizontalSize; ++h)	{        for(v = 1; v < mVerticalSize; ++v)		{            n = (	mWave2[v - 1][h    ] + 					mWave2[v + 1][h    ] + 					mWave2[v    ][h - 1] + 					mWave2[v    ][h + 1]) * 0.5f - mWave[v][h];            			n -= n * one_over_dampf;                        mWave[v][h] = n;						if(h == 1)			{				mWave2[v][0] = n;			}			else if(h == mHorizontalSize - 1)			{				mWave2[v][mHorizontalSize] = n;			}						if(v == 1)			{				mWave2[0][h] = n;			}			else if(v == mVerticalSize - 1)			{				mWave2[mVerticalSize][h] = n;			}        }    }        //swap wave1 and wave2        float **temp = mWave;    mWave = mWave2;    mWave2 = temp;}

Sorry to be a bother like this.

But thanks!

##### Share on other sites
Your loop goes up to mHorizontalSize and mVerticalSize yet you reference v+1 and h+1 - out of the array bounds. This might be causing your disturbance.

The other reason for your disturbance is that the simulation may be unstable. The (1 dimensional) wave equation is:

d2h/dt2 = c^2 d2h/dx2

(d2h/dt2 is the second derivative of h with respect to time, c is the wave speed).

Then this equation is discretised - second derivative terms can be written like this:

d2h(i)/dx2 = (h(i-1) + h(i+1) - 2*h(i))/dx^2

The d2h/dt2 term is done in the same way - a nieve implementation would need you to store 3 arrays of h (one for time t+dt, one for time t, one for time t-dt), but your implementation avoids this by overwriting one of the arrays as it goes.

Anyway, the thing is that your implementation is missing various terms like c, dt and dx - in effect you've set dx and dt to be 1 (with maybe a constant term I can't be bothered to work out). This is a bad thing to do, because the implementation is only stable if:

c^t dt^2 / dx^2 < 0.5

(this is known as the Courant-Friedrichs-Lewy or CFL criterion for stability), and since you've dropped those terms from your code, it's not immediately clear that this criterion is being met.

I think that you need read up a little bit on how that bit of code is derived from the maths - it's really not hard. I just checked and the GPG book does have it, but if you don't have that book there's stuff on the web - e.g.

http://www.osc.edu/education/su_programs/si/si1997/projects/wave.html

(I googled for wave equation discrete)

When you've rewritten that bit of code with the dx, dt and c terms, then you should choose the dt value such that the CFL criterion is satisfied - and in general is more than satisfied - for example choose dt such that

c^t dt^2 / h^2 = 0.05

since smaller dt values will also lead to a more accurate result.

I know it's tempting just to filch working physics code and plug it in... but unless you understand the basics behind it it tends to do wierd things!

Edit: Changed the link to one that remembers to include the dx/dt terms!!

[Edited by - MrRowl on September 5, 2005 4:14:03 PM]

##### Share on other sites
Thanks a lot MrRowl.

The thing is in our project I wasn't originally in charge with the physics but another guy in our group couldn't continue anymore so I had to take over and I honestly havn't got a clue why he did it like that thats why I tried playing around with it and posting here to find some answers.

And the mHorizontal size thing is correct because in the landscape my friend who is no longer with our group made the width and height of the arrays to be [mVerticalSize + 1][mHorizontalSize + 1] so if you check if something <= mVerticalSize then you're checking if something is < mVerticalSize + 1. Get it? I really want to change that because its pretty confusing sometimes and unnessasary.

Sorry I don't mean to talk bad of my friend who left, he's a good hardworking guy who didn't leave because he wanted to but because cercomstances forced him to.

Thanks again.

##### Share on other sites
Ok, I've implemented the velocities, times, etc. (got it from a tutorial) so I have:
float A = (3.0f * 0.10f / 1.0f) * (3.0f * 0.10f / 1.0f);float B = 2 - 4 * A;int   h, v;float one_over_dampf = 0.99;//Boundaries should not be updatedfor(h = 1; h < mHorizontalSize; ++h){	for(v = 1; v < mVerticalSize; ++v)	{		mWave[v][h] =	A * (	mWave2[v    ][h - 1]  + 					mWave2[v    ][h + 1]  + 					mWave2[v - 1][h    ]  + 					mWave2[v + 1][h    ]) + 				B * 	mWave2[v    ][h    ] - 					mWave [v    ][h    ];				mWave[v][h] *= one_over_dampf;				if(h == 1)		{			mWave[v][0] = mWave[v][1];		}		else if(h == mHorizontalSize - 1)		{			mWave[v][mHorizontalSize] = mWave[v][mHorizontalSize - 1];		}				if(v == 1)		{			mWave[0][h] = mWave[1][h];		}		else if(v == mVerticalSize - 1)		{			mWave[mVerticalSize][h] = mWave[mVerticalSize - 1][h];		}	}}       //swap wave1 and wave2      float **temp = mWave;   mWave = mWave2;   mWave2 = temp;

It still bounces. Sorry man, I tried a whole lot of diffrent things to stop the waves in their tracks, I've searched on the net and I even read that page you posted but no luck. I can understand if you are not going to post again. But thanks a lot so far for your help MrRowl.

##### Share on other sites
Ok, I've changed my if statement to this:

if(h == 1){	mWave2[v][0] = mWave2[v][1];}else if(h == mHorizontalSize - 1){	mWave2[v][mHorizontalSize] = mWave2[v][mHorizontalSize - 1];}if(v == 1){	mWave2[0][h] = mWave2[1][h];}else if(v == mVerticalSize - 1){	mWave2[mVerticalSize][h] = mWave2[mVerticalSize - 1][h];}

This prevented the bouncing and it still remained relatively stable but there was some bouncing back action but it was only 2 or 3 very small waves instead of 4 huge waves coming back. Any thoaghts?

##### Share on other sites
As you worked out - you need to make sure you apply the boundary condition to the correct array. It's probably simpler and "safer" to make a separate loop and apply it to mWave2 before your update since otherwise there's a danger that your boundary condition gets applied out of step at some of the edges (not got time to check if this is really the case).

Also - my suggestion for a boundary condition with mWave2[0][*] = mWave2[1][*] probably isn't good - sorry :). Some things to try:

---
mWave2[0][j] = 0; // i.e. simple fixed boundary
---
mWave2[0][j] = mWave2[1][j] + (mWave2[2][j] - mWave2[1][j]); // i.e. height gradient is continuous
---
// Gradually damp over a border zone - something like this
for (unsigned i = 0 ; i < 5 ; ++i)
mWave2[i][j] *= ((float) i) / 5;

and similarly for the other edges.

Sorry I don't have time to check if these are actually good ideas :) From a quick web search it seems like it's not completely trivial to totally eliminate reflected waves.

##### Share on other sites
I tried most of that stuff (and I hope I implemented them correctly) but they don't work.
I just want to know. What is it that makes a wave reflect? I saw that if you specify a range in the calculations where the damping stuff is 0 then the waves reflect but what should you make the damping to cancel the wave out?

##### Share on other sites
Can't you just make the damping really high near the borders?
The problem then is that any wave created inside these areas wouldn't spread a lot. But maybe for your purpose it's good enough.

There's also the idea of mirroring the waves. As anyone knows, this cancels a wave out, and I've been able to do it while experimenting with this code. Maybe it's your best shot, worth trying.

##### Share on other sites
How can you mirror a wave to cancel it out? I've googled for it but I found nothing relevant on the subject.

##### Share on other sites
I am having a siniliar issue with me sea game. the problem is, how do i create the waves in the first place? What direct x code will give me water? I am working with dev c++. Thanks.

##### Share on other sites
I'm also working in wxDevC++.
The thing is that DirectX doesn't have direct code for you to create water graphics, you'll need to write the stuff on your own, I saw this yesterday when I looked through the DirectX 8 SDK and saw the water example. I'm not sure about DirectX 9 though.

EDIT :
To create a wave you just set your wave's xy point (which you use to draw, ie. like my mWave var) to a big value, like mine is 80. Then the algorithm will sort itself out.

Anyway, I saw my lecturer yesterday and we've come to a cool discovery.

EDIT:
void CWave::update(){    float A = (3.0f * 0.2f / 1.0f) * (3.0f * 0.2f / 1.0f);    float B = 2 - 4 * A;    int   h, v;	    float one_over_dampf = 0.99f;        //Boundaries should not be updated    for(h = 1; h < mHorizontalSize; ++h)    {        for(v = 1; v < mVerticalSize; ++v)        {	    mWave[v][h] = 	A * (	mWave2[v    ][h - 1]  + 		        		mWave2[v    ][h + 1]  + 					mWave2[v - 1][h    ]  + 					mWave2[v + 1][h    ]) + 				B * 	mWave2[v    ][h    ] - 					mWave [v    ][h    ];	    //This is the filtering part	    if(mWave[v][h] < 0.00001f)		mWave[v][h] = 0.0f;	    else		mWave[v][h] *= one_over_dampf;	}    }        //swap wave1 and wave2        float **temp = mWave;    mWave = mWave2;    mWave2 = temp;}

This will prevent the bouncing but the waves runs out of energy real quick. So I tried playing around with the consept and found a few cool things but afterwards there is a weird thing going on in the middle, not that it goes wacked but its an irretating piece of wobbles (or small waves) that I really tried to filter out but that causes the waves to run out of energy even faster.

Any further thoaghts on the consept?

[Edited by - Last Attacker on September 9, 2005 9:48:50 AM]

##### Share on other sites
Sorry but I really can't see a noticeable difference between that piece of code and the one you had earlier :)

Anyways, back on my previous line of thought. In this sort of stuff I usually try to get as much inspiration from nature as possible; after all, nature got it right, and anything too different will feel artificial to us.
I think that a simplistic way to understand what happens when a wave hits a shore or a cliff is that it suddenly enters an small area (point of impact) wich not only reflects the wave, because it can't go on any more, but also has a high dampening value, because part of the wave's kinetic energy is released onto the rocks or whatever (that's how erosion works).

So, I really think you should have a border of 2 or 3 pixels that works as expected, except with a high dampening value. A small part of the waves will still reflect back, but that's how it happens naturaly anyway, and they're so small that they will die off quickly.

An improvement over this method would be accounting for the fact that obviously the border only reflects waves that are moving towards it, not perpendicular or away from it. It shouldn't be too hard to find the vector for the direction of the wave's movement for each pixel and, assuming the borders are lines along the X or Y axis, do something like this example for a left wall:
if (wave_velocity_x < 0) {
//Wave moving left, collision - high dampening
}

You could even make the dampening higher according to the wave's velocity on that axis, which would be even more realistic. If the border is irregular, and cannot be divided into simple line segments aligned with the X or Y axis, you'd have to play with the "surface"'s normal and the dot product. But maybe that's a bit over the top.

About mirroring the wave - I made it work in the middle of a "pond", not near a border that reflects waves. But it shouldn't be that hard. I'd imagine a construct like this for a left border:
 ___ ___ _________|   |   ||   |   || A | B |  C  ...|   |   ||___|___|_________

Areas B and C would be updating using the code you have. Area C goes on to the middle of the "pond". Area A would be special - it's constantly updated with the mirror image of area B. I imagine these 2 areas would be just 1 or 2 pixels thick. Given the behaviour of waves, my intuition tells me that if this doesn't work, areas A and B should be spaced by 1 pixel of the same as area C.

I'm sorry, that's a lot to digest! But I just can't help myself :) I hope it's helpful.

##### Share on other sites
Ok, so in my case the mirror image of edge at h=1 would be like this...?
mWave[v][1] = -mWave2[v][2];

Sorry, I'm not familiar with mirroring. I'm just guessing that the mirror image of a wave would be it's negative.

##### Share on other sites
What I mean by mirroring, in the example diagram above, would be mirroring over the X axis. As if there was a mirror between the 2 areas. I don't think it's possible to explain better than that :)

 ___ ___ _________|  v|v  ||ox | xo|| w | w |     ...|x  |  x||___|___|_________

## Create an account or sign in to comment

You need to be a member in order to leave a comment

## Create an account

Sign up for a new account in our community. It's easy!

Register a new account

• ## Partner Spotlight

• ### Forum Statistics

• Total Topics
627651
• Total Posts
2978406

• 10
• 12
• 22
• 13
• 33