#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},
{
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]);
}
}
}
spatialFrequencyForIndex(m));
}
}
+
+ fourierAmplitudeAt(0, 0) = {0, 0};
}