From: Alexander Schmidt <alex@treefish.org>
Date: Sat, 16 Jul 2016 22:34:21 +0000 (+0200)
Subject: Implemented generation of fourier amplitudes
X-Git-Url: http://git.treefish.org/~alex/seamulator.git/commitdiff_plain/e19dd20d04a2a92b3d91bd25a69210a2351b08c9

Implemented generation of fourier amplitudes
---

diff --git a/complexpair.h b/complexpair.h
new file mode 100644
index 0000000..9eea2a4
--- /dev/null
+++ b/complexpair.h
@@ -0,0 +1,3 @@
+#pragma once
+
+using ComplexPair = std::pair<std::complex<double>, std::complex<double>>;
diff --git a/sea.cpp b/sea.cpp
index 42772d4..30593f2 100644
--- 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));
+    }
+  }
+}
diff --git a/sea.h b/sea.h
index 48307ee..6e24ac3 100644
--- a/sea.h
+++ b/sea.h
@@ -1,7 +1,11 @@
 #pragma once
 
+#include <complex>
+#include <random>
+
 #include "seafwd.h"
 
+#include "complexpair.h"
 #include "watersurfacefwd.h"
 
 class Sea
@@ -11,5 +15,20 @@ class Sea
   void update();
 
  private:
+  static const double PHILLIPS_CONSTANT;
+  static const double GRAVITATIONAL_CONSTANT;
+
   WaterSurfacePtr m_surface;
+  double m_windDirection[2];
+  double m_windSpeed;
+  std::random_device m_randomDevice;
+  std::mt19937 m_randomGenerator;
+  std::normal_distribution<> m_normalDistribution;
+  std::vector<ComplexPair> m_fourierAmplitudes;
+
+  double phillipsSpectrum(double k_x, double k_y) const;
+  ComplexPair generateFourierAmplitude(double k_x, double k_y);
+  ComplexPair& fourierAmplitudeAt(int n, int m);
+  void generateFourierAmplitudes();
+  double spatialFrequencyForIndex(int n) const;
 };
diff --git a/seamulator.cpp b/seamulator.cpp
index 565b99a..1c1ff67 100644
--- a/seamulator.cpp
+++ b/seamulator.cpp
@@ -9,12 +9,12 @@
 #include "watersurface.h"
 
 const int LATTICE_SIZE{10};
-const double LATTICE_UNIT{1};
+const double LATTICE_EXTEND{10};
 const int INIT_WINDOW_POS_X{100};
 const int INIT_WINDOW_POS_Y{100};
 const int INIT_WINDOW_WIDTH{300};
 const int INIT_WINDOW_HEIGHT{300};
-const double INIT_VIEW_DISTANCE{LATTICE_SIZE * LATTICE_UNIT * 1.5};
+const double INIT_VIEW_DISTANCE{LATTICE_EXTEND * 1.5};
 const double INIT_VIEW_AZIMUTH{0};
 const double INIT_VIEW_ALTITUDE{M_PI / 4};
 const char WINDOW_TITLE[]{"seamulator"};
@@ -32,7 +32,7 @@ int main(int argc, char** argv)
 {
   std::srand(std::time(0));
 
-  surface = std::make_shared<WaterSurface>(LATTICE_SIZE, LATTICE_UNIT);
+  surface = std::make_shared<WaterSurface>(LATTICE_SIZE, LATTICE_EXTEND);
   sea = std::make_shared<Sea>(surface);
   seaView = std::make_unique<SeaView>(INIT_VIEW_DISTANCE, INIT_VIEW_AZIMUTH,
 				      INIT_VIEW_ALTITUDE);
diff --git a/watersurface.cpp b/watersurface.cpp
index c27e7db..a38ad57 100644
--- a/watersurface.cpp
+++ b/watersurface.cpp
@@ -2,9 +2,9 @@
 
 #include <GL/glut.h>
 
-WaterSurface::WaterSurface(int size, double unitLength) :
+WaterSurface::WaterSurface(int size, double extend) :
   m_size{size},
-  m_unitLength{unitLength}
+  m_extend{extend}
 {
   m_points.resize(size*size);
 }
@@ -24,14 +24,16 @@ int WaterSurface::size() const
   return m_size;
 }
 
-double WaterSurface::unitLength() const
+double WaterSurface::extend() const
 {
-  return m_unitLength;
+  return m_extend;
 }
 
 void WaterSurface::draw() const
 {
-  glScalef(m_unitLength, m_unitLength, 1.0f);
+  const double scaleFactor{m_extend / m_size};
+
+  glScalef(scaleFactor, scaleFactor, 1.0f);
   glTranslatef(-(float)(m_size - 1) / 2, -(float)(m_size - 1) / 2, 0);
 
   for (int y = 0; y < m_size - 1; ++y) {
diff --git a/watersurface.h b/watersurface.h
index 31e0e60..22d3c48 100644
--- a/watersurface.h
+++ b/watersurface.h
@@ -13,12 +13,12 @@ class WaterSurface
   SurfacePoint& at(int x, int y);
   const SurfacePoint& at(int x, int y) const;
   int size() const;
-  double unitLength() const;
+  double extend() const;
   void draw() const;
   void drawSingleTile(int x, int y) const;
 
  private:
   std::vector<SurfacePoint> m_points;
   int m_size;
-  double m_unitLength;
+  double m_extend;
 };