-#include "sea.h"
-
-#include <cmath>
-#include <cstdlib>
-#include <iostream>
-
-#include "watersurface.h"
-
-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(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};
-}