sea.cpp

   1 #include "sea.h"
   2
   3 #include <cmath>
   4 #include <cstdlib>
   5 #include <iostream>
   6
   7 #include "watersurface.h"
   8
   9 const double Sea::PHILLIPS_CONSTANT{0.0000003};
  10 const double Sea::GRAVITATIONAL_CONSTANT{9.8};
  11
  12 Sea::Sea(WaterSurfacePtr surface) :
  13   m_surface{surface},
  14   m_windDirection{1, 0},
  15   m_windSpeed{10},
  16   m_randomGenerator{m_randomDevice()},
  17   m_normalDistribution{0.0, 1.0}
  18 {
  19   m_fourierAmplitudes.resize(pow(m_surface->size(), 2));
  20   generateFourierAmplitudes();
  21
  22   fftw_init_threads();
  23   fftw_plan_with_nthreads(4);
  24
  25   m_fftwIn = (fftw_complex*)
  26     fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
  27   m_fftwOut = (fftw_complex*)
  28     fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
  29   m_fftwPlan = fftw_plan_dft_2d
  30     (m_surface->size(), m_surface->size(), m_fftwIn, m_fftwOut,
  31      FFTW_BACKWARD, FFTW_MEASURE);
  32
  33   m_startTime = std::chrono::system_clock::now();
  34 }
  35
  36 Sea::~Sea()
  37 {
  38   fftw_destroy_plan(m_fftwPlan);
  39   fftw_free(m_fftwIn); fftw_free(m_fftwOut);
  40 }
  41
  42 double Sea::getRuntime() const
  43 {
  44   auto timeNow = std::chrono::system_clock::now();
  45   auto durationMs =
  46     std::chrono::duration_cast<std::chrono::milliseconds>(timeNow - m_startTime);
  47
  48   return durationMs.count() / 1000.0;
  49 }
  50
  51 void Sea::update()
  52 {
  53   using namespace std::complex_literals;
  54
  55   const double runtime = getRuntime();
  56
  57   for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
  58     const int positiveM = (m + m_surface->size()) % m_surface->size();
  59
  60     for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
  61       const double k = sqrt(pow(spatialFrequencyForIndex(n), 2) +
  62                             pow(spatialFrequencyForIndex(m), 2));
  63       const double omega = sqrt(GRAVITATIONAL_CONSTANT * k);
  64
  65       std::complex<double> amplitude =
  66         fourierAmplitudeAt(n, m).first * exp(1i * omega * runtime) +
  67         fourierAmplitudeAt(n, m).second * exp(-1i * omega * runtime);
  68
  69       const int positiveN = (n + m_surface->size()) % m_surface->size();
  70       int fftwIndex = positiveM + positiveN * m_surface->size();
  71
  72       m_fftwIn[fftwIndex][0] = std::real(amplitude);
  73       m_fftwIn[fftwIndex][1] = std::imag(amplitude);
  74     }
  75   }
  76
  77   fftw_execute(m_fftwPlan);
  78
  79   for (int y = 0; y < m_surface->size(); ++y) {
  80     for (int x = 0; x < m_surface->size(); ++x) {
  81       m_surface->at(x, y)
  82         .setHeight(m_fftwOut[y + x * m_surface->size()][0]);
  83     }
  84   }
  85 }
  86
  87 double Sea::phillipsSpectrum(double k_x, double k_y) const
  88 {
  89   const double k = sqrt(pow(k_x, 2) + pow(k_y, 2));
  90   const double L = pow(m_windSpeed, 2) / GRAVITATIONAL_CONSTANT;
  91
  92   const double cosineFactor = pow((k_x / k) * m_windDirection[0] +
  93                                   (k_y / k) * m_windDirection[1], 2);
  94
  95   return PHILLIPS_CONSTANT * exp(-1 / pow(k * L, 2)) / pow(k, 4) *
  96     cosineFactor;
  97 }
  98
  99 ComplexPair Sea::generateFourierAmplitude(double k_x, double k_y)
 100 {
 101   std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
 102                              m_normalDistribution(m_randomGenerator));
 103
 104   return {cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2),
 105       std::conj(cDist) * sqrt(phillipsSpectrum(-k_x, -k_y)) / sqrt(2)};
 106 }
 107
 108 ComplexPair& Sea::fourierAmplitudeAt(int n, int m)
 109 {
 110   return m_fourierAmplitudes.at
 111     (n + m_surface->size()/2 +
 112      (m + m_surface->size()/2) * m_surface->size());
 113 }
 114
 115 double Sea::spatialFrequencyForIndex(int n) const
 116 {
 117   return 2 * M_PI * n / m_surface->size();
 118 }
 119
 120 void Sea::generateFourierAmplitudes()
 121 {
 122   for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
 123     for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
 124       fourierAmplitudeAt(n, m) =
 125         generateFourierAmplitude(spatialFrequencyForIndex(n),
 126                                  spatialFrequencyForIndex(m));
 127     }
 128   }
 129
 130   fourierAmplitudeAt(0, 0) = {0, 0};
 131 }