Implemented generation of fourier amplitudes

[seamulator.git] / sea.cpp

diff --git a/sea.cpp b/sea.cpp
index 42772d4d5534ba23e27bfbb2020c0c883ef766e9..30593f211f94f54cbf473799db34e8ccd2aae80e 100644 (file)
--- a/sea.cpp
+++ b/sea.cpp
@@ -1,11 +1,23 @@
 #include "sea.h"
 
+#include <cmath>
 #include <cstdlib>
+#include <iostream>
 
 #include "watersurface.h"
 
-Sea::Sea(WaterSurfacePtr surface) : m_surface{surface}
+const double Sea::PHILLIPS_CONSTANT{1};
+const double Sea::GRAVITATIONAL_CONSTANT{6.67408e-11};
+
+Sea::Sea(WaterSurfacePtr surface) :
+  m_surface{surface},
+  m_windDirection{1, 0},
+  m_windSpeed{10},
+  m_randomGenerator{m_randomDevice()},
+  m_normalDistribution{0.0, 1.0}
 {
+  m_fourierAmplitudes.resize(pow(m_surface->size(), 2));
+  generateFourierAmplitudes();
 }
 
 void Sea::update()
@@ -17,3 +29,47 @@ void Sea::update()
     }
   }
 }
+
+double Sea::phillipsSpectrum(double k_x, double k_y) const
+{
+  const double k = sqrt(pow(k_x, 2) + pow(k_y, 2));
+  const double L = pow(m_windSpeed, 2) / GRAVITATIONAL_CONSTANT;
+
+  const double cosineFactor = pow((k_x / k) * m_windDirection[0] +
+                                 (k_y / k) * m_windDirection[1], 2);
+
+  return PHILLIPS_CONSTANT * exp(-1 / pow(k * L, 2)) / pow(k, 4) *
+    cosineFactor;
+}
+
+ComplexPair Sea::generateFourierAmplitude(double k_x, double k_y)
+{
+  std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
+                            m_normalDistribution(m_randomGenerator));
+
+  return {cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2),
+      std::conj(cDist) * sqrt(phillipsSpectrum(-k_x, -k_y)) / sqrt(2)};
+}
+
+ComplexPair& Sea::fourierAmplitudeAt(int n, int m)
+{
+  return m_fourierAmplitudes.at
+    (n + m_surface->size()/2 +
+     (m + m_surface->size()/2) * m_surface->size());
+}
+
+double Sea::spatialFrequencyForIndex(int n) const
+{
+  return 2 * M_PI * n / m_surface->size();
+}
+
+void Sea::generateFourierAmplitudes()
+{
+  for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
+    for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
+      fourierAmplitudeAt(n, m) =
+       generateFourierAmplitude(spatialFrequencyForIndex(n),
+                                spatialFrequencyForIndex(m));
+    }
+  }
+}