]> git.treefish.org Git - seamulator.git/blob - src/sea.cpp
Fix compilation
[seamulator.git] / 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 Seamulator.  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 }