src/sea.cpp

   1 /**
   2  * Copyright (C) 2016  Alexander Schmidt
   3  *
   4  * This file is part of Seamulator.
   5  *
   6  * Seamulator is free software: you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation, either version 3 of the License, or
   9  * (at your option) any later version.
  10  *
  11  * Seamulator is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  * GNU General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with Foobar.  If not, see <http://www.gnu.org/licenses/>.
  18  */
  19
  20 #include "sea.h"
  21
  22 #include <cmath>
  23 #include <cstdlib>
  24 #include <iostream>
  25
  26 #include "watersurface.h"
  27
  28 const double Sea::GRAVITATIONAL_CONSTANT{9.8};
  29
  30 Sea::Sea(WaterSurfacePtr surface, double windSpeed, double magicConstant) :
  31   m_surface{surface},
  32   m_windDirection{1, 0},
  33   m_windSpeed{windSpeed},
  34   m_magicConstant{magicConstant},
  35   m_randomGenerator{m_randomDevice()},
  36   m_normalDistribution{0.0, 1.0}
  37 {
  38   m_fourierAmplitudes.resize(pow(m_surface->size() + 1, 2));
  39   generateFourierAmplitudes();
  40
  41   m_fftwIn = (fftw_complex*)
  42     fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
  43   m_fftwOut = (fftw_complex*)
  44     fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
  45   m_fftwPlan = fftw_plan_dft_2d
  46     (m_surface->size(), m_surface->size(), m_fftwIn, m_fftwOut,
  47      FFTW_BACKWARD, FFTW_MEASURE);
  48
  49   m_startTime = std::chrono::system_clock::now();
  50 }
  51
  52 Sea::~Sea()
  53 {
  54   fftw_destroy_plan(m_fftwPlan);
  55   fftw_free(m_fftwIn); fftw_free(m_fftwOut);
  56 }
  57
  58 double Sea::getRuntime() const
  59 {
  60   auto timeNow = std::chrono::system_clock::now();
  61   auto durationMs =
  62     std::chrono::duration_cast<std::chrono::milliseconds>(timeNow - m_startTime);
  63
  64   return durationMs.count() / 1000.0;
  65 }
  66
  67 void Sea::update()
  68 {
  69   using namespace std::complex_literals;
  70
  71   const double runtime = getRuntime();
  72
  73   for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
  74     const int positiveM = (m + m_surface->size()) % m_surface->size();
  75
  76     for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
  77       const double k = sqrt(pow(spatialFrequencyForIndex(n), 2) +
  78                             pow(spatialFrequencyForIndex(m), 2));
  79       const double omega = sqrt(GRAVITATIONAL_CONSTANT * k);
  80
  81       std::complex<double> amplitude =
  82         fourierAmplitudeAt(n, m) * exp(1i * omega * runtime) +
  83         std::conj(fourierAmplitudeAt(-n, -m)) * exp(-1i * omega * runtime);
  84
  85       const int positiveN = (n + m_surface->size()) % m_surface->size();
  86       int fftwIndex = positiveM + positiveN * m_surface->size();
  87
  88       m_fftwIn[fftwIndex][0] = std::real(amplitude);
  89       m_fftwIn[fftwIndex][1] = std::imag(amplitude);
  90     }
  91   }
  92
  93   fftw_execute(m_fftwPlan);
  94
  95   for (int y = 0; y < m_surface->size(); ++y) {
  96     for (int x = 0; x < m_surface->size(); ++x) {
  97       m_surface->at(x, y)
  98         .setHeight(m_fftwOut[y + x * m_surface->size()][0]);
  99     }
 100   }
 101 }
 102
 103 double Sea::phillipsSpectrum(double k_x, double k_y) const
 104 {
 105   const double k = sqrt(pow(k_x, 2) + pow(k_y, 2));
 106   const double L = pow(m_windSpeed, 2) / GRAVITATIONAL_CONSTANT;
 107
 108   const double cosineFactor = pow((k_x / k) * m_windDirection[0] +
 109                                   (k_y / k) * m_windDirection[1], 2);
 110
 111   return m_magicConstant * exp(-1 / pow(k * L, 2)) / pow(k, 4) *
 112     cosineFactor;
 113 }
 114
 115 std::complex<double>& Sea::fourierAmplitudeAt(int n, int m)
 116 {
 117   return m_fourierAmplitudes.at
 118     (n + m_surface->size()/2 +
 119      (m + m_surface->size()/2) * m_surface->size());
 120 }
 121
 122 double Sea::spatialFrequencyForIndex(int n) const
 123 {
 124   return 2 * M_PI * n / m_surface->size();
 125 }
 126
 127 void Sea::generateFourierAmplitudes()
 128 {
 129   for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
 130     const double k_y = spatialFrequencyForIndex(m);
 131
 132     for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
 133       const double k_x = spatialFrequencyForIndex(n);
 134
 135       std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
 136                                  m_normalDistribution(m_randomGenerator));
 137
 138       fourierAmplitudeAt(n, m) =
 139         cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2);
 140     }
 141   }
 142
 143   for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
 144     fourierAmplitudeAt(n, m_surface->size()/2) =
 145       fourierAmplitudeAt(n, -m_surface->size()/2);
 146   }
 147
 148   for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
 149     fourierAmplitudeAt(m_surface->size()/2, m) =
 150       fourierAmplitudeAt(-m_surface->size()/2, m);
 151   }
 152
 153   fourierAmplitudeAt(0, 0) = {0, 0};
 154 }