Wave Simulator with Java and C++

Here are the GitHub repositories: Java and C++

Introduction

This is a wave simulation engine which can be used to observe many wave-related phenomena. Diffraction, interference, reflection, refraction, etc. can be observed.

The engine is written in both Java and C++. A native Linux application is also written for the C++ version of the engine. The image above is a screenshot from the program. For the Java version, only a tutorial project is available.

There is also a C# version in the article titled Wave Simulator with C#. The C# version is currently very outdated compared to these versions.

Background

We are going to go through many concepts which are related to the simulation algorithm. Here is an image of an example pool.

Concepts

1. Pool

The simulation takes place inside a square area called pool. It is the entire region as shown in the image above.

2. Particle

The pool is made up of particles which are used to transmit energy across the pool. They have altitude, velocity, and mass. They are attracted to each other and want to stay in the same altitude that of neighbor particles. Note that the particles can never move laterally so the attraction force can only change their altitude.

The particles are similar to pixels on a monitor. If the size of a pool is 300x300 for example, then the number of particles in that pool will be 90000. They are also equally spaced like pixels.

There are different types of particles which behave differently.

2.1 Oscillator

An oscillator constantly vibrates at a certain frequency. It is not attracted to other particles and makes neighbor particles move up and down in the same frequency. The result is a sinusoidal wave propagating over the entire area. Two oscillators with different frequencies can be seen in the image above.

2.2 Static Particle

A static particle never moves. It is not attracted to neighbor particles and is not attractive itself. Obstacles, reflectors, cables and such can be made up of static particles. The horizontal yellow wall in the image above consists of static particles which are used to prevent waves of the two oscillators from mixing with each other.

3. Medium

The term medium is simply used to indicate a set of particles whose mass differs from that of the rest. A medium may be thought of as glass, air, water, etc. When waves enter a different medium, their speed and direction may change. There are two mediums as shown in the image above. They are heavier than the rest of the particles and hence slow down the waves and shorten their wavelength.

Calculation Steps

We are going to take a look at the most important function in the engine: calculateForces

In this function, forces exerted on each particle are calculated to determine the next position for each particle.

For each particle, the function does the following:

• Check if the particle is dynamic. Skip to the next particle otherwise.
• Sum the height of each dynamic neighbor particle. The variable which holds the sum is heights. The variable which holds the number of particles contributed to the sum is num_of_parts. For example, the following code checks if a dynamic neighbor in upper direction exists and adds its height to the sum accordingly.

if (index >= size && !vd_static[index - size]) {
            up_exists = true;
            heights += vd[index - size];
            num_of_parts++;
}

We are doing this to find the average height of the dynamic neighbor particles.

vd array holds the height of the particles, and vd_static array holds the information which determines if a particle is static or dynamic.

• Divide the heights by num_of_parts to obtain the average height of the dynamic neighbor particles which is stored in height_difference.
• Calculate the acceleration.

double acceleration = 0.0;
double height_difference = 0.0;

if (num_of_parts != 0) {
    heights /= num_of_parts;
    height_difference = vd[index] - heights;
    acceleration = -height_difference / vdm[index];
}

vdm array holds the mass of the particles.

As you can see, the Newton's second law of motion is used.

Force = Mass * Acceleration

Acceleration = Force / Mass

So the height_difference is directly considered to be a force. This is similar to the spring model wherein the Hooke's law states that the force needed to extend or compress the spring by some distance linearly changes with respect to that distance.

Force = Constant * Distance

Distance becomes height_difference in our model because the rest point of a particle is the average height of its neighbors. Constant is taken as 1.

• Add acceleration to the velocity of the particle.

  vdv[index] += acceleration;

  vdv array holds the velocity of the particles.

  Particles may lose some of their mechanical energy which is a sum of kinetic and potential energy over time depending on the efficiency of the system. There is a variable loss which is a ratio to determine the efficiency. The higher this value, the lower the efficiency. At each calculation, or iteration in another word, the mechanical energy is calculated as:

  Mechanical Energy *= (1 - Loss Ratio)
  
  Mechanical Energy = Kinetic Energy + Potential Energy

  Both the kinetic energy and potential energy are multiplied by the same value.

  The kinetic energy equation is as follows:

  Kinetic Energy = 0.5 * Mass * Velocity<sup>2</sup>

• Find the kinetic energy.

  double kinetic_energy = (0.5 * vdm[index] * pow(vdv[index], 2));

• Multiply the kinetic_energy by (1 - loss) and find the velocity for the resulting kinetic energy.

  vdv[index] = sqrt(2 * kinetic_energy * (1.0 - vdl[index]) / vdm[index]) * sgn(vdv[index]);

  vdl array holds the loss value for each particle.

  As velocity might have been a negative value, we multiplied the resulting velocity with the sign of the previous velocity.

  The potential energy equation is as follows:

  Potential Energy = 0.5 * Constant * Height Difference<sup>2</sup>

  This is originally the potential energy of a spring. Constant is again taken as 1 and Height Difference is the distance between the average height of the neighbor particles and the current particle's height.

• Find the potential energy.

  double potential_energy = (0.5 * pow(height_difference, 2.0));

• Multiply the potential energy with (1 - loss) and find the height for that potential energy.

  vd[index] += sqrt(2 * potential_energy * (1.0 - vdl[index])) * sgn(height_difference) 
               - height_difference;

  After this point, the oscillator-related calculations are done. The calculation involves putting a timestamp into a sine function getting amplitude as an output.

Using the Code

We are going to take a look at several methods, tips, and tricks. For more details, check out the manual PDF found in the data/help path of the wavesim_cpp repository.

Initialization

WaveEngine waveEngine = new WaveEngine();

The main thread and the co-threads are created and immediately put into a sleepy state. The engine doesn't paint or calculate anything until we call start(). It will stop as soon as we call stop().

Pool

Here are some configuration methods for the pool.

waveEngine.setSize(500);

Here, we specify the size of the pool as 500 in pixel units. Remember that the pool is always a square so both the width and height of the pool is set by this method.

The absorber can be enabled to eliminate echoes from the walls of the pool. Enabling the absorber increases the loss ratio close to the walls so that waves approaching the walls lose significant energy. The following code enables it.

waveEngine.setAbsorberEnabled(true);

Oscillators

The maximum number of oscillators which can run simultaneously is 9. This value can be changed by changing MAX_NUMBER_OF_OSCILLATORS. An unlimited number of oscillators could be supported but that would mean additional work in GUI part of the C++ version of the program.

The shortest way to define and activate an oscillator is as follows:

waveEngine.setOscillatorSource(0, OscillatorSource.PointSource);
waveEngine.setOscillatorLocation(0, 0, new Point(150, 150));
waveEngine.setOscillatorEnabled(0, true);

In the first line, we set the source type of the first oscillator as a point. It could also be a line or a moving point.

In the second line, we set the oscillator's primary location relative to the pool, not the window. For example, if the window in which the rendering takes place is 500x500 in size and the pool size is 300x300, then specifying (150,150) as a location would put the oscillator at the center of the pool.

In the last line, we set the oscillator as enabled so it emits energy.

The first parameter is always the index number of the particular oscillator. In setOscillatorLocation method and the other oscillator location related methods, the second parameter is the index of the location. 0 is primary and 1 is secondary. Secondary location is used for the line source or moving source wherein latter it defines the target point for the moving oscillator.

Many features of the oscillator like period, phase, amplitude, etc. can be changed. It is recommended to use a lower amplitude for line source oscillators as such the energy emitted from them are several magnitudes higher. With such an high energy, even the diffracted waves would be powerful and the directionality of line source oscillators wouldn't be realized.

Threading

The engine supports multi-threading. The main thread controls and gives missions to co-threads. Actual calculations and paintings are mostly done by the co-threads. The maximum number of co-threads is 32. This value can be changed by changing MAX_NUMBER_OF_THREADS. An unlimited number of threads could be supported but there were some problems which I couldn't deal with at that time.

We can set the number of co-threads as follows:

waveEngine.setNumberOfThreads(4);

Here, we specify the number of co-threads as 4. For the best performance, although not guaranteed, set this to the number of CPU cores. For example, if your CPU is a dual-core one, set it to 2, or, if it is a quad-core one, set it to 4.

Setting it to a lower value would cause less CPU usage. This is preferable if the iterations per second in runtime can reach the maximum value determined by the limiter. Now coming to the limiter, we can limit the number of iterations/calculation done every second by the following method:

waveEngine.setIterationsPerSecond(70);

Setting this to zero will disable the limiter and the co-threads will make calculations as fast as possible.

Rendering

All rendering/painting calculations are done on CPU. The graphical interface such as GDI or Cairo which are used to draw the final image may use GPU acceleration though but that is irrelevant to the engine's operation.

There is an FPS (frames per second) limiter which can be set by the following method.

waveEngine.setFramesPerSecond(30);

Here, we set the FPS to 30 so the engine will make 30 paintings every second. Note that the paintings are made by CPU, not GPU so increasing this value may result in serious performance issues. Setting this to zero will disable the limiter but is really not recommended.

To make the waves look stronger on the screen, increase the amplitude multiplier by the following method:

waveEngine.setAmplitudeMultiplier(10);

Here, we tell the engine to multiply the height values with 10 during the painting stage so the hard-to-see waves will look more clearly. This can also be thought of as adjusting the contrast.

waveEngine.setExtremeContrastEnabled(true);

The method above will basically set the amplitude multiplier to infinity.

The most important thing about rendering is the callback method. The RenderEventListener interface can be used to specify a callback method for rendering. Every time the engine makes a painting, it calls this method. It can be set as follows:

waveEngine.setRenderEventListener((byte[] bitmap_data) -> {
            System.arraycopy(bitmap_data, 0, data, 0, data.length);
            repaint();
 });

bitmap_data: Holds the 24-bit RGB value for each particle. This is the array that we would like to copy from.

data: This is the target buffer array which is used by an image such as BufferedImage. The data array could be obtained by the following code:

BufferedImage bufImage = new BufferedImage(size, size, BufferedImage.TYPE_3BYTE_BGR);
Raster raster = bufImage.getRaster();
DataBufferByte dataBuffer = (DataBufferByte) raster.getDataBuffer();
byte[] data = dataBuffer.getData();

data can be modified directly to make changes on the pixels of the bufImage.

One thing to note here is that the color format for the engine is 24-bit RGB while we selected color format for the bufImage is 24-bit BGR. We see that the red and blue components are swapped. In such a case, we could overcome the potential issues by swapping the red and blue components of the colors while setting the crest, trough, and static particles colors.

So in this case, if want to set the crest color to red, instead of setting it to (255,0,0), we set it to (0,0,255) so that it will look red as we desired. There is no method to change the color format of the engine from RGB to BGR.

Crest, trough, and static particle colors can be respectively set by the following methods:

waveEngine.setCrestColor(new Color((byte)255, (byte)0, (byte)0));
waveEngine.setTroughColor(new Color((byte)0, (byte)0, (byte)255));
waveEngine.setStaticColor(new Color((byte)0, (byte)255, (byte)0));

Here, we make the crests look red, trough look blue, and static particles look green.

Locking and Unlocking

To read/write particle attributes such as velocity and mass, we need to use lock() and unlock() methods. Here is an example code from the Java project.

// We are going to modify particle data which requires a call to lock().
// Between calls to lock() and unlock(), only getData() must be called,
// otherwise the thread which called another method will hang.
// Because of the reasons mentioned above we call getSize() before
// calling lock().
int size = waveEngine.getSize();
waveEngine.lock();

// Get the fixity data which determines if a particle is static or
// dynamic. There are size*size elements in the returned array which
// is the number of particles in the pool.
// The method returns an array whose type is determined by the
// ParticleAttribute. If it is Fixity, the type will be byte[],
// otherwise the type will be double[].
byte fixities[] = (byte[]) waveEngine.getData(ParticleAttribute.Fixity);

// Set a rectangle of static particles to the left.
for (int x = 10; x < 30; x++) {
      for (int y = 120; y < 220; y++) {
            fixities[x + size * y] = 1;
      }
}
        
// Let's also add a different medium to the top-right.
// Sadly, the different medium will be visible by only its effects on
// the incident waves (no colors involved).
double masses[] = (double[]) waveEngine.getData(ParticleAttribute.Mass);
        
// Set a rectangle of heavier particles to the top.
for (int x = 180; x < size; x++) {
      for (int y = 0; y < 150; y++) {
            masses[x + size * y] = 0.4;
      }
}
        
// After calling unlock(), do not read/write the data obtained
// by getData() call. Doing so will either result in a program error or
// calculation errors.
waveEngine.unlock();

Points of Interest

I learned a lot of things while working on these projects and I don't regret spending my time on them.

Some time ago, I decided to give Linux OS a try and installed & tried many distributions. I actually enjoyed my time and decided to dive deeper.

One day, I remembered that I had this Wave Simulator with C# article years ago, so why not writing a native Linux application which will be much more functional and educative for my programming and physics knowledge, I asked to myself. So I have made this engine with its program, learned tons of things about technologies that I have never known before such as gcc, g++, glade, cairo, gtk, pango, and many more. The user interface design with glade was fun except a few bugs. I was really impressed by the neat and robust nature of the UI format.

I also decided to write the Java version of the engine.

I really like waves and think that the waves are one of the most fundamental things in the universe.

History

• 20^th September, 2018: Initial version