]> git.treefish.org Git - seamulator.git/blob - sea.cpp
Basically it's working
[seamulator.git] / 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 }