Skip to content

FFT: wrap up

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;
			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);
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;

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!

FFT: Bit Reversing and butterfly

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).

y_0 = x_0 + x_1 \,

y_1 = x_0 - x_1. \,

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.

Fast Fourier Transform: Math

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


560558_10151845490857299_251833683_nNatz 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

a mac and web build will be up soon 😀

implamentaion: changing to unity due to visual studio not working :@

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;
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);

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;


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;


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;
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);
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;
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_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

implamentaion: evaluate waves using our height displacement and normal function

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 !

implementation: evaluateing the wave height, displacement and normal.

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 :

Kx = 2πn







but because we are using indices like a lot of the formulas we have to change it ao it now turns into :

Kx =









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:

and now we just add the height complex we made before:

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;