X-Git-Url: http://git.treefish.org/~alex/seamulator.git/blobdiff_plain/31e7aa23c4e95b5e37f33938e9bb47ea82a3f4d9..f2e44f70d56b60e6eef78656ca4b558917474276:/sea.cpp

diff --git a/sea.cpp b/sea.cpp
index 42772d4..85191d4 100644
--- a/sea.cpp
+++ b/sea.cpp
@@ -1,19 +1,135 @@
 #include "sea.h"
 
+#include <cmath>
 #include <cstdlib>
+#include <iostream>
 
 #include "watersurface.h"
 
-Sea::Sea(WaterSurfacePtr surface) : m_surface{surface}
+const double Sea::PHILLIPS_CONSTANT{0.0000001};
+const double Sea::GRAVITATIONAL_CONSTANT{9.8};
+
+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() + 1, 2));
+  generateFourierAmplitudes();
+
+  m_fftwIn = (fftw_complex*)
+    fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
+  m_fftwOut = (fftw_complex*)
+    fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
+  m_fftwPlan = fftw_plan_dft_2d
+    (m_surface->size(), m_surface->size(), m_fftwIn, m_fftwOut,
+     FFTW_BACKWARD, FFTW_MEASURE);
+
+  m_startTime = std::chrono::system_clock::now();
+}
+
+Sea::~Sea()
 {
+  fftw_destroy_plan(m_fftwPlan);
+  fftw_free(m_fftwIn); fftw_free(m_fftwOut);
+}
+
+double Sea::getRuntime() const
+{
+  auto timeNow = std::chrono::system_clock::now();
+  auto durationMs =
+    std::chrono::duration_cast<std::chrono::milliseconds>(timeNow - m_startTime);
+
+  return durationMs.count() / 1000.0;
 }
 
 void Sea::update()
 {
+  using namespace std::complex_literals;
+
+  const double runtime = getRuntime();
+
+  for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
+    const int positiveM = (m + m_surface->size()) % m_surface->size();
+
+    for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
+      const double k = sqrt(pow(spatialFrequencyForIndex(n), 2) +
+			    pow(spatialFrequencyForIndex(m), 2));
+      const double omega = sqrt(GRAVITATIONAL_CONSTANT * k);
+
+      std::complex<double> amplitude =
+       	fourierAmplitudeAt(n, m) * exp(1i * omega * runtime) +
+       	std::conj(fourierAmplitudeAt(-n, -m)) * exp(-1i * omega * runtime);
+
+      const int positiveN = (n + m_surface->size()) % m_surface->size();
+      int fftwIndex = positiveM + positiveN * m_surface->size();
+
+      m_fftwIn[fftwIndex][0] = std::real(amplitude);
+      m_fftwIn[fftwIndex][1] = std::imag(amplitude);
+    }
+  }
+
+  fftw_execute(m_fftwPlan);
+
   for (int y = 0; y < m_surface->size(); ++y) {
     for (int x = 0; x < m_surface->size(); ++x) {
       m_surface->at(x, y)
-	.setHeight(((double)std::rand()/(double)RAND_MAX));
+	.setHeight(m_fftwOut[y + x * m_surface->size()][0]);
+    }
+  }
+}
+
+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;
+}
+
+std::complex<double>& 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) {
+    const double k_y = spatialFrequencyForIndex(m);
+
+    for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
+      const double k_x = spatialFrequencyForIndex(n);
+
+      std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
+				 m_normalDistribution(m_randomGenerator));
+
+      fourierAmplitudeAt(n, m) =
+	cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2);
     }
   }
+
+  for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
+    fourierAmplitudeAt(n, m_surface->size()/2) =
+      fourierAmplitudeAt(n, -m_surface->size()/2);
+  }
+
+  for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
+    fourierAmplitudeAt(m_surface->size()/2, m) =
+      fourierAmplitudeAt(-m_surface->size()/2, m);
+  }
+
+  fourierAmplitudeAt(0, 0) = {0, 0};
 }