Basically it's working

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 63f645f322b95023485476aef13ab23b23b3fa3f..986edd410546422456b809c55a9484c70ce283fa 100644 (file)
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,12 +2,20 @@ cmake_minimum_required(VERSION 2.8)
 
 project(seamulator)
 
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14")
+
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}")
+
 find_package(OpenGL REQUIRED)
 find_package(GLUT REQUIRED)
-include_directories(${OPENGL_INCLUDE_DIRS} ${GLUT_INCLUDE_DIRS})
+find_package(fftw3 REQUIRED)
+
+include_directories(${OPENGL_INCLUDE_DIRS} ${GLUT_INCLUDE_DIRS}
+                   ${FFTW_INCLUDES})
 
 add_executable(seamulator seamulator.cpp
                          sea.cpp watersurface.cpp surfacepoint.cpp
                          seaview.cpp)
 
-target_link_libraries(seamulator ${OPENGL_LIBRARIES} ${GLUT_LIBRARY})
+target_link_libraries(seamulator ${OPENGL_LIBRARIES} ${GLUT_LIBRARY}
+                                ${FFTW_LIBRARIES} ${FFTW_LIBRARIES_THREADS})

diff --git a/Findfftw3.cmake b/Findfftw3.cmake
new file mode 100644 (file)
index 0000000..00c3401
--- /dev/null
+++ b/Findfftw3.cmake
@@ -0,0 +1,22 @@
+# - Find FFTW
+# Find the native FFTW includes and library
+#
+#  FFTW_INCLUDES    - where to find fftw3.h
+#  FFTW_LIBRARIES   - List of libraries when using FFTW.
+#  FFTW_FOUND       - True if FFTW found.
+
+if (FFTW_INCLUDES)
+  # Already in cache, be silent
+  set (FFTW_FIND_QUIETLY TRUE)
+endif (FFTW_INCLUDES)
+
+find_path (FFTW_INCLUDES fftw3.h)
+
+find_library (FFTW_LIBRARIES NAMES fftw3)
+
+# handle the QUIETLY and REQUIRED arguments and set FFTW_FOUND to TRUE if
+# all listed variables are TRUE
+include (FindPackageHandleStandardArgs)
+find_package_handle_standard_args (FFTW DEFAULT_MSG FFTW_LIBRARIES FFTW_INCLUDES)
+
+mark_as_advanced (FFTW_LIBRARIES FFTW_INCLUDES)

diff --git a/sea.cpp b/sea.cpp
index 30593f211f94f54cbf473799db34e8ccd2aae80e..b76a9ddf6a900eeb380d9cdc0077a51d020f18db 100644 (file)
--- a/sea.cpp
+++ b/sea.cpp
@@ -6,8 +6,8 @@
 
 #include "watersurface.h"
 
-const double Sea::PHILLIPS_CONSTANT{1};
-const double Sea::GRAVITATIONAL_CONSTANT{6.67408e-11};
+const double Sea::PHILLIPS_CONSTANT{0.0000003};
+const double Sea::GRAVITATIONAL_CONSTANT{9.8};
 
 Sea::Sea(WaterSurfacePtr surface) :
   m_surface{surface},
@@ -18,14 +18,68 @@ Sea::Sea(WaterSurfacePtr surface) :
 {
   m_fourierAmplitudes.resize(pow(m_surface->size(), 2));
   generateFourierAmplitudes();
+
+  fftw_init_threads();
+  fftw_plan_with_nthreads(4);
+
+  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).first * exp(1i * omega * runtime) +
+       fourierAmplitudeAt(n, m).second * 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]);
     }
   }
 }
@@ -72,4 +126,6 @@ void Sea::generateFourierAmplitudes()
                                 spatialFrequencyForIndex(m));
     }
   }
+
+  fourierAmplitudeAt(0, 0) = {0, 0};
 }

diff --git a/sea.h b/sea.h
index 6e24ac3bd3a5f859198ab6329015601c7bdb3d73..20f5d638d33479fd21d5b8daf0626e591d9926ac 100644 (file)
--- a/sea.h
+++ b/sea.h
@@ -1,8 +1,11 @@
 #pragma once
 
+#include <chrono>
 #include <complex>
 #include <random>
 
+#include <fftw3.h>
+
 #include "seafwd.h"
 
 #include "complexpair.h"
@@ -12,6 +15,9 @@ class Sea
 {
  public:
   Sea(WaterSurfacePtr surface);
+  ~Sea();
+  Sea(const Sea&) = delete;
+  Sea& operator=(const Sea&) = delete;
   void update();
 
  private:
@@ -25,10 +31,14 @@ class Sea
   std::mt19937 m_randomGenerator;
   std::normal_distribution<> m_normalDistribution;
   std::vector<ComplexPair> m_fourierAmplitudes;
+  fftw_complex *m_fftwIn, *m_fftwOut;
+  fftw_plan m_fftwPlan;
+  std::chrono::time_point<std::chrono::system_clock> m_startTime;
 
   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;
+  double getRuntime() const;
 };

diff --git a/seamulator.cpp b/seamulator.cpp
index 1c1ff671e9d29d91f3f25aad3a0ad2d00f93619e..c5f7936c652c9aed72afb23747476b15b6a9b457 100644 (file)
--- a/seamulator.cpp
+++ b/seamulator.cpp
@@ -8,12 +8,12 @@
 #include "seaview.h"
 #include "watersurface.h"
 
-const int LATTICE_SIZE{10};
+const int LATTICE_SIZE{128};
 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 int INIT_WINDOW_POS_X{50};
+const int INIT_WINDOW_POS_Y{50};
+const int INIT_WINDOW_WIDTH{800};
+const int INIT_WINDOW_HEIGHT{600};
 const double INIT_VIEW_DISTANCE{LATTICE_EXTEND * 1.5};
 const double INIT_VIEW_AZIMUTH{0};
 const double INIT_VIEW_ALTITUDE{M_PI / 4};