#include "watersurface.h"
-const double Sea::PHILLIPS_CONSTANT{0.0000003};
+const double Sea::PHILLIPS_CONSTANT{0.0000002};
const double Sea::GRAVITATIONAL_CONSTANT{9.8};
Sea::Sea(WaterSurfacePtr surface) :
m_randomGenerator{m_randomDevice()},
m_normalDistribution{0.0, 1.0}
{
- m_fourierAmplitudes.resize(pow(m_surface->size(), 2));
+ m_fourierAmplitudes.resize(pow(m_surface->size() + 1, 2));
generateFourierAmplitudes();
- fftw_init_threads();
- fftw_plan_with_nthreads(4);
-
m_fftwIn = (fftw_complex*)
fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
m_fftwOut = (fftw_complex*)
const double omega = sqrt(GRAVITATIONAL_CONSTANT * k);
std::complex<double> amplitude =
- fourierAmplitudeAt(n, m).first * exp(1i * omega * runtime) +
- fourierAmplitudeAt(n, m).second * exp(-1i * omega * runtime);
+ fourierAmplitudeAt(n, m) * exp(1i * omega * runtime) +
+ std::conj(fourierAmplitudeAt(-n, -m)) * exp(-1i * omega * runtime);
const int positiveN = (n + m_surface->size()) % m_surface->size();
int fftwIndex = positiveM + positiveN * m_surface->size();
cosineFactor;
}
-ComplexPair Sea::generateFourierAmplitude(double k_x, double k_y)
-{
- std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
- m_normalDistribution(m_randomGenerator));
-
- return {cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2),
- std::conj(cDist) * sqrt(phillipsSpectrum(-k_x, -k_y)) / sqrt(2)};
-}
-
-ComplexPair& Sea::fourierAmplitudeAt(int n, int m)
+std::complex<double>& Sea::fourierAmplitudeAt(int n, int m)
{
return m_fourierAmplitudes.at
(n + m_surface->size()/2 +
void Sea::generateFourierAmplitudes()
{
for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
+ const double k_y = spatialFrequencyForIndex(m);
+
for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
+ const double k_x = spatialFrequencyForIndex(n);
+
+ std::complex<double> cDist(m_normalDistribution(m_randomGenerator),
+ m_normalDistribution(m_randomGenerator));
+
fourierAmplitudeAt(n, m) =
- generateFourierAmplitude(spatialFrequencyForIndex(n),
- spatialFrequencyForIndex(m));
+ cDist * sqrt(phillipsSpectrum(k_x, k_y)) / sqrt(2);
}
}
+ for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
+ fourierAmplitudeAt(n, m_surface->size()/2) =
+ fourierAmplitudeAt(n, -m_surface->size()/2);
+ }
+
+ for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
+ fourierAmplitudeAt(m_surface->size()/2, m) =
+ fourierAmplitudeAt(-m_surface->size()/2, m);
+ }
+
fourierAmplitudeAt(0, 0) = {0, 0};
}
#include "seafwd.h"
-#include "complexpair.h"
#include "watersurfacefwd.h"
class Sea
std::random_device m_randomDevice;
std::mt19937 m_randomGenerator;
std::normal_distribution<> m_normalDistribution;
- std::vector<ComplexPair> m_fourierAmplitudes;
+ std::vector<std::complex<double>> m_fourierAmplitudes;
fftw_complex *m_fftwIn, *m_fftwOut;
fftw_plan m_fftwPlan;
std::chrono::time_point<std::chrono::system_clock> m_startTime;
double phillipsSpectrum(double k_x, double k_y) const;
- ComplexPair generateFourierAmplitude(double k_x, double k_y);
- ComplexPair& fourierAmplitudeAt(int n, int m);
+ std::complex<double>& fourierAmplitudeAt(int n, int m);
void generateFourierAmplitudes();
double spatialFrequencyForIndex(int n) const;
double getRuntime() const;