13,705,955 members
alternative version

#### Stats

24K views
51 bookmarked
Posted 24 May 2015
Licenced CPOL

# Wave Simulator with C#

, 24 May 2015
Experimental wave simulation program made using VS 2010 C#

## Introduction

This is an experimental wave simulation program. It is used to produce waves in a particle pool. Water waves can be generated for example. Obstacles can be implemented. Reflection, diffraction and interaction of waves can be simulated.

## Background

I wanted to generate waves in computer environment using my own ideas. I will not go into deep physics since you can make a deep research on waves of course. I will only explain my own technique of wave generation. Here are the basic concepts.

### 1 - Square Pool

This is where the waves are generated. It is filled with particles. It's that greenish-blue square surface above.

### 2 - Particle

Particles resemble pixels. They are arranged in order and make up the whole pool. Each of them have height, velocity, acceleration, mass and sustainability (anti-damping). They move up or down only. Every particle is attracted to the other 8 neighbor particles in terms of height. In the illustration above, red one is the main particle and blue ones are neighbor particles. This attraction force is the only reason a particle moves naturally.

### 3 - Oscillator

It is used to generate sinusoid waves. It actively moves a certain particle up and down so neighbor particles are attracted to that particle and they move up and down too. This makes up a propagating wave. In the first image above, two oscillators are running at the bottom near the horizontal center.

### 4 - Obstacle

It is a static particle. It doesn't move. It is not affected by any force. It doesn't attract other particles. In the first image, yellow lines at the bottom are obstacles.

### How does it work?

The wave engine doesn't use any sophisticated algorithm. It simply uses the following.

`Force = Mass * Acceleration;`

All particles attract each other. They just want to stay close to each other in terms of height. Remember that they only change their altitude. They want to be at the same altitude. Attraction force is proportional to the height difference between particles.

When the engine makes calculations, it checks out every single particle. If a particle is an oscillator or static, the engine doesn't make force calculations on that particle and continues with the next. If a particle is normal dynamic one, calculations are made. It starts with the acceleration.

`Acceleration = (H - h) / m;`

`H` is the average height of the eight neighbor(surrounding) particles.

`h` is the height of the surrounded particle.

`m` is the mass.

Another alike formula is the following.

`Acceleration = Math.Sign(H - h) * Math.Pow(Math.Abs(H - h), P) / m;`

`P` is the power.

So now it gained an acceleration due to height difference. There is also a resistance.

This formula counteracts the acceleration. It somehow simulates damping and prevents further wave propagation.

`Acceleration -= Velocity / Sustainability;`

After that, the engine limits the acceleration to a prespecified value. So the absolute acceleration can't exceed the limit. `Limit` is necessary for eliminating both an uncertainty and overflow.

```if (Acceleration > Limit)
Acceleration = Limit;
else if (Acceleration < -Limit)
Acceleration = -Limit;```

Now the acceleration feeds velocity.

`Velocity += Acceleration;`

The velocity also feeds height but it doesn't fully feed it at once.

`Height += Velocity / ActionResolution;`

Limit the height values.

```if (Height > Limit)
Height = Limit;
else if (Height < -Limit)
Height = -Limit;```

There is also a value called `ActionResolution` to prevent positive feedback. This is important to understand because if particles make all movements at once, a chaos is unavoidable. So what is the problem?

1 - Imagine two particles which want to get to the same altitude but are in different altitudes.

`d` is the distance between particle and origin so these particles have a distance of `2d` between them. What happens in this case? They first gain a lot of acceleration. Acceleration is reduced to a degree due to damping. Acceleration feeds velocity. Now it's the time for feeding the height. If we move these particles at once, the next case will most likely be something like this.

They got furher instead of getting closer. They will never stop. They would stop if we wouldn't let them move at once. They have to go on their way with small steps. Every step must lead to a distance shorter than the previous one. This would prevent the situation above and they finally would stop at the same altitude.

For the last step, all of the particles shift to the origin point as a whole. Otherwise, they may wander away from the origin since there is not any force which pulls them to the origin point. To do this, all of the particles' height values are summed up and the resultant is divided by the number of particles and hence the average height of the particles is found. We subtract this value from the height of each particle.

```int Total_Height = 0;

for (int index = 0; index < Height_Array.Length; index++)
Total_Height += Height_Array[index];

for (int index = 0; index < Height_Array.Length; index++)
Height_Array[index] -= Total_Height / Height_Array.Length;```

This is the end of calculations. After all, the wave engine depicts the particles and renders it on a control.

## Using the code

There is a class called `WaveEngine.cs` to do all calculations and rendering. Implement it to your project first.

```using System;
using System.Drawing;
using System.Drawing.Imaging;
using System.Windows;
using System.Windows.Forms;

namespace WaveSimulator
{
public class WaveEngine
{
Bitmap bmp; // This will be our canvas.
BitmapData bd; // This will be used to modify the RGB pixel array of the "bmp".

// "vd" means "vertex data"
float[] vd; // Height map
float[] vdv; // Velocity map
float[] vda; // Acceleration map
float[] vds; // Sustainability map. Sustainability can be thought as anti-damping.
bool[] vd_static; // Static particle map. Particles which will act like a obstacle or wall.

float mass = 0.1f; // Mass of each particle. It is the same for all particles.
float limit = 500f; // Maximum absolute height a particle can reach.
float action_resolution = 20f; // Resolution of movement of particles.
float sustain = 1000f; // Anti-damping. Propagation range increases by increasing this variable. Minimum is 1f.
int delay = 1; // Time-out in milliseconds for force calculations.
float phase1 = 0f; // Current phase value of oscillator1.
float phase2 = 0f; // Current phase value of oscillator2.
float freq1 = 0.2f; // Phase changing rate of oscillator1 per calculation. Frequency increases by increasing this variable.
float freq2 = 0.2f; // Phase changing rate of oscillator2 per calculation. Frequency increases by increasing this variable.
float power = 1.0f; // Power of the output force exerted on each particle. Natural value is 1.0f

BufferedGraphics bufgraph; // Double-buffered graphics for rendering. It minimizes flickering.
BufferedGraphicsContext bufgcont; // Will be used to initialize bufgraph.

bool work_now = false; // True = Thread must make calculations now, False = Thread must sleep now.

bool highcont = false; // High contrast drawing.

bool disposing = false; // It will be true once the termination starts.

bool osc1active = false; // Is oscillator1 is turned on?
bool osc2active = false; // Is oscillator2 is turned on?

int osc1point = 0; // Location of the oscillator1 in the wave pool. It is an index value.
int osc2point = 0; // Location of the oscillator2 in the wave pool. It is an index value.

int size = 200; // Size of the wave pool. It indicates both the width and height since the pool will always be a square.

Color color1 = Color.Black; // Color of the crest or trough.
Color color2 = Color.Cyan; // Color of the crest or trough.

Color colorstatic = Color.Yellow; // Color of the static particles.

// These variables are used for edge absorbtion. It is used for eliminating reflection from window boundaries.
int absorb_offset = 10; // Offset from each window boundary where the sustainability starts to decrease.
float min_sustain = 2f; // The lowest sustainability value. They are located at the boundaries.
bool edge_absorbtion = true; // If true, the particles near the boundaries will have low sustainability.

Control control; // This will be the control where the engine runs and renders on.

public enum ParticleAttribute
{
Height = 1,
Velocity = 2,
Acceleration = 4,
Sustainability = 8,
Fixity = 16,
All = 32,
}

public float Mass
{
get { return mass; }
set
{
if (value > 0f)
{
mutex.WaitOne();
mass = value;
mutex.ReleaseMutex();
}
}
}

public float Limit
{
get { return limit; }
set
{
if (value > 0f)
{
mutex.WaitOne();
limit = value;
mutex.ReleaseMutex();
}
}
}

public float ActionResolution
{
get { return action_resolution; }
set
{
if (value >= 1f)
{
mutex.WaitOne();
action_resolution = value;
mutex.ReleaseMutex();
}
}
}

public float Sustainability
{
get { return sustain; }
set
{
if (value >= 1f)
{
mutex.WaitOne();
sustain = value;
setSustain();
mutex.ReleaseMutex();
}
}
}
public int Delay
{
get { return delay; }
set
{
if (value >= 0)
{
mutex.WaitOne();
delay = value;
mutex.ReleaseMutex();
}
}
}
public float PhaseRate1
{
get { return freq1; }
set
{
if (value > 0f && value < Math.PI * 2f)
{
mutex.WaitOne();
freq1 = value;
mutex.ReleaseMutex();
}
}
}
public float PhaseRate2
{
get { return freq2; }
set
{
if (value > 0f && value < Math.PI * 2f)
{
mutex.WaitOne();
freq2 = value;
mutex.ReleaseMutex();
}
}
}
public float Power
{
get { return power; }
set
{
if (power > 0f)
{
mutex.WaitOne();
power = value;
mutex.ReleaseMutex();
}
}
}
public int Size
{
get { return size; }
set
{
if (size >= 1f)
{
mutex.WaitOne();
size = value;
setPool();
mutex.ReleaseMutex();
}
}
}
public float EdgeSustainability
{
get { return min_sustain; }
set
{
if (value >= 1f)
{
mutex.WaitOne();
min_sustain = value;
setSustain();
mutex.ReleaseMutex();
}
}
}
public int AbsorbtionOffset
{
get { return absorb_offset; }
set
{
if (value > 0 && value < size / 2)
{
mutex.WaitOne();
absorb_offset = value;
setSustain();
mutex.ReleaseMutex();
}
}
}
public Color Color1
{
get { return color1; }
set
{
mutex.WaitOne();
color1 = value;
mutex.ReleaseMutex();
}
}
public Color Color2
{
get { return color2; }
set
{
mutex.WaitOne();
color2 = value;
mutex.ReleaseMutex();
}
}
public Color ColorStatic
{
get { return colorstatic; }
set
{
mutex.WaitOne();
colorstatic = value;
mutex.ReleaseMutex();
}
}
public bool HighContrast
{
get { return highcont; }
set
{
mutex.WaitOne();
highcont = value;
mutex.ReleaseMutex();
}
}
public bool EdgeAbsorbtion
{
get { return edge_absorbtion; }
set
{
mutex.WaitOne();
edge_absorbtion = value;
setSustain();
mutex.ReleaseMutex();
}
}

public bool Oscillator1Active
{
get { return osc1active; }
set
{
mutex.WaitOne();
osc1active = value;
setSustain();
mutex.ReleaseMutex();
}
}

public bool Oscillator2Active
{
get { return osc2active; }
set
{
mutex.WaitOne();
osc2active = value;
setSustain();
mutex.ReleaseMutex();
}
}

public Point Oscillator1Position
{
get { return new Point(osc1point % size, (int)Math.Floor((float)osc1point / (float)size)); }
set
{
if (value.X + value.Y * size < size * size)
{
mutex.WaitOne();
osc1point = value.X + value.Y * size;
setSustain();
mutex.ReleaseMutex();
}
}
}

public Point Oscillator2Position
{
get { return new Point(osc2point % size, (int)Math.Floor((float)osc2point / (float)size)); }
set
{
if (value.X + value.Y * size < size * size)
{
mutex.WaitOne();
osc2point = value.X + value.Y * size;
setSustain();
mutex.ReleaseMutex();
}
}
}

/// <summary>
/// Initializes the WaveEngine
/// </summary>
/// <param name="control">The control where the engine renders on.</param>
public WaveEngine(Control control)
{
this.control = control;
control.Resize += new EventHandler(control_Resize);
setPool();
mutex = new Mutex();
{

while (!disposing)
{
try
{
while (work_now)
{
mutex.WaitOne();
int beginning = System.Environment.TickCount;
while (System.Environment.TickCount - beginning < delay)
CalculateForces();
generatebitmap();
bufgraph.Graphics.DrawImage(bmp, 0, 0, control.ClientSize.Width, control.ClientSize.Height);
bufgraph.Render();
mutex.ReleaseMutex();
}
}
catch
{ work_now = false; mutex.ReleaseMutex(); }
}

});

ForceCalcT.Start();

}

void control_Resize(object sender, EventArgs e)
{
{
mutex.WaitOne();
if (bufgraph != null)
bufgraph.Dispose();

if (bufgcont != null)
bufgcont.Dispose();

bufgcont = new BufferedGraphicsContext();

bufgraph = bufgcont.Allocate(control.CreateGraphics(), control.ClientRectangle);
mutex.ReleaseMutex();
});

}

/// <summary>
/// Sets particles' specified attribute(s) to a specified value in a specified rectangular area.
/// </summary>
/// <param name="rect">Rectangular area which contains particles.</param>
/// <param name="value">Value to set the particles to.</param>
/// <param name="partatt">Attribute(s) that will be set.</param>
public void SetParticles(Rectangle rect, float value, ParticleAttribute partatt)
{
mutex.WaitOne();

if (rect.X < 0)
rect.X = 0;

if (rect.Y < 0)
rect.Y = 0;

if (rect.Width + rect.X > size)
rect.Width -= (rect.X + rect.Width) - size;

if (rect.Height + rect.Y > size)
rect.Height -= (rect.Y + rect.Height) - size;

bool xh = false, xv = false, xa = false, xs = false, xf = false;
// Let's see which attributes we are gonna deal with.
if ((ParticleAttribute.All & partatt) == ParticleAttribute.All)
{
xh = true; xv = true; xa = true; xs = true; xf = true;
}
else
{
if ((ParticleAttribute.Height & partatt) == ParticleAttribute.Height)
xh = true;
if ((ParticleAttribute.Velocity & partatt) == ParticleAttribute.Velocity)
xv = true;
if ((ParticleAttribute.Acceleration & partatt) == ParticleAttribute.Acceleration)
xa = true;
if ((ParticleAttribute.Sustainability & partatt) == ParticleAttribute.Sustainability)
xs = true;
if ((ParticleAttribute.Fixity & partatt) == ParticleAttribute.Fixity)
xf = true;
}

for (int y = rect.Y * size; y < rect.Y * size + rect.Height * size; y += size)
{
for (int x = rect.X; x < rect.X + rect.Width; x++)
{
if (xh)
vd[x + y] = value;
if (xv)
vdv[x + y] = value;
if (xa)
vda[x + y] = value;
if (xs)
vds[x + y] = value;
if (xf)
vd_static[x + y] = Convert.ToBoolean(value);
}
}
mutex.ReleaseMutex();
}

/// <summary>
/// Gives a float array of specified attribute of particles in a specified rectangular area.
/// </summary>
/// <param name="rect">Rectangular area which contains particles.</param>
/// <param name="partatt">Attribute whose array will be given. Only one attribute can be specified and "All" cannot be specified.</param>
public float[] GetParticles(Rectangle rect, ParticleAttribute partatt)
{
float[] result = new float[1];

bool xh = false, xv = false, xa = false, xs = false, xf = false;

if ((int)partatt == 1 || (int)partatt == 2 || (int)partatt == 4 || (int)partatt == 8 || (int)partatt == 16)
{
mutex.WaitOne();

if (rect.X < 0)
rect.X = 0;

if (rect.Y < 0)
rect.Y = 0;

if (rect.Width + rect.X > size)
rect.Width -= (rect.X + rect.Width) - size;

if (rect.Height + rect.Y > size)
rect.Height -= (rect.Y + rect.Height) - size;

result = new float[rect.Width * rect.Height];

if (partatt == ParticleAttribute.Height)
xh = true;
if (partatt == ParticleAttribute.Velocity)
xv = true;
if (partatt == ParticleAttribute.Acceleration)
xa = true;
if (partatt == ParticleAttribute.Sustainability)
xs = true;
if (partatt == ParticleAttribute.Fixity)
xf = true;

int r = 0;
for (int y = rect.Y * size; y < rect.Y * size + rect.Height * size; y += size)
{
for (int x = rect.X; x < rect.X + rect.Width; x++)
{
if (xh)
result[r] = vd[x + y];
if (xv)
result[r] = vdv[x + y];
if (xa)
result[r] = vda[x + y];
if (xs)
result[r] = vds[x + y];
if (xf)
result[r] = Convert.ToSingle(vd_static[x + y]);
r += 1;
}
}
mutex.ReleaseMutex();
}

return result;
}
/// <summary>
/// Starts the force calculation.
/// </summary>
public void Run()
{
work_now = true;
}

/// <summary>
/// Suspends the force calculation indefinitely.
/// </summary>
public void Stop()
{
work_now = false;
}

public void Dispose()
{
work_now = false;
disposing = true;
{
mutex.WaitOne();
bmp.Dispose();
mutex.Close();
});
}

void CalculateForces()
{

float total_height = 0;// This will be used to shift the height center of the whole particle system to the origin.

// This loop calculates the forces exerted on the particles.
for (int index = 0; index < vd.Length; index += 1)
{
// If this is a static particle, it will not move at all. Continue with the next particle.
if (vd_static[index])
{
vd[index] = 0;
vdv[index] = 0;
vda[index] = 0;
continue;
}

if (index == osc1point && osc1active)
{
// This is where the oscillator1 is located. It is currently active.
// So this particle only serves as an oscillator for neighbor particles.
// It will not be affected by any forces. It will just move up and down.
vdv[index] = 0;
vda[index] = 0;
vd[index] = limit * (float)Math.Sin(phase1);
phase1 += freq1;
if (phase1 >= 2f * (float)Math.PI)
phase1 -= (float)Math.PI * 2f;

continue;
}

if (index == osc2point && osc2active)
{
vdv[index] = 0;
vda[index] = 0;
vd[index] = limit * (float)Math.Sin(phase2);
phase2 += freq2;
if (phase2 >= 2f * (float)Math.PI)
phase2 -= (float)Math.PI * 2f;

continue;
}

// So this particle is neither an oscillator nor static. So let's calculate the force.

// Reset the acceleration. We do this because acceleration dynamically changes with the force.
vda[index] = 0;

// Sum up all the height values so we will find the average height of the system.
// This doesn't contribute to the force calculation. It is immaterial.
total_height += vd[index];

// Now we will find out the average height of the 8 neighbor particles.
// So we will know where the current particle will be attracted to.

// "heights" is the sum of all the height values of neighbor particles.
float heights = 0;
// "num_of_part" is the number of particles which contributed to the "heights".
int num_of_part = 0;

//// UP
if (!(index >= 0 && index < size)) // Make sure that this point is not on a boundary.
{
if (!vd_static[index - size]) // Make sure that the neighbor particle is not static.
{
heights += vd[index - size];

num_of_part += 1;
}
}

//// UPPER-RIGHT
if (!((index + 1) % size == 0 || (index >= 0 && index < size)))
{
if (!vd_static[index - size + 1])
{
heights += vd[index - size + 1];

num_of_part += 1;
}
}

//// RIGHT
if (!((index + 1) % size == 0))
{
if (!vd_static[index + 1])
{
heights += vd[index + 1];

num_of_part += 1;
}
}

//// LOWER-RIGHT
if (!((index + 1) % size == 0 || (index >= (size * size) - size && index < (size * size))))
{
if (!vd_static[index + size + 1])
{
heights += vd[index + size + 1];

num_of_part += 1;
}
}

//// DOWN
if (!(index >= (size * size) - size && index < (size * size)))
{
if (!vd_static[index + size])
{
heights += vd[index + size];

num_of_part += 1;
}
}

//// LOWER-LEFT
if (!(index % size == 0 || (index >= (size * size) - size && index < (size * size))))
{
if (!vd_static[index + size - 1])
{
heights += vd[index + size - 1];

num_of_part += 1;
}
}

//// LEFT
if (!(index % size == 0))
{
if (!vd_static[index - 1])
{
heights += vd[index - 1];

num_of_part += 1;
}
}

// UPPER-LEFT

if (!(index % size == 0 || (index >= 0 && index < size)))
{
if (!vd_static[index - size - 1])
{
heights += vd[index - size - 1];

num_of_part += 1;
}
}

if (num_of_part != 0)
{
heights /= (float)num_of_part;

if (power != 1.0f)
vda[index] += Math.Sign(heights - vd[index]) * (float)Math.Pow(Math.Abs(vd[index] - heights), power) / mass;
else
vda[index] += -(vd[index] - heights) / mass;
}

// Damping takes place.
vda[index] -= vdv[index] / vds[index];

// Don't let things go beyond their limit.
// This makes sense. It eliminates a critic uncertainty.
if (vda[index] > limit)
vda[index] = limit;
else if (vda[index] < -limit)
vda[index] = -limit;

}
// Now we have finished with the force calculation.

// Origin height is zero. So "shifting" is the distance between the system average height and the origin.
float shifting = -total_height / (float)vd.Length;

// We are taking the final steps in this loop
for (int index = 0; index < vd.Length; index += 1)
{

// Acceleration feeds velocity. Don't forget that we took care of the damping before.
vdv[index] += vda[index];

// Here is the purpose of "action_resolution":
// It is used to divide movements.
// If the particle goes along the road at once, a chaos is most likely unavoidable.
if (vd[index] + vdv[index] / action_resolution > limit)
vd[index] = limit;
else if (vd[index] + vdv[index] / action_resolution <= limit && vd[index] + vdv[index] / action_resolution >= -limit)
vd[index] += vdv[index] / action_resolution; // Velocity feeds height.
else
vd[index] = -limit;

// Here is the last step on shifting the whole system to the origin point.
vd[index] += shifting;

}

}

void generatebitmap()
{
if (bmp == null || bmp.Width != size)
bmp = new Bitmap(size, size, PixelFormat.Format24bppRgb); // 24 bit RGB. No need for alpha channel.

// Get the bitmap data of "bmp".
bd = bmp.LockBits(new Rectangle(0, 0, size, size), System.Drawing.Imaging.ImageLockMode.WriteOnly, System.Drawing.Imaging.PixelFormat.Format24bppRgb);

IntPtr ptr = bd.Scan0; // Get the address of the first line in "bd"
int bytes = bd.Stride * bd.Height; // "Stride" gives the size of a line in bytes.
byte[] rgbdata = new byte[bytes];

// It's time for the coloration of the height.
for (int index = 0; index < vd.Length; index++)
{
// Brightness. This value is the 'brightness' of the height.
// Now we see why "limit" makes sense.
byte bright = (byte)(((float)vd[index] + limit) / (float)((limit * 2f) / 255f));

if (vd_static[index])
{
rgbdata[index * 3] = ColorStatic.B;
rgbdata[index * 3 + 1] = ColorStatic.G;
rgbdata[index * 3 + 2] = ColorStatic.R;
}
else
{

if (highcont)
{

byte red_average = (byte)((float)(color1.R + color2.R) / 2f);
byte green_average = (byte)((float)(color1.G + color2.G) / 2f);
byte blue_average = (byte)((float)(color1.B + color2.B) / 2f);

if (vd[index] > 0)
{
rgbdata[index * 3] = color1.B;
rgbdata[index * 3 + 1] = color1.G;
rgbdata[index * 3 + 2] = color1.R;
}

else if (vd[index] < 0)
{
rgbdata[index * 3] = color2.B;
rgbdata[index * 3 + 1] = color2.G;
rgbdata[index * 3 + 2] = color2.R;
}
else if (vd[index] == 0)
{
rgbdata[index * 3] = blue_average;
rgbdata[index * 3 + 1] = green_average;
rgbdata[index * 3 + 2] = red_average;
}
}
else
{
float brightr1 = (float)bright / 255f;
float brightr2 = 1f - (float)bright / 255f;
rgbdata[index * 3] = (byte)((float)color1.B * brightr1 + (float)color2.B * brightr2);
rgbdata[index * 3 + 1] = (byte)((float)color1.G * brightr1 + (float)color2.G * brightr2);
rgbdata[index * 3 + 2] = (byte)((float)color1.R * brightr1 + (float)color2.R * brightr2);
}

}
}

// At last, we overwrite and release the bitmap data.
System.Runtime.InteropServices.Marshal.Copy(rgbdata, 0, ptr, bytes);
bmp.UnlockBits(bd);
}

/// <summary>
/// Sets sustainability of each particle.
/// </summary>
void setSustain()
{
if (edge_absorbtion)
{
// We will fill "vds" array with "sustain" then we will deal with elements near to window boundaries.

// Since we want the sustainability to decrease towards the edges, "min_sustain" can't be bigger than "sustain".
if (min_sustain > sustain)
{
min_sustain = 1.0f; // even "sustain" can't be less than 1.0f so this is a reliable value.
}

// Sustainability reduction fields should not mix with each other. So the maximum offset is the middle-screen.
if (absorb_offset >= size / 2)
{
absorb_offset = size / 2 - 1;
}

// This value is sustainability decreasion rate per row/column. The decreasion is linear.
float dec = (sustain - min_sustain) / (float)absorb_offset;
// This one stores the current sustainability.
float cur = min_sustain;

// First, we fill "vds" array with "sustain".
for (int i = 0; i < vds.Length - 1; i++)
vds[i] = sustain;

// This loop sets up the sustainability values for the top.
for (int off = 0; off <= absorb_offset; off++)
{
// Process each row/column from the edge to the offset.
for (int x = off; x < size - off; x++)
{
// Process each sustainability element in the current row/column
vds[x + off * size] = cur;
}
cur += dec;
}

cur = sustain; // Reset the current sustainability.

// This loop sets up the sustainability values for the bottom.
for (int off = 0; off <= absorb_offset; off++)
{
for (int x = absorb_offset - off; x < size - (absorb_offset - off); x++)
{
vds[x + off * size + size * (size - absorb_offset - 1)] = cur;
}
cur -= dec;
}

cur = sustain;

// This loop sets up the sustainability values for the left.
for (int off = 0; off <= absorb_offset; off++)
{
for (int x = absorb_offset - off; x < size - (absorb_offset - off); x++)
{
vds[x * size + (absorb_offset - off)] = cur;
}
cur -= dec;
}

cur = sustain;

// This loop sets up the sustainability values for the right.
for (int off = 0; off <= absorb_offset; off++)
{
for (int x = absorb_offset - off; x < size - (absorb_offset - off); x++)
{
vds[x * size + off + size - absorb_offset - 1] = cur;
}
cur -= dec;
}
}
else
{
// The only thing to do is to fill "vds" array with "sustain" in this case.
for (int i = 0; i < vds.Length; i++)
vds[i] = sustain;
}
}

/// <summary>
/// Initializes the wave pool system.
/// </summary>
void setPool()
{

if (bufgraph != null)
bufgraph.Dispose();

if (bufgcont != null)
bufgcont.Dispose();

bufgcont = new BufferedGraphicsContext();

bufgraph = bufgcont.Allocate(control.CreateGraphics(), control.ClientRectangle);

vd = new float[size * size];

vdv = new float[size * size];

vda = new float[size * size];

vd_static = new bool[size * size];

vds = new float[size * size];

setSustain();

}

}
}```

Now here is how we could use it in a `Form`. Let's start with this.

```using System;
using System.Drawing;
using System.Windows.Forms;

namespace WaveSimulator
{
public partial class MainForm : Form
{
public MainForm()
{
InitializeComponent();
}
}
}```

Declare an instance of the class in the global scope.

`WaveEngine we;`

Init `we` during the initialization. Specify `this` in the parameter.

`we = new WaveEngine(this);`

Now we may run it where ever we want it to run. Let's run it just after the initialization.

`we.Run();`

We also need to dispose the `we` while terminating the program. The program will be still running after the termination otherwise. The following code works for me.

```private void MainForm_FormClosing(object sender, FormClosingEventArgs e)
{
we.Dispose();
}```

That's it. Execute the program. You will see a window like this.

Here is the final code.

```using System;
using System.Drawing;
using System.Windows.Forms;

namespace WaveSimulator
{
public partial class MainForm : Form
{
WaveEngine we;

public MainForm()
{
InitializeComponent();

we = new WaveEngine(this);
we.Run();
}

private void MainForm_FormClosing(object sender, FormClosingEventArgs e)
{
we.Dispose();
}
}
}```

There is not any action currently. We need to change an attribute of a particle in order to see some action. We may simply achieve this by activating an oscillator. At most two oscillators can operate simultaneously.

We could click on the window to activate an oscillator at the mouse position. This one would work fine.

```private void MainForm_MouseDown(object sender, MouseEventArgs e)
{

we.Oscillator1Position = new Point((int)(((float)e.X / (float)ClientSize.Width) * (float)we.Size), (int)(((float)e.Y / (float)ClientSize.Height) * (float)we.Size));

we.Oscillator1Active = true;

}```

Execute the program and click anywhere on the window. You will see the waves. To make it a bit more interesting, add the following command after initializing `we`.

`we.EdgeAbsorbtion = false;`

Waves should bounce back from boundaries of the window now.

Since the class and its commentary are pretty self-explanatory, we will not go through the every single step. You may download the source code or demo. The project fully uses the class. Beside the main window, it contains the following.

An `Adjustments` window to control all of the properties. There is also a `1DView` window to see any cross-sectional area of the pool. It provides a side view of the pool.

## Points of Interest

When I was about to start making this, I thought it wouldn't work or it would just barely work but I was surprised when I realized that how simple its way of working could be. I was thinking about springs which could pull every single particle to its origin point but it was failure. Something was always going wrong. Most of the time I encountered an annoying scene like the following. I call it "chaos" as it fairly feels that way.

Reason of this situation was positive feedback due to low movement resolution but it took too much time for me to realize that. Two particles were first attracting each other and gaining a lot of acceleration and were moving so sharp. Finally they had more distance between each other than before. They had to gain greater acceleration in this case. So they were getting far away from each other in every movement instead of getting close. This leads to the chaos.

I was also having a problem with something simple-looking but yet hard-to-overcome. It was absorbtion of waves to prevent reflection. In a pool with an homogenous sustainability, waves simply bounce back from the edges of the window. It was first fun to see but soon I wanted to eliminate it. I tried to make particles near the edge static but it didn't work a bit. I tried a lot but all of them failed. After some time, I really tried to think so hard on this: "What could absorb waves? Friction? Of course not ju... wait... wait a minute... why not friction? High and smooth fricton on boundaries... Yes! Why not to give it a try?..". It worked but absorbtion parameters were highly dependent to particles' attributes.

There is another annoying thing that I'm yet to understand. When one particle attracts the other eight neighbor particles and also the particle itself is attracted by them, we need to consider the distance between each particle. Right? Unfortunately, considering the distance makes things go wrong somehow.

Above is an illustration which shows both possible distances between neighbor particles. One of them is 1 unit and the other is square root of 2 unit. Regarding to distance, particles which are square-root-of-two units away from each other should attract each other less while particles which are one unit away from each other should attract each other more so distance is inversely proportional to attraction force. Taking this situation into account makes things go wrong. I'm dividing the resultant force by the distance but it doesn't actually work.

## History

1 - Deleted `QueueUserWorkItem` function from properties and many functions. They used to use threads because I didn't want these operations to suspend the UI but these threaded operations often leads to unexpected results.

```WaveEngine we = new WaveEngine();

we.Size = 100;

int x = we.Size;```

What is `x` after all? It is not 100. If this code is executed without any interruption, `x` will not be equal to 100. This is because of the threaded property setting. The thread can't finish the operation before we retrieve `Size`. If one wants not to suspend the UI, a user-provided thread is necessary.

2 - Deleted `SwitchStaticDynamic` function. `SetParticles` function will be used instead.

3 - Added "Fixity" to `ParticleAttribute` enum. This makes it possible to use `SetParticles` function to convert particles between dynamic and static.

Minor fixes and improvements.

Important Note: I was astounded by todays modern computers' CPU power when I tried my program on one of them. My own computer is almost 10 years old with 1 GB RAM and 2.01 GHz Dualcore AMD CPU. I tried my program on a desktop computer equipped with Intel i7 3.4 GHz Quadcore CPU. It was an extremely fast simulation compared to the simulation speed in my computer. It was very fast but I wanted to slow it down with a smooth FPS. I adjusted the `Delay` but it didn't give me a smooth FPS. I finally increased the `Size` so number of particles has been increased. This increased the load on CPU and finally the rendering became smooth and slow. `Size` was great and the speed was still satisfying. So increase `Size` if your CPU is very fast and/or you think the simulation is very fast. Another option is to increase the movement resolution.

## Share

 Turkey
I'm a materials engineer, interested in simulations, programming, physics, biology, and space.

## You may also be interested in...

 First Prev Next
 My vote of 5 gregunger4-Mar-16 13:32 gregunger 4-Mar-16 13:32
 Re: My vote of 5 Mustafa Sami Salt15-Mar-16 11:11 Mustafa Sami Salt 15-Mar-16 11:11
 My vote of 5 D V L9-Sep-15 0:58 D V L 9-Sep-15 0:58
 Re: My vote of 5 Mustafa Sami Salt11-Sep-15 6:19 Mustafa Sami Salt 11-Sep-15 6:19
 My vote of 5 Yvan Rodrigues27-May-15 12:13 Yvan Rodrigues 27-May-15 12:13
 Re: My vote of 5 Mustafa Sami Salt27-May-15 23:06 Mustafa Sami Salt 27-May-15 23:06
 Distance going wrong FrankSzendzielarz26-May-15 11:36 FrankSzendzielarz 26-May-15 11:36
 Re: Distance going wrong Mustafa Sami Salt20-Jun-15 22:10 Mustafa Sami Salt 20-Jun-15 22:10
 My vote of 4 Santhakumar Munuswamy @ Chennai24-May-15 7:33 Santhakumar Munuswamy @ Chennai 24-May-15 7:33
 Re: My vote of 4 Mustafa Sami Salt24-May-15 8:34 Mustafa Sami Salt 24-May-15 8:34
 I like it Kenneth Haugland24-May-15 6:58 Kenneth Haugland 24-May-15 6:58
 Re: I like it Mustafa Sami Salt24-May-15 7:23 Mustafa Sami Salt 24-May-15 7:23
 Re: I like it Kenneth Haugland24-May-15 9:36 Kenneth Haugland 24-May-15 9:36
 My vote of 4 bobfox24-May-15 6:48 bobfox 24-May-15 6:48
 Re: My vote of 4 Mustafa Sami Salt24-May-15 7:04 Mustafa Sami Salt 24-May-15 7:04
 Last Visit: 31-Dec-99 18:00     Last Update: 23-Sep-18 19:38 Refresh 1