diff --git a/shaders/surface.vert b/shaders/surface.vert
index 4f8e1ad6dd95aed93a2ec1f3a27dde78efdb66bc..c855b14abbb4c42dd58fe144879581eb54bc6595 100644
--- a/shaders/surface.vert
+++ b/shaders/surface.vert
@@ -13,6 +13,9 @@ uniform int in_WavesNum;
 out vec3 world_pos;
 out vec3 normal;
 
+const float pi = 3.14159265358979323846;
+const float g = 9.82;
+
 struct Wave {
     float L; // Wavelength (distance between waves in world space)
     float A; // Amplitude (height from surface to wave crest)
@@ -24,18 +27,21 @@ layout(std140) uniform WaveBuffer {
     Wave waves[16];
 };
 
-/* Get the offset for a single wave */
-float get_gerstner1(Wave wave, float x, float z, float t) {
-    float w = 2.0/wave.L; // Frequency
-    float phase_const = wave.S * w; // Phase constant
+/* Get the offset for a single wave along the y-axis. */
+float get_gerstner_y(Wave wave, float x, float z, float t) {
+    // float w = 2.0/wave.L;
+    float w = sqrt(g*(2*pi/wave.L));
+    float phase_const = wave.S * w;
 
     return wave.A * sin(dot(normalize(wave.D) * w, vec2(x, z)) + t * phase_const);
 }
 
-float get_gerstner2(Wave wave, float dir, float x, float z, float t) {
-    float w = 2.0/wave.L; // Frequency
-    float phase_const = wave.S * w; // Phase constant
-    float Q = 0.3/(w*wave.A);
+/* Get the offset for a single wave in the xz-plane. */
+float get_gerstner_xz(Wave wave, float dir, float x, float z, float t) {
+    // float w = 2.0/wave.L;
+    float w = sqrt(g*(2*pi/wave.L));
+    float phase_const = wave.S * w;
+    float Q = 0.3/(w*wave.A*in_WavesNum);
     
     return Q * wave.A * dir * cos(dot(normalize(wave.D) * w, vec2(x, z)) + t * phase_const);
 }
@@ -45,10 +51,10 @@ vec3 get_waves(const vec3 pos, const float t) {
     vec3 offset = vec3(0.0, 0.0, 0.0);
 
     for (int i = 0; i < in_WavesNum; i++) {
-        offset.y += get_gerstner1(waves[i], pos.x, pos.z, t);
+        offset.y += get_gerstner_y(waves[i], pos.x, pos.z, t);
 
-        offset.x += get_gerstner2(waves[i], normalize(waves[i].D).x, pos.x, pos.z, t);
-        offset.z += get_gerstner2(waves[i], normalize(waves[i].D).y, pos.x, pos.z, t);
+        offset.x += get_gerstner_xz(waves[i], normalize(waves[i].D).x, pos.x, pos.z, t);
+        offset.z += get_gerstner_xz(waves[i], normalize(waves[i].D).y, pos.x, pos.z, t);
     }
     
     return offset + vec3(pos.x, 0.0, pos.z);
@@ -74,11 +80,9 @@ vec3 compute_normal(vec3 pos, float t) {
 
 void main(void)
 {
-    //vec3 offset_pos = vec3(in_Position.x, in_Position.y + 0.05 * sin(10.0 * in_Position.x + time), in_Position.z);
     vec3 offset_pos = get_waves(in_Position, time);
     world_pos = offset_pos;
     normal = compute_normal(in_Position, time);
-    // normal = in_Normal;
 
 	gl_Position = projectionMatrix * modelToWorldToView * vec4(offset_pos, 1.0);
 }
diff --git a/src/main.cpp b/src/main.cpp
index d7cd13ac12ea62fcbf5ba16ad1b2646e664a837d..72f1b11d2f5add71e3f6dfbbaa7878090110e5bd 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -11,6 +11,7 @@
 #include "VectorUtils4.h"
 #include "waves.h"
 #include <vector>
+#include <random>
 // uses framework OpenGL
 // uses framework Cocoa
 
@@ -33,9 +34,13 @@ struct Object {
 
 
 struct Surface : public Object {
-    static int const MAX_WAVES {16};
     static GLuint const BINDING_POINT {0};
     char const* UNIFORM_BLOCK = "WaveBuffer";
+    static int const MAX_WAVES {16};
+    static constexpr float const MEDIAN_WAVELENGTH {2.6};
+    static constexpr float const MEDIAN_AMPLITUDE {0.002};
+    static constexpr float const MEDIAN_SPEED {0.3};
+    static constexpr float const MEDIAN_DIR {3.14/4.0};
 
     struct Wave {
         GLfloat L; // Wavelength (distance between waves in world space)
@@ -48,10 +53,43 @@ struct Surface : public Object {
 
     GLuint ubo;
     GLuint block_index;
+
     std::vector<Wave> waves {};
+    // Direction of the wind [0-2π]
+    float wind_dir {3.14/3.0};
 
     Surface() = default;
 
+    Wave get_wave() {
+        std::random_device rd;
+        std::mt19937 gen(rd());
+
+        std::uniform_real_distribution<float> rd_wavelength(
+            0.5 * MEDIAN_WAVELENGTH, 2 * MEDIAN_WAVELENGTH
+        );
+
+        std::uniform_real_distribution<float> rd_amplitude(
+            0.5 * MEDIAN_AMPLITUDE, 2 * MEDIAN_AMPLITUDE
+        );
+
+        std::uniform_real_distribution<float> rd_speed(
+            0.5 * MEDIAN_SPEED, 2 * MEDIAN_SPEED
+        );
+
+        std::uniform_real_distribution<float> rd_dir(
+            wind_dir - MEDIAN_DIR, wind_dir + MEDIAN_DIR
+        );
+
+        float const dir {rd_dir(rd)};
+        return {
+            rd_wavelength(rd),
+            rd_amplitude(rd),
+            rd_speed(rd),
+            {},
+            vec2(cos(dir), sin(dir))
+        };
+    }
+
     Surface(Model* model, GLuint program)
         : Object{model, program} {
 
@@ -74,7 +112,9 @@ struct Surface : public Object {
         // // Link the UBO to the shader’s uniform block
         glUniformBlockBinding(program, block_index, BINDING_POINT);
 
-        waves.push_back({0.3f, 0.1f, 0.157f, {}, vec2{1.0f, 0.0f}});
+        for (int i = 0; i < MAX_WAVES; i++) {
+            waves.push_back(get_wave());
+        }
     }
 
     void draw() const