From e19dd20d04a2a92b3d91bd25a69210a2351b08c9 Mon Sep 17 00:00:00 2001 From: Alexander Schmidt Date: Sun, 17 Jul 2016 00:34:21 +0200 Subject: [PATCH] Implemented generation of fourier amplitudes --- complexpair.h | 3 +++ sea.cpp | 58 +++++++++++++++++++++++++++++++++++++++++++++++- sea.h | 19 ++++++++++++++++ seamulator.cpp | 6 ++--- watersurface.cpp | 12 +++++----- watersurface.h | 4 ++-- 6 files changed, 91 insertions(+), 11 deletions(-) create mode 100644 complexpair.h diff --git a/complexpair.h b/complexpair.h new file mode 100644 index 0000000..9eea2a4 --- /dev/null +++ b/complexpair.h @@ -0,0 +1,3 @@ +#pragma once + +using ComplexPair = std::pair, std::complex>; diff --git a/sea.cpp b/sea.cpp index 42772d4..30593f2 100644 --- a/sea.cpp +++ b/sea.cpp @@ -1,11 +1,23 @@ #include "sea.h" +#include #include +#include #include "watersurface.h" -Sea::Sea(WaterSurfacePtr surface) : m_surface{surface} +const double Sea::PHILLIPS_CONSTANT{1}; +const double Sea::GRAVITATIONAL_CONSTANT{6.67408e-11}; + +Sea::Sea(WaterSurfacePtr surface) : + m_surface{surface}, + m_windDirection{1, 0}, + m_windSpeed{10}, + m_randomGenerator{m_randomDevice()}, + m_normalDistribution{0.0, 1.0} { + m_fourierAmplitudes.resize(pow(m_surface->size(), 2)); + generateFourierAmplitudes(); } void Sea::update() @@ -17,3 +29,47 @@ void Sea::update() } } } + +double Sea::phillipsSpectrum(double k_x, double k_y) const +{ + const double k = sqrt(pow(k_x, 2) + pow(k_y, 2)); + const double L = pow(m_windSpeed, 2) / GRAVITATIONAL_CONSTANT; + + const double cosineFactor = pow((k_x / k) * m_windDirection[0] + + (k_y / k) * m_windDirection[1], 2); + + return PHILLIPS_CONSTANT * exp(-1 / pow(k * L, 2)) / pow(k, 4) * + cosineFactor; +} + +ComplexPair Sea::generateFourierAmplitude(double k_x, double k_y) +{ + std::complex 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) +{ + return m_fourierAmplitudes.at + (n + m_surface->size()/2 + + (m + m_surface->size()/2) * m_surface->size()); +} + +double Sea::spatialFrequencyForIndex(int n) const +{ + return 2 * M_PI * n / m_surface->size(); +} + +void Sea::generateFourierAmplitudes() +{ + for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) { + for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) { + fourierAmplitudeAt(n, m) = + generateFourierAmplitude(spatialFrequencyForIndex(n), + spatialFrequencyForIndex(m)); + } + } +} diff --git a/sea.h b/sea.h index 48307ee..6e24ac3 100644 --- a/sea.h +++ b/sea.h @@ -1,7 +1,11 @@ #pragma once +#include +#include + #include "seafwd.h" +#include "complexpair.h" #include "watersurfacefwd.h" class Sea @@ -11,5 +15,20 @@ class Sea void update(); private: + static const double PHILLIPS_CONSTANT; + static const double GRAVITATIONAL_CONSTANT; + WaterSurfacePtr m_surface; + double m_windDirection[2]; + double m_windSpeed; + std::random_device m_randomDevice; + std::mt19937 m_randomGenerator; + std::normal_distribution<> m_normalDistribution; + std::vector m_fourierAmplitudes; + + double phillipsSpectrum(double k_x, double k_y) const; + ComplexPair generateFourierAmplitude(double k_x, double k_y); + ComplexPair& fourierAmplitudeAt(int n, int m); + void generateFourierAmplitudes(); + double spatialFrequencyForIndex(int n) const; }; diff --git a/seamulator.cpp b/seamulator.cpp index 565b99a..1c1ff67 100644 --- a/seamulator.cpp +++ b/seamulator.cpp @@ -9,12 +9,12 @@ #include "watersurface.h" const int LATTICE_SIZE{10}; -const double LATTICE_UNIT{1}; +const double LATTICE_EXTEND{10}; const int INIT_WINDOW_POS_X{100}; const int INIT_WINDOW_POS_Y{100}; const int INIT_WINDOW_WIDTH{300}; const int INIT_WINDOW_HEIGHT{300}; -const double INIT_VIEW_DISTANCE{LATTICE_SIZE * LATTICE_UNIT * 1.5}; +const double INIT_VIEW_DISTANCE{LATTICE_EXTEND * 1.5}; const double INIT_VIEW_AZIMUTH{0}; const double INIT_VIEW_ALTITUDE{M_PI / 4}; const char WINDOW_TITLE[]{"seamulator"}; @@ -32,7 +32,7 @@ int main(int argc, char** argv) { std::srand(std::time(0)); - surface = std::make_shared(LATTICE_SIZE, LATTICE_UNIT); + surface = std::make_shared(LATTICE_SIZE, LATTICE_EXTEND); sea = std::make_shared(surface); seaView = std::make_unique(INIT_VIEW_DISTANCE, INIT_VIEW_AZIMUTH, INIT_VIEW_ALTITUDE); diff --git a/watersurface.cpp b/watersurface.cpp index c27e7db..a38ad57 100644 --- a/watersurface.cpp +++ b/watersurface.cpp @@ -2,9 +2,9 @@ #include -WaterSurface::WaterSurface(int size, double unitLength) : +WaterSurface::WaterSurface(int size, double extend) : m_size{size}, - m_unitLength{unitLength} + m_extend{extend} { m_points.resize(size*size); } @@ -24,14 +24,16 @@ int WaterSurface::size() const return m_size; } -double WaterSurface::unitLength() const +double WaterSurface::extend() const { - return m_unitLength; + return m_extend; } void WaterSurface::draw() const { - glScalef(m_unitLength, m_unitLength, 1.0f); + const double scaleFactor{m_extend / m_size}; + + glScalef(scaleFactor, scaleFactor, 1.0f); glTranslatef(-(float)(m_size - 1) / 2, -(float)(m_size - 1) / 2, 0); for (int y = 0; y < m_size - 1; ++y) { diff --git a/watersurface.h b/watersurface.h index 31e0e60..22d3c48 100644 --- a/watersurface.h +++ b/watersurface.h @@ -13,12 +13,12 @@ class WaterSurface SurfacePoint& at(int x, int y); const SurfacePoint& at(int x, int y) const; int size() const; - double unitLength() const; + double extend() const; void draw() const; void drawSingleTile(int x, int y) const; private: std::vector m_points; int m_size; - double m_unitLength; + double m_extend; }; -- 2.39.5