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