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{1};
  10 const double Sea::GRAVITATIONAL_CONSTANT{6.67408e-11};
  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
  23 void Sea::update()
  24 {
  25   for (int y = 0; y < m_surface->size(); ++y) {
  26     for (int x = 0; x < m_surface->size(); ++x) {
  27       m_surface->at(x, y)
  28         .setHeight(((double)std::rand()/(double)RAND_MAX));
  29     }
  30   }
  31 }
  32
  33 double Sea::phillipsSpectrum(double k_x, double k_y) const
  34 {
  35   const double k = sqrt(pow(k_x, 2) + pow(k_y, 2));
  36   const double L = pow(m_windSpeed, 2) / GRAVITATIONAL_CONSTANT;
  37
  38   const double cosineFactor = pow((k_x / k) * m_windDirection[0] +
  39                                   (k_y / k) * m_windDirection[1], 2);
  40
  41   return PHILLIPS_CONSTANT * exp(-1 / pow(k * L, 2)) / pow(k, 4) *
  42     cosineFactor;
  43 }
  44
  45 ComplexPair Sea::generateFourierAmplitude(double k_x, double k_y)
  46 {
  47   std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
  48                              m_normalDistribution(m_randomGenerator));
  49
  50   return {cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2),
  51       std::conj(cDist) * sqrt(phillipsSpectrum(-k_x, -k_y)) / sqrt(2)};
  52 }
  53
  54 ComplexPair& Sea::fourierAmplitudeAt(int n, int m)
  55 {
  56   return m_fourierAmplitudes.at
  57     (n + m_surface->size()/2 +
  58      (m + m_surface->size()/2) * m_surface->size());
  59 }
  60
  61 double Sea::spatialFrequencyForIndex(int n) const
  62 {
  63   return 2 * M_PI * n / m_surface->size();
  64 }
  65
  66 void Sea::generateFourierAmplitudes()
  67 {
  68   for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
  69     for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
  70       fourierAmplitudeAt(n, m) =
  71         generateFourierAmplitude(spatialFrequencyForIndex(n),
  72                                  spatialFrequencyForIndex(m));
  73     }
  74   }
  75 }