/**
 * Copyright (C) 2016  Alexander Schmidt
 *
 * This file is part of Seamulator.
 *
 * Seamulator is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Seamulator is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with Seamulator.  If not, see <http://www.gnu.org/licenses/>.
 */

#include "sea.h"

#include <cmath>
#include <cstdlib>
#include <iostream>

#include "watersurface.h"

const double Sea::GRAVITATIONAL_CONSTANT{9.8};

Sea::Sea(WaterSurfacePtr surface, double windSpeed, double magicConstant) :
  m_surface{surface},
  m_windDirection{1, 0},
  m_windSpeed{windSpeed},
  m_magicConstant{magicConstant},
  m_randomGenerator{m_randomDevice()},
  m_normalDistribution{0.0, 1.0}
{
  m_fourierAmplitudes.resize(pow(m_surface->size() + 1, 2));
  generateFourierAmplitudes();

  m_fftwIn = (fftw_complex*)
    fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
  m_fftwOut = (fftw_complex*)
    fftw_malloc(sizeof(fftw_complex) * pow(m_surface->size(), 2));
  m_fftwPlan = fftw_plan_dft_2d
    (m_surface->size(), m_surface->size(), m_fftwIn, m_fftwOut,
     FFTW_BACKWARD, FFTW_MEASURE);

  m_startTime = std::chrono::system_clock::now();
}

Sea::~Sea()
{
  fftw_destroy_plan(m_fftwPlan);
  fftw_free(m_fftwIn); fftw_free(m_fftwOut);
}

double Sea::getRuntime() const
{
  auto timeNow = std::chrono::system_clock::now();
  auto durationMs =
    std::chrono::duration_cast<std::chrono::milliseconds>(timeNow - m_startTime);

  return durationMs.count() / 1000.0;
}

void Sea::update()
{
  using namespace std::complex_literals;

  const double runtime = getRuntime();

  for (int m = -m_surface->size()/2; m < m_surface->size()/2; ++m) {
    const int positiveM = (m + m_surface->size()) % m_surface->size();

    for (int n = -m_surface->size()/2; n < m_surface->size()/2; ++n) {
      const double k = sqrt(pow(spatialFrequencyForIndex(n), 2) +
			    pow(spatialFrequencyForIndex(m), 2));
      const double omega = sqrt(GRAVITATIONAL_CONSTANT * k);

      std::complex<double> amplitude =
       	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();

      m_fftwIn[fftwIndex][0] = std::real(amplitude);
      m_fftwIn[fftwIndex][1] = std::imag(amplitude);
    }
  }

  fftw_execute(m_fftwPlan);

  for (int y = 0; y < m_surface->size(); ++y) {
    for (int x = 0; x < m_surface->size(); ++x) {
      m_surface->at(x, y)
	.setHeight(m_fftwOut[y + x * m_surface->size()][0]);
    }
  }
}

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 m_magicConstant * exp(-1 / pow(k * L, 2)) / pow(k, 4) *
    cosineFactor;
}

std::complex<double>& 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) {
    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) =
	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};
}