now all i have to do is show you how i changed my Evaluate waves most stuff is very similar to the previous evaluate waves but has some more code to help implement the FFT. Also due to me converting it into c# in unity some things also had to be changed.
firstly i had to make two index variables the original only had one which was index into vertices but i need one for buffers aswell so

index  = m_prime * N + n_prime;	// index into buffers
index1 = m_prime * Nplus1 + n_prime;	// index into vertices


also i have got rid of my Height_Displacement_and_Normal function and implemented it into this function ( found it worked better)

so at the top i basically just placed the function up there but it has also changed but it still very similar to its previous edition

for (int m_prime = 0; m_prime < N; m_prime++)
{
kz = Mathf.PI * (2.0f * m_prime - N) / m_length;

for (int n_prime = 0; n_prime < N; n_prime++)
{
kx = Mathf.PI*(2 * n_prime - N) / m_length;
len = Mathf.Sqrt(kx * kx + kz * kz);
index = m_prime * N + n_prime;

Vector2 c = InitSpectrum(t, n_prime, m_prime);

m_heightBuffer[1,index].x = c.x;
m_heightBuffer[1,index].y = c.y;

m_slopeBuffer[1,index].x = -c.y*kx;
m_slopeBuffer[1,index].y = c.x*kx;

m_slopeBuffer[1,index].z = -c.y*kz;
m_slopeBuffer[1,index].w = c.x*kz;

if (len < 0.000001f)
{
m_displacementBuffer[1,index].x = 0.0f;
m_displacementBuffer[1,index].y = 0.0f;
m_displacementBuffer[1,index].z = 0.0f;
m_displacementBuffer[1,index].w = 0.0f;
}
else
{
m_displacementBuffer[1,index].x = -c.y * -(kx/len);
m_displacementBuffer[1,index].y = c.x * -(kx/len);
m_displacementBuffer[1,index].z = -c.y * -(kz/len);
m_displacementBuffer[1,index].w = c.x * -(kz/len);
}
}
}


and after we call the Height_Displacement_and_Normal function and before we we do the evaluate waves we have to call the performfft function on the arrays of vecotr 2 and 3s we made:

m_fourier.PeformFFT(0, m_heightBuffer);
m_fourier.PeformFFT(0, m_slopeBuffer);
m_fourier.PeformFFT(0, m_displacementBuffer);


but i now cant call the function due to it not existing anymore so i had to make arrays to store the into as you can see above with m_displacementBuffer, m_slopeBuffer and m_heightBuffer

so setting up our normal, height and displacements are a little different to last time

// height
vertices[index1].y = m_heightBuffer[1, index].x * sign;

// displacement
vertices[index1].x = m_position[index1].x + m_displacementBuffer[1, index].x * lambda * sign;
vertices[index1].z = m_position[index1].z + m_displacementBuffer[1, index].z * lambda * sign;

// normal
n = new Vector3(-m_slopeBuffer[1, index].x * sign, 1.0f, -m_slopeBuffer[1, index].z * sign);
n.Normalize();

normals[index1].x =  n.x;
normals[index1].y =  n.y;
normals[index1].z =  n.z;



and the three if statements are also a little different

if (n_prime == 0 && m_prime == 0)
{
vertices[index1 + N + Nplus1 * N].y = m_heightBuffer[1, index].x * sign;

vertices[index1 + N + Nplus1 * N].x = m_position[index1 + N + Nplus1 * N].x + m_displacementBuffer[1, index].x * lambda * sign;
vertices[index1 + N + Nplus1 * N].z = m_position[index1 + N + Nplus1 * N].z + m_displacementBuffer[1, index].z * lambda * sign;

normals[index1 + N + Nplus1 * N].x =  n.x;
normals[index1 + N + Nplus1 * N].y =  n.y;
normals[index1 + N + Nplus1 * N].z =  n.z;
}
if (n_prime == 0)
{
vertices[index1 + N].y = m_heightBuffer[1, index].x * sign;

vertices[index1 + N].x = m_position[index1 + N].x + m_displacementBuffer[1, index].x * lambda * sign;
vertices[index1 + N].z = m_position[index1 + N].z + m_displacementBuffer[1, index].z * lambda * sign;

normals[index1 + N].x =  n.x;
normals[index1 + N].y =  n.y;
normals[index1 + N].z =  n.z;
}
if (m_prime == 0)
{
vertices[index1 + Nplus1 * N].y = m_heightBuffer[1, index].x * sign;

vertices[index1 + Nplus1 * N].x = m_position[index1 + Nplus1 * N].x + m_displacementBuffer[1, index].x * lambda * sign;
vertices[index1 + Nplus1 * N].z = m_position[index1 + Nplus1 * N].z + m_displacementBuffer[1, index].z * lambda * sign;

normals[index1 + Nplus1 * N].x =  n.x;
normals[index1 + Nplus1 * N].y =  n.y;
normals[index1 + Nplus1 * N].z =  n.z;
}


and after all this we just apply our fft formula to the mesh we made:

m_mesh.vertices = vertices;
m_mesh.normals = normals;
m_mesh.RecalculateBounds();


and were all done! it now all works. better than i thought it would but there are still other techniques o could try and implement but for now ill leave it at this!

this part of my research project has to be the hardest ive done so far with little sources explaining how this works without me knowing much about complicated uni grade math. so this one required some help from people who have already done this, they explained how it works and gave me some source code and links to read.

a blogger by the username Scrawk gave me some code but couldn’t explain how it works so i had to search far and wide to fine a simple explanation on how bit reversing and butterfly tables work.

i discovered that most formulas of FFT have to at least do some type of Bit reversal. for the bit reversal i found this website which explains in great detail what bit reversal does and what it is

it basically does what it says and reverses bits example the binary number 110 will now become 011. there is a lot more than that but its irreverent to the research so i recommend reading reading the page if you want to know more why do we do Bit reversal in FFT formulas?

well lets look at this pic i found from this website this pic shows an example of the time domain decomposition used in the FFT. In this example, a 16 point signal is decomposed through four separate stages. The first stage breaks the 16 point signal into two signals each consisting of 8 points. The second stage decomposes the data into four signals of 4 points. This pattern continues until there are N signals composed of a single point. An interlaced decomposition is used each time a signal is broken in two, that is, the signal is separated into its even and odd numbered samples.

once you look at the structure it becomes clear (apparently). The decomposition is nothing more than a reordering of the samples in the signal

this pic shows  the rearrangement pattern required On the left, the sample numbers of the original signal are listed along with their binary equivalents. On the right, the rearranged sample numbers are listed, also along with their binary equivalents. binary numbers are the reversals of each other! so a bit reversal is a lot cheaper and easier to do.

for butterfly diagrams the best place i could find to find some information on it was Wikipedia.

basically what a butterfly is is a portion of the computation that combines the results of smaller  discrete Fourier transform (DFTs) into a larger DFT or vice versa.

the butterfly diagram is commonly used in the cooley-turkey algorithm where a DFT of size N is recursively broken down into smaller transforms of size M where r is the size of radix of the transform.   These smaller DFTs are then combined via size-r butterflies, which themselves are DFTs of size r that are performed m times on corresponding outputs of the sub-transforms . these dfts are then pre-multiplied by roots of unity (known as twiddle factors).

lets say we have a radix-2 Cooley–Tukey algorithm, the butterfly is simply a DFT of size-2 that takes two inputs (x0x1) (corresponding outputs of the two sub-transforms) and gives two outputs (y0y1) by the formula (not including twiddle factors).  If one draws the data-flow diagram for this pair of operations, the (x0x1) to (y0y1) lines cross and resemble the wings of a butterfly hence the name…

the  picture below will illustrate this The butterfly can also be used to improve the randomness of large arrays of partially random numbers, by bringing every 32 or 64 bit word into causal contact with every other word through a desired hashing algorithm, so that a change in any one bit has the possibility of changing all the bits in the large array.

after some studying i under stand bit reversals a lot better and butterfly a little more hopefully i will understand it more before project is due.

lets Start with the equation from previous posts that described the wave height at the location n and m in terms of our indices and we re arrange it so it comes: this makes it easy for our two-dimensional discrete Fourier transform to become two One-dimensional discrete Fourier transforms.

they should look something like this now we want the fast Fourier transform to generate a wave height field at discrete points. this is done with this equation: but as i look down these formulas they seem to change this formula to : the h” equation is then changed by letting Lx=N and Lz=M so it becomes: now we need h” to be split into sums for over the even indices and the sum over the odd indices provided that N>_ 2

this is done by: (please note that this is just one equation not three. the image just shows how it changes )

we can makes this equation shorter if we let in the last line then the equation now becomes : now we want to eliminate the -1x from the equation we can do this by making a new function called h”‘. we wont be getting rid of the -1z which will come back once we do a second transform. so we have taken a size N DFT and split it into two size N/2 DFTs where the second has been multiplied by a twiddle  factor. Additionally, we have observed that the value at x+N2 can be obtained by reusing the intermediate computations used at x .

this equation should make a butterfly diagram like this: now if we want to do one for N≥4 then we can split h′′′ a second time. making this butterfly diagram: few things to note. Our x values run from − (N/2) through (N/2) −1 , so we’ve used the identity, Tx+N2N=−TxN , in the diagrams above. Our implementation will translate our x values by +(N/2) , so we will eliminate the negative. Additionally, in our implementation we can precompute our values for T . Lastly, the input values in the diagrams above are in bit-reversed order

now that i have the math i will find out how to implement it and show you in my next post. Natz and Boltz is a 2.5D puzzle platformer made for windows and mac which is being developed for our end of school assignment at the Academy of Interactive Entertainment. In the game you play as a young girl named Natz and her mech named Boltz who are trying to escape a factory which boltz has been locked up in.  All the puzzles are timed based and not very complicated due to only having about 3 months to develop this game. Our team (Team Freefall) consist of 6 artist and 4 programmers who all worked together to make this wonderful game. I worked on the Camera to make it and the programming of the puzzles.  The camera isn’t like any normal side scroller camera where the camera is fixed to the player, my camera follow the player like its a balloon tied to the player this gives the game a little more depth  and gives the game a unique feel . this was done by simply lerping the cameras position once the player has moved a certain distance from the camera. The camera will also zoom in depending on which room the player is in  to help the player know where to go or see the puzzle more easier. Our game also came runner up at the Sydney Independent game  festival which was a huge honor for us seeing as we beat games which had been in development for about a year and ours at the time was only about a month and a half into ours. We also had the privileged to go to the Australian EB expo as pert of our schools booth to show of our game to the public and it was a huge success with the kids and the occasional adult who loved the art and the game play Some even came back to play it more!  You can play test our game now for windows at  the games website natznboltz.com

a mac and web build will be up soon 😀

Due to the amount of time i have left to do this project and the amount of problems i am having with visual studio

i have decided to change to to unity to be my engine instead of making my own. the code is basically the same just i cant use the complex library i made so i just have to use vector 2

in this post ill show how i converted my c++ to c# which isn’t hard and what i changed to make this work

ill start with some new variables that i have to make things easier in unity

public Material m_mat;
public int numGridsX = 2;
public int numGridsZ = 2;
public int N = 64;
public float length = 64;
public float waveAmp = 0.0002f; // phillips spectrum parameter -- affects heights of waves
public Vector2 windSpeed = new Vector2(32.0f,32.0f);

GameObject[] m_oceanGrid;
Mesh m_mesh;
int Nplus1;
Vector2 m_windDirection;
FourierCPU m_fourier;// THIS IS FOR WHEN WE DO THE FAST FOURIER TRANSFORM EQUATION WILL EXPLAIN LATER
Vector2[,] heightBuffer;
Vector4[,] slopeBuffer, displacementBuffer;
Vector2[] spectrum, spectrum_conj;
Vector3[] position;
float[] dispersionTable;
Texture2D fresnelLookUp;


my my constructor had to be changed for unity as well and now looks like this:

void Start()
{
Nplus1 = N+1;

fourier = new FourierCPU(N);//again this if for fft

windDirection = new Vector2(m_windSpeed.x, m_windSpeed.y);
windDirection.Normalize();

dispersionTable = new float[Nplus1*Nplus1];
//this is where the Dispersion formula for hTilde has been moved this makes it more efficient and get higher fps
for (int m_prime = 0; m_prime &lt; Nplus1; m_prime++)
{
for (int n_prime = 0; n_prime &lt; Nplus1; n_prime++)
{
int index = m_prime * Nplus1 + n_prime;
dispersionTable[index] = Dispersion(n_prime,m_prime);
}
}

heightBuffer = new Vector2[2,N*N];
slopeBuffer = new Vector4[2,N*N];
displacementBuffer = new Vector4[2,N*N];

spectrum = new Vector2[Nplus1*Nplus1];
spectrum_conj = new Vector2[Nplus1*Nplus1];
position = new Vector3[Nplus1*Nplus1];

mesh = MakeMesh(Nplus1);

oceanGrid = new GameObject[numGridsX*numGridsZ];
//this renders my ocean plain but the meshes of the plain are made later will explain how meshes are rendered  later in this post
for(int x = 0; x &lt; numGridsX; x++)
{
for(int z = 0; z &lt; numGridsZ; z++)
{
int idx = x + z * numGridsX;

m_oceanGrid[idx] = new GameObject("Ocean grid " + idx.ToString());
m_oceanGrid[idx].renderer.material = m_mat;
m_oceanGrid[idx].GetComponent&lt;MeshFilter&gt;().mesh = m_mesh;
m_oceanGrid[idx].transform.Translate(new Vector3(x * m_length - numGridsX*m_length/2, 0.0f, z * length - numGridsZ*length/2));
}
}

Random.seed = 0;

Vector3[] vertices = m_mesh.vertices;

for (int m_prime = 0; m_prime &lt; Nplus1; m_prime++)
{
for (int n_prime = 0; n_prime &lt; Nplus1; n_prime++)
{
int index = m_prime * Nplus1 + n_prime;

m_spectrum[index] = GetSpectrum( n_prime,  m_prime);

m_spectrum_conj[index] = GetSpectrum(-n_prime, -m_prime);
m_spectrum_conj[index].y *= -1.0f;

m_position[index].x = vertices[index].x =  n_prime * m_length/N;
m_position[index].y = vertices[index].y =  0.0f;
m_position[index].z = vertices[index].z =  m_prime * m_length/N;

}
}

m_mesh.vertices = vertices;
m_mesh.RecalculateBounds();

CreateFresnelLookUp();
}


dispersion pretty much looks the same just now we have to use some Mathf variables

float Dispersion(int n_prime, int m_prime)
{
float w_0 = 2.0f * Mathf.PI / 200.0f;
float kx = Mathf.PI * (2 * n_prime - N) / m_length;
float kz = Mathf.PI * (2 * m_prime - N) / m_length;
return Mathf.Floor(Mathf.Sqrt(GRAVITY * Mathf.Sqrt(kx * kx + kz * kz)) / w_0) * w_0;
}


same with the phillips fomrula

float PhillipsSpectrum(int n_prime, int m_prime)
{
Vector2 k = new Vector2(Mathf.PI * (2 * n_prime - N) / length, Mathf.PI * (2 * m_prime - N) / length);
float length  = k.magnitude;
if (k_length &lt; 0.000001f) return 0.0f;

float length2 = length  * length;
float length4 = length2 * length2;

k.Normalize();

float k_dot_w   = Vector2.Dot(k, windDirection);
float k_dot_w2  = k_dot_w * k_dot_w * k_dot_w * k_dot_w * k_dot_w * k_dot_w;

float w_length  = windSpeed.magnitude;
float L         = w_length * w_length / GRAVITY;
float L2        = L * L;

float damping   = 0.001f;
float l2        = L2 * damping * damping;

return waveAmp * Mathf.Exp(-1.0f / (length2 * L2)) / k_length4 * k_dot_w2 * Mathf.Exp(-length2 * l2);
}


i also have to redo the Gaussian Random number function

Vector2 GaussianRandomVariable()
{
float x1, x2, w;
do
{
x1 = 2.0f * Random.value - 1.0f;
x2 = 2.0f * Random.value - 1.0f;
w = x1 * x1 + x2 * x2;
}
while ( w &gt;= 1.0f );

w = Mathf.Sqrt((-2.0f * Mathf.Log(w)) / w);
return new Vector2(x1 * w, x2 * w);
}


now instead of a complex we now have to use vector 2

Vector2 hTilde_0(int n_prime, int m_prime)
{
Vector2 r = GaussianRandomVariable();
return r * Mathf.Sqrt(PhillipsSpectrum(n_prime, m_prime) / 2.0f);
}
[\code]
in hTilde i have decided to move the dispersion formula to a new location  and stored in in a float called dispersionTable

so the fomrmula now looks a little different

Vector2 hTilde(float t, int n_prime, int m_prime)
{
int index = m_prime * Nplus1 + n_prime;

float dis = dispersionTable[index] * t;

float cos = Mathf.Cos(dis);
float sin = Mathf.Sin(dis);

float c0a = spectrum[index].x*cos - spectrum[index].y*sin;
float c0b = spectrum[index].x*sin + spectrum[index].y*cos;

float c1a = spectrum_conj[index].x*cos - spectrum_conj[index].y*-sin;
float c1b = spectrum_conj[index].x*-sin + spectrum_conj[index].y*cos;

return new Vector2(c0a+c1a, c0b+c1b);
}


now the normals, height and displacements are at the moment being reimplemented  while i do the evaluate waves for the fast Fourier transform  so when i show you how to do that very soon when i do the FFT implementation which isn't much more

also in unity i have created a function that stores the Fresnel term in a texture that will be used to reflect the sky of the ocean surface to give a more ocean feel. i am aware that its more common to do this in a shader by using a approximation but i have yet to learn how to use unity's shader and i'm still not that good at it so i am just going to put in in a function.  After talking with a guy named scrawk (not sure if that's his real name...) he told me how to the full math for calculating the Fresnel saying i only need to supply the refractive index of a material and it will give the correct Fresnel values. and because the math is a bit complicated that using a texture as a look up table would be wise

so here is the code he gave me:

void CreateFresnelLookUp()
{
float nSnell = 1.34f; //Refractive index of water

m_fresnelLookUp = new Texture2D(512, 1, TextureFormat.Alpha8, false);
m_fresnelLookUp.filterMode = FilterMode.Bilinear;
m_fresnelLookUp.wrapMode = TextureWrapMode.Clamp;
m_fresnelLookUp.anisoLevel = 0;

for(int x = 0; x &lt; 512; x++)
{
float fresnel = 0.0f;
float costhetai = (float)x/511.0f;
float thetai = Mathf.Acos(costhetai);
float sinthetat = Mathf.Sin(thetai)/nSnell;
float thetat = Mathf.Asin(sinthetat);

if(thetai == 0.0f)
{
fresnel = (nSnell - 1.0f)/(nSnell + 1.0f);
fresnel = fresnel * fresnel;
}
else
{
float fs = Mathf.Sin(thetat - thetai) / Mathf.Sin(thetat + thetai);
float ts = Mathf.Tan(thetat - thetai) / Mathf.Tan(thetat + thetai);
fresnel = 0.5f * ( fs*fs + ts*ts );
}

m_fresnelLookUp.SetPixel(x, 0, new Color(fresnel,fresnel,fresnel,fresnel));
}

m_fresnelLookUp.Apply();

m_mat.SetTexture("_FresnelLookUp", m_fresnelLookUp);
}


this is how i made a mesh:

Mesh MakeMesh(int size)
{

Vector3[] vertices = new Vector3[size*size];
Vector2[] texcoords = new Vector2[size*size];
Vector3[] normals = new Vector3[size*size];
int[] indices = new int[size*size*6];

for(int x = 0; x &lt; size; x++)
{
for(int y = 0; y &lt; size; y++)
{
Vector2 uv = new Vector3( (float)x / (float)(size-1), (float)y / (float)(size-1) );
Vector3 pos = new Vector3(x, 0.0f, y);
Vector3 norm = new Vector3(0.0f, 1.0f, 0.0f);

texcoords[x+y*size] = uv;
vertices[x+y*size] = pos;
normals[x+y*size] = norm;
}
}

int num = 0;
for(int x = 0; x &lt; size-1; x++)
{
for(int y = 0; y &lt; size-1; y++)
{
indices[num++] = x + y * size;
indices[num++] = x + (y+1) * size;
indices[num++] = (x+1) + y * size;

indices[num++] = x + (y+1) * size;
indices[num++] = (x+1) + (y+1) * size;
indices[num++] = (x+1) + y * size;
}
}

Mesh mesh = new Mesh();

mesh.vertices = vertices;
mesh.uv = texcoords;
mesh.triangles = indices;
mesh.normals = normals;

return mesh;
}


and i think were done next post ill finally start implementing the fast Fourier transform stuff

let me start by giving the code first then explaining it.

so here is the code:

void Ocean::evaluateWaves(float t)

{
float lambda = -1.0;
int index;
vector2 x;
vector2 d;
complex_vector_normal h_d_and_n;
for (int m_prime = 0; m_prime &lt; N; m_prime++)

{
for (int n_prime = 0; n_prime &lt; N; n_prime++) {
index = m_prime * Nplus1 + n_prime;

x = vector2(vertices[index].x, vertices[index].z);

h_d_and_n = Height_Displacement_and_Normal(x, t);

vertices[index].y = h_d_and_n.h.a;

vertices[index].x = vertices[index].ox + lambda*h_d_and_n.D.x;
vertices[index].z = vertices[index].oz + lambda*h_d_and_n.D.y;

vertices[index].nx = h_d_and_n.n.x;
vertices[index].ny = h_d_and_n.n.y;
vertices[index].nz = h_d_and_n.n.z;

if (n_prime == 0 &amp;&amp; m_prime == 0) {
vertices[index + N + Nplus1 * N].y = h_d_and_n.h.a;

vertices[index + N + Nplus1 * N].x = vertices[index + N + Nplus1 * N].ox + lambda*h_d_and_n.D.x;
vertices[index + N + Nplus1 * N].z = vertices[index + N + Nplus1 * N].oz + lambda*h_d_and_n.D.y;

vertices[index + N + Nplus1 * N].nx = h_d_and_n.n.x;
vertices[index + N + Nplus1 * N].ny = h_d_and_n.n.y;
vertices[index + N + Nplus1 * N].nz = h_d_and_n.n.z;
}
if (n_prime == 0) {
vertices[index + N].y = h_d_and_n.h.a;

vertices[index + N].x = vertices[index + N].ox + lambda*h_d_and_n.D.x;
vertices[index + N].z = vertices[index + N].oz + lambda*h_d_and_n.D.y;

vertices[index + N].nx = h_d_and_n.n.x;
vertices[index + N].ny = h_d_and_n.n.y;
vertices[index + N].nz = h_d_and_n.n.z;
}
if (m_prime == 0) {
vertices[index + Nplus1 * N].y = h_d_and_n.h.a;

vertices[index + Nplus1 * N].x = vertices[index + Nplus1 * N].ox + lambda*h_d_and_n.D.x;
vertices[index + Nplus1 * N].z = vertices[index + Nplus1 * N].oz + lambda*h_d_and_n.D.y;

vertices[index + Nplus1 * N].nx = h_d_and_n.n.x;
vertices[index + Nplus1 * N].ny = h_d_and_n.n.y;
vertices[index + Nplus1 * N].nz = h_d_and_n.n.z;
}
}
}
}


the three if statements are done so we can tile this seamlessly .  We’ve defined the size of the grid as Nplus1 so when we evaluate the wave heights, we only evaluate those heights on the N by N grid.

The last row and column are reserved to equal the wave heights of the first row and column, respectively. If we were dealing only with points, we wouldn’t need this extra row and column, but we want to fill that row and column with primitives, so those vertices have been added to draw the triangles. The first corner is a special case for adding the primitives to the last corner.

Since the height of the last vertices match those of the first, we can fit the second tile next to the first, and it will appear endless without any funny looking tilling.

if you look in the code you’ll see that lambda is -1. this  so so that the waves peaks are sharpened and the troughs are broadened. is it was 1 then it would be the opposite and would look strange. this is the only reason it is -1. and that’s basically it not a huge entry but it’ll do. next well finally start rendering this to get some pretty stuff on our screens !

remember that long and complicated equation we did  for the wave height? We define k⃗ =(kx,kz)

well we are now implementing that as well as our displacement and normals.

lets create some variables shall we to make things easier.

complex height(0.0f, 0.0f);
vector2 Displacement(0.0f, 0.0f);
vector3 normal(0.0f, 0.0f, 0.0f);

k = vector2(kx, kz);
k_length = k.length();
k_dot_x = k * x;


kx and kz are defined as : but because we are using indices like a lot of the formulas we have to change it ao it now turns into : Lx and Lz are the length of the plain we are using.

so implementation will look something like this

for (int m_prime = 0; m_prime &lt; N; m_prime++)
{
kz = 2.0f * M_PI * (m_prime - N / 2.0f) / length;
for (int n_prime = 0; n_prime &lt; N; n_prime++)
{
kx = 2.0f * M_PI * (n_prime - N / 2.0f) / length;
}

}


now we have to do the exp(ik⃗ ⋅x⃗ )

which is the Euler’s formula again

c = complex(cos(k_dot_x), sin(k_dot_x));


now we just do the furmula like this:

htilde_c = hTilde(t, n_prime, m_prime) * c;


where n and m prime are using the Kx znd Kz and we are timing  them by the Euler’s formula:

height = height + htilde_c;


now we have to do the normals

which is done like this which in code looks like this :

normal = normal + vector3(-kx * htilde_c.b, 0.0f, -kz * htilde_c.b);

normal = (vector3(0.0f, 1.0f, 0.0f) - normal).unit();


now for the displacement which has the wave height formula in it towards the end so all we pretty much have to do is do the

k⃗  \ k  part of it and multiply it by the wave height:

Displacement = Displacement + vector2(kx / k_length * htilde_c.b, kz / k_length * htilde_c.b);
//now we just have to return all three at the end :
complex_vector_normal cvn;
cvn.h = height;
cvn.D = Displacement;
cvn.n = normal;
return cvn; and were done!


full code:

complex_vector_normal Ocean::Height_Displacement_and_Normal(vector2 x, float t)
{
complex height(0.0f, 0.0f);
vector2 Displacement(0.0f, 0.0f);
vector3 normal(0.0f, 0.0f, 0.0f);
complex c, res, htilde_c;
vector2 k;
float kx, kz, k_length, k_dot_x;

for (int m_prime = 0; m_prime &lt; N; m_prime++)
{
kz = 2.0f * M_PI * (m_prime - N / 2.0f) / length;
for (int n_prime = 0; n_prime &lt; N; n_prime++)
{
kx = 2.0f * M_PI * (n_prime - N / 2.0f) / length;
k = vector2(kx, kz);

k_length = k.length();
k_dot_x = k * x;

c = complex(cos(k_dot_x), sin(k_dot_x));
htilde_c = hTilde(t, n_prime, m_prime) * c;

height = height + htilde_c;

normal = normal + vector3(-kx * htilde_c.b, 0.0f, -kz * htilde_c.b);

if (k_length &lt; 0.000001) continue;
Displacement = Displacement + vector2(kx / k_length * htilde_c.b, kz / k_length * htilde_c.b);
}
}

normal = (vector3(0.0f, 1.0f, 0.0f) - normal).unit();

complex_vector_normal cvn;
cvn.h = height;
cvn.D = Displacement;
cvn.n = normal;
return cvn;
}


ifest this year was a great day for my team ( team free fall) seeing a wide range of games from creative puzzle games to a creepy horror game called echo.

one of my favorite games at i fest was a game for iPhone and Ipad called TDTDTD by SMG which mixed a typical tower defense game mixed with xcom soldiers which you could strategically  move around the map to help the towers you place around the map to hold of the creative fish like aliens from destroying your base. i highly recommend it to anyone who has a iPhone or Ipad.

other games such particualars by Seethrough studios which came first in ifest as the best game was also one of my favorites due to its art direction and was also very calming to play. I’m interested to see where this game goes in the future

our game a got a lot of excellent feedback about the art and game play which was pleasing. they loved it so much that we came second for best game at ifest which is a great honor seeing as we beat games that have been in production for about a year or more and ours has only been in development  for about 2 – 3 months. so you can expect our team was a little over the moon when we got given this prize.

we are now even more excited for if we go to the eb expo in Sydney with our school to present our game at the AIE booth  to get even more people playing with it and giving us feedback so we can improve our game even more.

i

okay now we have to implement this equation into your code as well as this one so lets begin

ξr and ξi are independent Gaussian random variables so i went looking on how to do a random Gaussian variable equation and found third part source code due to time i decided not to look into this equation much but if i get time i will look into it and do a post about it. the code goes like this

float uniformRandomVariable()
{
return (float)rand()/RAND_MAX;
}

complex gaussianRandomVariable()
{
float x1, x2, w;
do
{
x1 = 2.f * uniformRandomVariable() - 1.f;
x2 = 2.f * uniformRandomVariable() - 1.f;
w = x1 * x1 + x2 * x2;
}
while ( w &gt;= 1.f );
w = sqrt((-2.f * log(w)) / w);
return complex(x1 * w, x2 * w);
}


also we finally get to use the previous equations

Ph(k) is the Phillips spectrum

so if we want to we could do the equation like this :

first we do the Gaussian values

complex r = gaussianRandomVariable();


then we just type the equation

return r * sqrt(phillips(n_prime, m_prime) / 2.0f);


we have the /2.0f because the equation also has a 1 over  sqrt of 2 in it so we just put at the end seeing as 1/X is well x.

now we go to the second and slightly more messy equation. but we can make it a lot cleaner by using a different equation that does the same thing

which i found during my research  which is the Euler’s Formula.

this equation uses the dispersion w(k) which we also made before

lets make a little function that does the dispersion equation for us to make the code a little neater.

float dis = dispersion(n_prime, m_prime) * t;


we want to use Euler’s Formula with complex numbers so we have to use this equation :  we want to cos and sin  our dispersion  so lets do that

float cos_ = cos(dis);
float sin_ = sin(dis);


and now we just do a part of  the Euler formula just the sine and cosine

complex c0(cos_,  sin_);
complex c1(cos_, -sin_);


before we continue we must also make a struct to hold variable for our vertex ocean which this equation uses

struct vertex_ocean
{

GLfloat   a,   b,   c; // htilde0
GLfloat  _a,  _b,  _c; // htilde0mk conjugate

};


we then get the get the vertices from the struct  which are for the vertex buffer object which ill get to later on

int index = m_prime * Nplus1 + n_prime;

complex htilde0(vertices[index].a,  vertices[index].b);
complex htilde0mkconj(vertices[index]._a, vertices[index]._b);


now we just finish of the Euler formula and we are done for this part:

return htilde0 * c0 + htilde0mkconj*c1;


sorry if this part was difficult to read this has been the hardest part of my research assignment  (thus far) for me and required a lot of research I’ve added links to the Euler formula if you wish to look into it a bit more as well

full code:

struct vertex_ocean
{

GLfloat   a,   b,   c; // htilde0
GLfloat  _a,  _b,  _c; // htilde0mk conjugate

};

complex Ocean::hTilde_0(int n_prime, int m_prime)

{
complex r = gaussianRandomVariable();
return r * sqrt(phillips(n_prime, m_prime) / 2.0f);
}

complex Ocean::hTilde(float t, int n_prime, int m_prime)

{
int index = m_prime * Nplus1 + n_prime;

complex htilde0(vertices[index].a,  vertices[index].b);
complex htilde0mkconj(vertices[index]._a, vertices[index]._b);

float dis = dispersion(n_prime, m_prime) * t;

float cos_ = cos(dis);
float sin_ = sin(dis);

complex c0(cos_,  sin_);
complex c1(cos_, -sin_);

return htilde0 * c0 + htilde0mkconj*c1;
}


okay the Phillips spectrum is the formula that helps us give us waves of random heights by using factors such as wind strength and direction stronger wind bigger waves essentially

the formula is as so : the w⃗   stands for the direction of the wind and L id the highest possible wave that can be generated by the wind speed.

now we have to expand this formula so we can easily under stand this to be implemented to our indices on our plain. lets just do the  k first by doing a vector 2
vector2 k(M_PI * (2 * n_prime – N) / length,
M_PI * (2 * m_prime – N) / length);

simple enough also note that this was also done in the  dispersion equation

also in the header file we have to put two variables

float spectrumPara; which we use to effect the height of the waves and is pretty much the A in the formula.

vector2 wind; which we will use to represent the wind parameters.

okay now we just got to make a couple of variables to make the equations easier to read these are not necessary but would be wise  to put in.

float k_length  = k.length();
if (k_length &lt; 0.000001) return 0.0;

float k_length2 = k_length  * k_length;
float k_length4 = k_length2 * k_length2;

float k_dot_w   = k.unit() * wind.unit();
float k_dot_w2  = k_dot_w * k_dot_w;

float w_length  = wind.length();
float L         = w_length * w_length / g;
float L2        = L * L;

float damping   = 0.001;
float l2        = L2 * damping * damping;


now we just make the formula by:

return spectrumPara * exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2 * exp(-k_length2 * l2);


the formula isn’t really that hard once you find out how to expand it to make it easier to read

next ill do the h~′0 and the complex conjugate, h~′0∗

full code :

float Ocean::phillips(int n_prime, int m_prime) {
vector2 k(M_PI * (2 * n_prime - N) / length,
M_PI * (2 * m_prime - N) / length);
float k_length  = k.length();
if (k_length &lt; 0.000001) return 0.0;

float k_length2 = k_length  * k_length;
float k_length4 = k_length2 * k_length2;

float k_dot_w   = k.unit() * wind.unit();
float k_dot_w2  = k_dot_w * k_dot_w;

float w_length  = wind.length();
float L         = w_length * w_length / g;
float L2        = L * L;

float damping   = 0.001;
float l2        = L2 * damping * damping;

return spectrumPara * exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2 * exp(-k_length2 * l2);
}