BasicOceanSystem.cpp

#include "Rendering/IGraphicsDevice.h"
#include "Rendering/ICapabilities.h"
#include "Foundation/Logging/Logger.h"
 
#include "BasicOceanSystem.h"
 
#include "Components/Appearance/MaterialComponent.h"
#include "Components/Geometry/AdaptivePlanarGridComponent.h"
#include "Components/Core/SceneComponent.h"
#include "Components/Core/CameraComponent.h"
#include "Components/Behavior/ReflectionComponent.h"
 
#include "Systems/Core/TransformSystem.h"
#include "Systems/Geometry/AdaptivePlanarGridSystem.h"
 
#include "Resources/MaterialManager.h"
#include "Resources/TextureManager.h"
#include "Resources/Texture.h"
 
#include "Services/Variables.h"
#include "Services/Time.h"
#include "Services/TaskManager.h"
 
#include "EntityStore.h"
 
#include "Utilities/Parsing.h"
 
#include "Context.h"
 
#define _USE_MATH_DEFINES
#include <cmath>
#include <glm/glm.hpp>
#include <algorithm>
#include <limits>
#include <numbers>
#include <random>
 
#include <cassert>
#if !defined(__EMSCRIPTEN__) && !defined(__APPLE__)
#include <xmmintrin.h>
#endif
 
namespace {
  using namespace Cogs::Core;
  Cogs::Logging::Log logger = Cogs::Logging::getLogger("BasicOceanSystem");
 
  const std::string animateKey = "basic-ocean.animate";
  const std::string eightbitKey = "basic-ocean.8bit";
 
  void initMaterialInstanceCallback(MaterialInstance* instance, Cogs::ComponentModel::Entity* /*container*/, void* data)
  {
    auto * oceanData = static_cast<BasicOceanData*>(data);
    assert(oceanData->waveVariant != nullptr);
    instance->setVariant("Encoding", oceanData->encoding);
    instance->setVariant("Waves", oceanData->waveVariant);
    instance->setVariant("BeliefSystem", oceanData->beliefSystem);
    instance->setVariant("LightModel", oceanData->lightModelVariant);
    instance->setVariant("Reflection", oceanData->reflectionVariant);
  }
 
  using std::complex;
 
  // Reverse the lower logN bits of i
  inline size_t reverseBits16(size_t i, size_t log2N)
  {
    assert(log2N <= 16);
 
    size_t t = i;
    t = ((0xAAAAu & t) >> 1) | ((t << 1) & 0xAAAAu);  // 1010101010101010
    t = ((0xCCCCu & t) >> 2) | ((t << 2) & 0xCCCCu);  // 1100110011001100
    t = ((0xF0F0u & t) >> 4) | ((t << 4) & 0xF0F0u);  // 1111000011110000
    t = ((0xFF00u & t) >> 8) | ((t << 8) & 0xFF00u);  // 1111111100000000
    return t >> (16 - log2N);
  }
 
#if 0
  const float twoPi = float(2.0 * M_PI);
  // Calculate complex twiddle-factor e^{-2\pi n/N}
  inline complex<float> twiddleFactor(const int n, const int N)
  {
    const float arg = (-twoPi * float(n)) / float(N);
    return complex<float>(cos(arg), sin(arg));
  }
#endif
 
 
#if !defined(__EMSCRIPTEN__) && !defined(__APPLE__)
  void inverseRadix2MajorIndexTransposeAlignedSoALog4SSE1(float* dstReal,
    float* dstImag,
    const float* srcReal,
    const float* srcImag,
    const float factor,
    const float scale,
    const size_t log2N)
  {
    size_t N = size_t(1) << log2N;
    size_t Nq = N / 4;
 
    // first pass, also handle transpose and shuffle 
    for (size_t j = 0; j < N / 2; j++) {
      size_t index0 = 2 * j;
      size_t index1 = index0 + 1;
      size_t reversedIndex0 = reverseBits16(index0, log2N);
      size_t reversedIndex1 = reverseBits16(index1, log2N);
      __m128 mm_scale = _mm_set1_ps(scale);
      for (size_t i = 0; i < Nq; i++) {
        size_t srcIx0 = (4 * i * N + reversedIndex0);
        __m128 aReal02 = _mm_unpacklo_ps(_mm_load_ss(srcReal + srcIx0 + 0 * N), _mm_load_ss(srcReal + srcIx0 + 2 * N));
        __m128 aReal13 = _mm_unpacklo_ps(_mm_load_ss(srcReal + srcIx0 + 1 * N), _mm_load_ss(srcReal + srcIx0 + 3 * N));
        __m128 aReal = _mm_mul_ps(mm_scale, _mm_unpacklo_ps(aReal02, aReal13));
 
        __m128 aImag02 = _mm_unpacklo_ps(_mm_load_ss(srcImag + srcIx0 + 0 * N), _mm_load_ss(srcImag + srcIx0 + 2 * N));
        __m128 aImag13 = _mm_unpacklo_ps(_mm_load_ss(srcImag + srcIx0 + 1 * N), _mm_load_ss(srcImag + srcIx0 + 3 * N));
        __m128 aImag = _mm_mul_ps(mm_scale, _mm_unpacklo_ps(aImag02, aImag13));
 
        size_t srcIx1 = (4 * i * N + reversedIndex1);
        __m128 bReal02 = _mm_unpacklo_ps(_mm_load_ss(srcReal + srcIx1 + 0 * N), _mm_load_ss(srcReal + srcIx1 + 2 * N));
        __m128 bReal13 = _mm_unpacklo_ps(_mm_load_ss(srcReal + srcIx1 + 1 * N), _mm_load_ss(srcReal + srcIx1 + 3 * N));
        __m128 bReal = _mm_mul_ps(mm_scale, _mm_unpacklo_ps(bReal02, bReal13));
 
        __m128 bImag02 = _mm_unpacklo_ps(_mm_load_ss(srcImag + srcIx1 + 0 * N), _mm_load_ss(srcImag + srcIx1 + 2 * N));
        __m128 bImag13 = _mm_unpacklo_ps(_mm_load_ss(srcImag + srcIx1 + 1 * N), _mm_load_ss(srcImag + srcIx1 + 3 * N));
        __m128 bImag = _mm_mul_ps(mm_scale, _mm_unpacklo_ps(bImag02, bImag13));
 
        __m128 pReal = _mm_add_ps(aReal, bReal);
        __m128 pImag = _mm_add_ps(aImag, bImag);
        size_t dstIx0 = (index0 * N + 4 * i);
        _mm_store_ps(dstReal + dstIx0, pReal);
        _mm_store_ps(dstImag + dstIx0, pImag);
 
        __m128 qReal = _mm_sub_ps(aReal, bReal);
        __m128 qImag = _mm_sub_ps(aImag, bImag);
        size_t dstIx1 = (index1 * N + 4 * i);
        _mm_store_ps(dstReal + dstIx1, qReal);
        _mm_store_ps(dstImag + dstIx1, qImag);
      }
    }
 
    // remaining passes
    for (size_t l = 1; l < log2N; l++) {
      int blockSize = 2 << l;
      for (size_t j = 0; j < N / 2; j++) {
        size_t blockNumber = j >> l;
        size_t blockIndex = j & ((1 << l) - 1);
        size_t index0 = (blockSize * blockNumber + blockIndex);
        size_t index1 = (index0 + blockSize / 2);
        float* real0 = dstReal + index0 * N;
        float* imag0 = dstImag + index0 * N;
        float* real1 = dstReal + index1 * N;
        float* imag1 = dstImag + index1 * N;
 
        int twiddleIndex = int(blockIndex << (log2N - l - 1));  // of N
        const float twiddleArg = (factor * twiddleIndex) / N;
        __m128 wReal = _mm_set1_ps(cos(twiddleArg));
        __m128 wImag = _mm_set1_ps(sin(twiddleArg));
        for (size_t i = 0; i < Nq; i++) {
          __m128 aReal = _mm_load_ps(real0);
          __m128 bReal = _mm_load_ps(real1);
 
          __m128 aImag = _mm_load_ps(imag0);
          __m128 bImag = _mm_load_ps(imag1);
 
          __m128 cReal = _mm_sub_ps(_mm_mul_ps(wReal, bReal), _mm_mul_ps(wImag, bImag));
          __m128 cImag = _mm_add_ps(_mm_mul_ps(wReal, bImag), _mm_mul_ps(wImag, bReal));
 
          _mm_store_ps(real0, _mm_add_ps(aReal, cReal)); real0 += 4;
          _mm_store_ps(imag0, _mm_add_ps(aImag, cImag)); imag0 += 4;
 
          _mm_store_ps(real1, _mm_sub_ps(aReal, cReal)); real1 += 4;
          _mm_store_ps(imag1, _mm_sub_ps(aImag, cImag)); imag1 += 4;
        }
      }
    }
  }
#endif
 
  void inverseRadix2MajorIndexTranspose(float* dstReal,
    float* dstImag,
    const size_t dstStride,
    const float* srcReal,
    const float* srcImag,
    const size_t srcStride,
    const float factor,
    const float scale,
    const size_t log2N)
  {
    size_t N = (size_t)1 << log2N;
 
    // first pass, also handle transpose and shuffle 
    for (size_t j = 0; j < N / 2; j++) {
      size_t index0 = 2 * j;
      size_t index1 = index0 + 1;
      size_t reversedIndex0 = reverseBits16(index0, log2N);
      size_t reversedIndex1 = reverseBits16(index1, log2N);
      for (size_t i = 0; i < N; i++) {
        size_t srcIx0 = srcStride * (i * N + reversedIndex0);
        size_t srcIx1 = srcStride * (i * N + reversedIndex1);
        size_t dstIx0 = dstStride * (index0 * N + i);
        size_t dstIx1 = dstStride * (index1 * N + i);
        complex<float> a = scale * complex<float>(srcReal[srcIx0], srcImag[srcIx0]);
        complex<float> b = scale * complex<float>(srcReal[srcIx1], srcImag[srcIx1]);
        complex<float> p = a + b;
        complex<float> q = a - b;
        dstReal[dstIx0] = p.real();  dstImag[dstIx0] = p.imag();
        dstReal[dstIx1] = q.real();  dstImag[dstIx1] = q.imag();
      }
    }
 
    // remaining passes
    for (size_t l = 1; l < log2N; l++) {
      int blockSize = 2 << l;
      for (size_t j = 0; j < N / 2; j++) {
        size_t blockNumber = j >> l;
        size_t blockIndex = j & ((1 << l) - 1);
        size_t index0 = (blockSize * blockNumber + blockIndex);
        size_t index1 = (index0 + blockSize / 2);
        int twiddleIndex = int(blockIndex << (log2N - l - 1));  // of N
        const float twiddleArg = (factor * twiddleIndex) / N;
        complex<float> w(cos(twiddleArg), sin(twiddleArg));
        for (size_t i = 0; i < N; i++) {
          size_t srcDstIx0 = dstStride * (index0 * N + i);
          size_t srcDstIx1 = dstStride * (index1 * N + i);
          complex<float> a = complex<float>(dstReal[srcDstIx0], dstImag[srcDstIx0]);
          complex<float> b = complex<float>(dstReal[srcDstIx1], dstImag[srcDstIx1]);
          complex<float> c = w * b;
          complex<float> p = a + c;
          complex<float> q = a - c;
          dstReal[srcDstIx0] = p.real();  dstImag[srcDstIx0] = p.imag();
          dstReal[srcDstIx1] = q.real();  dstImag[srcDstIx1] = q.imag();
        }
      }
    }
  }
 
  void fastGenericFourierTransform2D(Cogs::Core::Context* /*context*/,
    uint8_t* scratch,
    uint8_t* dstReal,
    uint8_t* dstImag,
    const size_t dstStride,
    const uint8_t* srcReal,
    const uint8_t* srcImag,
    const size_t srcStride,
    const float factor,
    const float scale,
    const size_t log2N)
  {
    // We assume that we can cast these pointers to float pointers (i.e., multiple of four bytes)
    assert((reinterpret_cast<size_t>(dstReal) & 0x3) == 0);
    assert((reinterpret_cast<size_t>(dstImag) & 0x3) == 0);
    assert((reinterpret_cast<size_t>(srcReal) & 0x3) == 0);
    assert((reinterpret_cast<size_t>(srcImag) & 0x3) == 0);
    assert((dstStride & 0x3) == 0);
    assert((srcStride & 0x3) == 0);
 
    float* scratchf = reinterpret_cast<float*>(32 * ((reinterpret_cast<size_t>(scratch) + 31) / 32));
    float* dstRealf = reinterpret_cast<float*>(dstReal);
    float* dstImagf = reinterpret_cast<float*>(dstImag);
    const float* srcRealf = reinterpret_cast<const float*>(srcReal);
    const float* srcImagf = reinterpret_cast<const float*>(srcImag);
    size_t dstStridef = dstStride / sizeof(float);
    size_t srcStridef = srcStride / sizeof(float);
 
#if !defined(__EMSCRIPTEN__) && !defined(__APPLE__)
    size_t N = size_t(1) << log2N;
    if ((N >= 4) &&
      ((N & 3) == 0) &&
      (srcStridef == 1) &&
      (dstStridef == 1) &&
      ((reinterpret_cast<size_t>(dstReal) & 0xf) == 0) &&
      ((reinterpret_cast<size_t>(dstImag) & 0xf) == 0) &&
      ((reinterpret_cast<size_t>(srcReal) & 0xf) == 0) &&
      ((reinterpret_cast<size_t>(srcImag) & 0xf) == 0))
    {
      inverseRadix2MajorIndexTransposeAlignedSoALog4SSE1(scratchf, scratchf + N * N,
        srcRealf, srcImagf,
        factor, scale, log2N);
 
      inverseRadix2MajorIndexTransposeAlignedSoALog4SSE1(dstRealf, dstImagf,
        scratchf, scratchf + N * N,
        factor, 1.f, log2N);
      return;
    }
#endif
 
    inverseRadix2MajorIndexTranspose(scratchf, scratchf + 1, 2,
      srcRealf, srcImagf, srcStridef,
      factor, scale, log2N);
 
    inverseRadix2MajorIndexTranspose(dstRealf, dstImagf, dstStridef,
      scratchf, scratchf + 1, 2,
      factor, 1.f, log2N);
    /*
    inverseRadix2MajorIndexTranspose(scratchf, scratchf + N*N, 1,
                                     srcRealf, srcImagf, srcStridef,
                                     factor, log2N);
 
    inverseRadix2MajorIndexTranspose(dstRealf, dstImagf, dstStridef,
                                     scratchf, scratchf + N*N, 1,
                                     1.f, log2N);
                                     */
  }
 
  inline size_t fastGenericFourierTransform2DScratchsize(size_t log2N) { return sizeof(float) * (size_t(2) << (log2N << 1)) + 64; }
 
  inline void fastInverseFourierTransform2D(Cogs::Core::Context* context,
    std::vector<uint8_t>& scratch,
    Cogs::Core::ComplexArray& dst,
    const Cogs::Core::ComplexArray& src,
    const float scale,
    const size_t log2N)
  {
    CpuInstrumentationScope(SCOPE_GEOMETRY, "OceanSystem::fastInverseFourierTransform2D");
 
    scratch.resize(fastGenericFourierTransform2DScratchsize(log2N));
    fastGenericFourierTransform2D(context,
      reinterpret_cast<uint8_t*>(scratch.data()),
      reinterpret_cast<uint8_t*>(dst.real()),
      reinterpret_cast<uint8_t*>(dst.imag()),
      sizeof(float),
      reinterpret_cast<const uint8_t*>(src.real()),
      reinterpret_cast<const uint8_t*>(src.imag()),
      sizeof(float),
      float(2.0 * 3.14159265358979323846),
      scale,
      log2N);
  }
 
  // Linear ocean wave theory
  // ------------------------
  //
  // We assume that wave amplitudes are small, and the surface elevation z of a
  // wave traveling along x is
  //
  //     z = a sin( k x - omega t),
  //
  // where
  //
  //     omega = 2 pi f = 2 pi / T   is the wave frequency in rad/s,
  //     f                           is the frequency in Hz,
  //     T                           is the wave period in seconds,
  //     k     = 2 pi / L            is the wave number (in rad/meter I guess)
  //     L                           is the wave length in meters.
  //
  // Here, the wave frequency is described in two domains, omega is from the
  // _time domain_  and k is from the _spatial domain_. These two domains are
  // related by the dispersion relation,
  //
  //    omega^2 = g k tanh(k d),
  //
  // where
  //
  //    g                            is the acceleration of gravity, and
  //    d                            is the water depth.
  //
  // In deep water (d > L/4), we can use the following approximation
  //
  //    omega^2 = g k,
  //
  // and in shallow water (d < L/11), we can use the following approximation
  //
  //     omega^2 = g k^2 d.
  //
  // Phase velocity is the speed at which the wave propagates,
  //
  //     c = omega / k = L / T.
  //
  // The deep and shallow water approximations give
  //
  //     c = sqrt(g/k) = g/omega   in deep-water, and
  //     c = sqrt(g d)             in shallow water.
  //
  // Significant wave height H is given by
  //
  //     H = 4 stddev(z),
  //
  // i.e., four times the standard deviation of the surface displacement.
  //
  //
  // Von Gerstner waves
  // ------------------
  //
  // Linear waves describe waves with sinusoidal shapes, which is appropriate in
  // calm weather. However, when wave steepness increases, the crests become
  // sharper and troughs flatter.
  //
  // This can be modeled by trochoids, that is, a fixed point on a circle as the
  // circle spins,
  //
  //     x = a theta - b sin( theta )
  //     y = a       - b cos( theta ),
  //
  // in other words, add spatial displacement in addition to vertical
  // displacement. To classify trochoids, we have prolate (folding, a/b < 1),
  // common (sharp tip, a/b=1), and prolate (smoothed tip, a/b > 1). When our
  // trochoids are prolate, our water surface folds, which is a hint that we
  // may have reached the limits of the validity of our model.
  //
  // Our sea is a composition of waves, and summing over several wave numbers,
  //
  //     X(u,v,t) = u + \sum_j\sum_i i/|ij| a(ij) sin(dot(ij,uv) - omega(ij)t + phi(ij)),
  //     Y(u,v,t) = v + \sum_j\sum_i j/|ij| a(ij) sin(dot(ij,uv) - omega(ij)t + phi(ij)),
  //     Z(u,v,t) =   - \sum_j\sum_i        a(ij) cos(dot(ij,uv) - omega(ij)t + phi(ij)).
  //
  // To get analytical normal vectors,
  //     dX(u,v,t)/du =                1 + \sum_j\sum_i i^2/|ij| a(ij) cos(dot(ij,uv) - omega(ij)t + phi(ij)),
  //     dX(u,v,t)/dv = dY(u,v,t)/du =     \sum_j\sum_i  ij/|ij| a(ij) cos(dot(ij,uv) - omega(ij)t + phi(ij)),
  //     dY(u,v,t)/dv =                1 + \sum_j\sum_i j^2/|ij| a(ij) cos(dot(ij,uv) - omega(ij)t + phi(ij)),
  //     dZ(u,v,t)/du =                    \sum_j\sum_i        i a(ij) sin(dot(ij,uv) - omega(ij)t + phi(ij)).
  //     dZ(u,v,t)/dv =                    \sum_j\sum_i        j a(ij) sin(dot(ij,uv) - omega(ij)t + phi(ij)).
  //
  // Wave spectrum
  // --------------
  //
  //
  // Swell periods
  // -------------
  //
  // Notes about wave periods according to magicseaweed surf forecasting site:
  //
  // 1-4   seconds: Weak and unsurfable. Sea will look lumpy and bumpy, but
  //                you'll struggle to see individual waves.
  // 5-6   seconds: You will see the odd weak rideable wave face if you're very
  //                desperate.
  // 7-9   seconds: Ok for areas that don't get great waves, will in general be
  //                weaker and jumbled up without clear sets.
  // 10-12 seconds: Swells that often be starting to head away from the storms
  //                that created them, and an travel the open ocean for some
  //                distance. Often creates good quality surf.
  // 13-15 seconds: Definitely groundswell, normally created some considerable
  //                distance from the beach by powerful storms. They most often
  //                arrive without the storm that created them. These swells
  //                will have more defined sets and look a lot more lined up
  //                than lower period swells.
  // 16+   seconds: Extremely powerful swells generated by distant storms often
  //                traveling the breadth of the largest oceans.
 
 
  // z(u,v,t) = \sum_i\sum_j a(i,j)cos(<ij,uv> - omega(k)t + phi(k))
  // dz(u,v,t)/du = \sum\sum -ia(i,j)sin(..)
  // dz(u,v,t)/dv = \sum\sum -ja(i,j)sin(..)
 
  using std::complex;
  using std::sqrt;
  using std::min;
  using std::max;
  using std::pow;
  using std::exp;
  using std::tgamma;
  using std::sin;
  using std::cos;
  using std::acos;
  using std::log;
  using std::polar;
 
  using glm::vec2;
  using glm::vec3;
  using glm::vec4;
  using glm::ivec2;
  using glm::length;
  using glm::dot;
  using glm::cross;
 
  const float pi = std::acos(-1.f);
 
  const float twoPi = float(2.0 * glm::pi<float>());
 
  const float g = 9.80665f; // gravity of earth
 
  float waveSpectrumPhillips(const float omega,
    const float alpha = 0.0081f)
  {
    const float omega2 = omega * omega;
    const float omega4 = omega2 * omega2;
    const float omega5 = omega4 * omega;
 
    return (alpha * g * g) / omega5;
  }
 
  float waveSpectrumPiersonMoskowitz(const float omega,
    const float omega_p, // frequency of the spectrum peak
    const float alpha = 0.0081f)
  {
    const float omega_p_over_omega = omega_p / omega;
    const float omega_p_over_omega4 = (omega_p_over_omega * omega_p_over_omega) * (omega_p_over_omega * omega_p_over_omega);
 
    return waveSpectrumPhillips(omega, alpha) * exp(float(-5.0 / 4.0) * omega_p_over_omega4);
  }
 
  float dispersionLonguetHiggins(const float omega,
    const float omega_p,
    const float U10,
    const float windWaveAngle) // in [-pi,pi]
  {
 
    // Sharpness of directional spreading (Mitsuyasu et al., 1975)
    float mu = omega <= omega_p ? 5.f : -2.5f;
    float s_omega = 11.5f * pow(g / (omega_p * U10), 2.5f) * pow(omega / omega_p, mu);
    assert(std::isfinite(s_omega));
 
    // Normalization factor of dispersion relation
 
    // Part of this is the ratio gamma(x+1)/gamma(x+0.5). Evaluating this
    // directly is begging for overflow since gamma rapidly becomes'
    // ridiculously large, e.g. gamma(30)=8.8x10^30. However, since we are
    // interested in the ratio, the following asymptotic series,
    //
    // gamma(J+1/2)/gamma(J) = sqrt(gamma)*(1 - 1/8J + 1/128J^2
    //                                      5/1024J^3 - 21/32768J^4 + ... )
    //
    // is handy, as if we define J = s_omega + 0.5, we can use it directly
    // to evaluate the ration.
    float J = s_omega + 0.5f;
    float gammaRatio = sqrt(J) * (1.f - 1.f / (8.f * J) + 1.f / (128.f * J * J) + 5.f / (1024.f * J * J * J) - 21.f / (32768.f * J * J * J * J));
    float N_s_omega = float(1.0 / (2.0 * sqrt(glm::pi<float>()))) * (gammaRatio);
    assert(std::isfinite(N_s_omega));
 
    //    x!(x / ((1 / 2 (2 x + 1))!) + 1 / (2 (1 / 2 (2 x + 1))!))
 
        // Wind-wave direction angle halved
    //    float theta_half = 0.5f*(waveDirection);
 
        // Dispersion relation (Longuet-Higgins, 1962)
    float D = N_s_omega * pow(max(0.0f, cos(0.5f * windWaveAngle)), 2.f * s_omega);
    assert(std::isfinite(D));
    return D;
  }
 
 
 
  void createDirectionalWaveSpectrum(std::vector<float>& E,
    size_t N,
    float L,
    float /*omega_pz*/,
    float windSpeed,
    float windDirection,
    float /*significantWaveHeight*/,
    float dominantWavePeriod)
  {
    if (windDirection >= pi) {
      windDirection -= twoPi;
    }
    if (windSpeed < 2.f) {
      windSpeed = 2.f;
    }
 
    float dominantAngularVelocity = float(2.0 * glm::pi<float>()) / dominantWavePeriod;
    float sumS = 0.f;
 
    glm::vec2 windDir(glm::cos(windDirection), glm::sin(windDirection));
 
    for (size_t j = 0; j < N; j++) {
      for (size_t i = 0; i < N; i++) {
        if ((i == 0) && (j == 0)) {
          continue;
        }
        // [0,N/2]     -> [0,N/2]
        // [N/2+1,N-1] -> [-N/2+1,-1]
        ivec2 ij(i <= N / 2 ? int(i) : int(i) - int(N),
          j <= N / 2 ? int(j) : int(j) - int(N));
 
        const vec2 K = (twoPi / L) * vec2(ij);  // K:     Wave vector
        float k = length(K);                  // k:     Wave number
        float angularVelocity = sqrt(g * k);    // omega: Angular velocity
        //float waveDirection = std::atan2(K.y, K.x); // theta_w: Wave direction
 
        float cosWindWaveAngle = std::max(-1.f, std::min(1.f, (1.f / k) * dot(K, windDir)));
        float windWaveAngle = std::acos(cosWindWaveAngle);
        assert(std::isfinite(windWaveAngle));
 
        const float S = waveSpectrumPiersonMoskowitz(angularVelocity,
          dominantAngularVelocity,
          0.0081f);
        const float D = dispersionLonguetHiggins(angularVelocity,
          dominantAngularVelocity,
          windSpeed,
          windWaveAngle);
 
        const float chainFactor = (1.f / (2.f * k)) * sqrt(g / k);
 
        sumS += chainFactor * S;
        E[N * j + i] = chainFactor * S * D;
      }
    }
    E[0] = 0.f;
 
    // Adjust spectrum to match significant wave height
    float w = 1.f / sumS;// (significantWaveHeight*significantWaveHeight) / (sumS);
    for (size_t i = 0; i < N * N; i++) {
      E[i] = w * E[i];
    }
 
  }
 
  void createRandomizedWaveSpectrumInstance(Cogs::Core::ComplexArray& H,
    const std::vector<float>& E,
    const unsigned int seed,
    const size_t N)
  {
    std::minstd_rand eng(seed);
    std::uniform_int_distribution<int> dist(0, RAND_MAX);
 
    for (size_t j = 0; j < N; j++) {
      for (size_t i = 0; i < N; i++) {
#if 0
        float eta0 = 0.f;
        float eta1 = 0.f;
        for (int q = 0; q < 12; q++) {
          eta0 += float(dist(eng)) / float(RAND_MAX);
          eta1 += float(dist(eng)) / float(RAND_MAX);
        }
        eta0 -= 12.f;
        eta1 -= 12.f;
        H.store(N * j + i, sqrt(E[N * j + i] / 2.f) * std::complex<float>(eta0, eta1));
#else
 
        // Check the sanity of this trying to comply with normal distribution
 
        float U0 = float(dist(eng) + 1) / (float(RAND_MAX) + 2);
        float U1 = (2.f * float(glm::pi<float>())) * (float(dist(eng)) / (float(RAND_MAX) + 1));
 
        // Box-Muller transform to create normal distributed numbers, mean 0, var 1.
        complex<float> eta = sqrt(-2.f * log(U0)) * complex<float>(cos(U1), sin(U1));
 
        H.store(N * j + i, sqrt(E[N * j + i] / 2.f) * eta);
#endif
      }
    }
  }
 
  void fastSinCos(float& sin, float& cos, const float x)
  {
    float r = (2.0f / glm::pi<float>()) * x;
 
    int k = int(r);
    float x_ = r - k;
    float y_ = 1.f - x_;
 
    float vx = x_ * x_ * (x_ * (2.f - 33.f / 20.f) + (33.f / 20.f - 3.f)) + 1.f;
    float vy = y_ * y_ * (y_ * (2.f - 33.f / 20.f) + (33.f / 20.f - 3.f)) + 1.f;
 
    float c = ((k & 1) == 0) ? vx : vy;
    float s = ((k & 1) == 0) ? vy : vx;
 
    cos = (((k + 1) & 2) == 0) ? c : -c;
    sin = ((k & 2) == 0) ? s : -s;
  }
 
  complex<float> fastPolar(float rho, float theta)
  {
    float cos, sin;
 
    fastSinCos(sin, cos, theta);
 
    return complex<float>(rho * cos, rho * sin);
  }
 
  void disperseWaves(Cogs::Core::ComplexArray& dZdu,
    Cogs::Core::ComplexArray& dZdv,
    Cogs::Core::ComplexArray& Ax,
    Cogs::Core::ComplexArray& Ay,
    Cogs::Core::ComplexArray& Az,
    const Cogs::Core::ComplexArray& a,
    const float t,
    const float L,
    const size_t N,
    const size_t start,
    const size_t end)
  {
    for (size_t j = start; j < end; j++) {
      for (size_t i = 0; i < N; i++) {
        const vec2 K = (twoPi / L) * vec2(i <= N / 2 ? int(i) : int(i) - int(N),
          j <= N / 2 ? int(j) : int(j) - int(N)); // K:     Wave vector
        float k = length(K);                                                    // k:     Wave number
        vec2 kn = (1.f / k) * K;
        float omega = sqrt(g * k);                                                // omega: Angular velocity
 
        complex<float> A = a.load(N * j + i) * fastPolar(1.f, omega * t);
 
        dZdu.store(N * j + i, K.x * A);
        dZdv.store(N * j + i, K.y * A);
        Ax.store(N * j + i, kn.x * A);
        Ay.store(N * j + i, kn.y * A);
        Az.store(N * j + i, A);
      }
    }
    dZdu.store(0, complex<float>(0.f));
    dZdv.store(0, complex<float>(0.f));
    Ax.store(0, complex<float>(0.f));
    Ay.store(0, complex<float>(0.f));
    Az.store(0, complex<float>(0.f));
  }
 
  void encodeTextureData(vec4* tangents,
    vec4* real,
    vec4* imag,
    const Cogs::Core::ComplexArray& dZdu,
    const Cogs::Core::ComplexArray& dZdv,
    const Cogs::Core::ComplexArray& Dx,
    const Cogs::Core::ComplexArray& Dy,
    const Cogs::Core::ComplexArray& Dz,
    const size_t begin,
    const size_t end)
  {
    for (size_t i = begin; i < end; i++) {
      const size_t ix = i;
 
      glm::vec3 re = glm::vec3(Dx.load(ix).real(), Dy.load(ix).real(), Dz.load(ix).real());
      glm::vec3 im = glm::vec3(Dx.load(ix).imag(), Dy.load(ix).imag(), Dz.load(ix).imag());
      glm::vec4 ta = glm::vec4(dZdu.load(ix).real(), dZdu.load(ix).imag(), dZdv.load(ix).real(), dZdv.load(ix).imag());
 
      real[ix] = vec4(re, 0.f);
      imag[ix] = vec4(im, 0.f);
      tangents[ix] = ta;
    }
  }
 
  void encodeTextureDataTangents(glm::u8vec4* tangents,
    glm::u8vec4* real,
    float& magnitude0Out,
    float& magnitude1Out,
    const Cogs::Core::ComplexArray& dZdu,
    const Cogs::Core::ComplexArray& dZdv,
    const Cogs::Core::ComplexArray& Dx,
    const Cogs::Core::ComplexArray& Dy,
    const Cogs::Core::ComplexArray& Dz,
    const float magnitude0In,
    const float magnitude1In,
    const size_t begin,
    const size_t end)
  {
    const float scale0 = 0.5f * 255.f / magnitude0In;
    const float scale1 = 0.5f * 255.f / magnitude1In;
 
    float magnitude0 = 0.f;
    float magnitude1 = 0.f;
    for (size_t i = begin; i < end; i++) {
      const size_t ix = i;
 
      glm::vec3 re = glm::vec3(Dx.load(ix).real(), Dy.load(ix).real(), Dz.load(ix).real());
      glm::vec3 im = glm::vec3(Dx.load(ix).imag(), Dy.load(ix).imag(), Dz.load(ix).imag());
      glm::vec4 ta = glm::vec4(dZdu.load(ix).real(), dZdu.load(ix).imag(), dZdv.load(ix).real(), dZdv.load(ix).imag());
 
      magnitude0 = std::max(std::max(std::max(std::abs(re.x), std::abs(re.y)),
        std::max(std::abs(re.z), std::abs(im.x))),
        std::max(std::max(std::abs(im.y), std::abs(im.z)),
          magnitude0));
 
      magnitude1 = std::max(std::max(std::max(std::abs(ta.x), std::abs(ta.y)),
        std::max(std::abs(ta.z), std::abs(ta.w))),
        magnitude1);
 
      real[ix] = vec4(glm::clamp(scale0 * (re + glm::vec3(magnitude0In)), glm::vec3(0.f), glm::vec3(255.f)), 0.f);
      tangents[ix] = glm::clamp(scale1 * (ta + vec4(magnitude1In)), glm::vec4(0.f), glm::vec4(255.f));
    }
    magnitude0Out = magnitude0;
    magnitude1Out = magnitude1;
  }
 
 
  void encodeTextureDataQuux(vec4* pos,
    vec4* nrm,
    const Cogs::Core::ComplexArray& dZdu,
    const Cogs::Core::ComplexArray& dZdv,
    const Cogs::Core::ComplexArray& Dx,
    const Cogs::Core::ComplexArray& Dy,
    const Cogs::Core::ComplexArray& Dz,
    const float L,
    const float significantWaveHeight,
    const size_t begin,
    const size_t end,
    const size_t m)
  {
    const size_t N = (size_t)1 << m;
    const auto M = N - 1; // assume N is a power of two.
    const auto s = 1.f / L;
    const auto l = L / N;
    const auto h = significantWaveHeight;
    //const auto worldSpaceMeter = 1.f;
 
    // Displacement function is
    //
    //          +-                 -+
    //          |  u + Ix(s*u, s*v) |
    // P(u,v) = |  v + Iy(s*u, s*v) |
    //          | -h * Rz(s*u, s*v) |
    //          +-                 -+
    //
    // where s = 1.f / fftTileExtent and h is significant wave height,
    // and the two partial derivatives is thus
    //
    //              +-                         -+
    //              |  1 + s * dIx(s*u, s*v)/du |
    // dP(u,v)/du = |      s * dIy(s*u, s*v)/du |
    //              | -s * h * dRz(s*u, s*v)/du |
    //              +-                         -+
    //
    //              +-                         -+
    //              |      s * dIx(s*u, s*v)/dv |
    // dP(u,v)/dv = |  1 + s * dIy(s*u, s*v)/dv |
    //              | -s * h * dRz(s*u, s*v)/dv |
    //              +-                         -+
    //
    // We have dRz/du and dRz/dv explicitly. dIx/du|dv and dIy/du|dv is
    // approximated using forward differences,
    //
 
    for (size_t j = begin; j < end; j++) {
      size_t jp = (j + 1) & M;
 
      for (size_t i = 0; i < N; i++) {
        size_t ip = (i + 1) & M;
 
        float x = Dx.load((j << m) + i).imag();
        float y = Dy.load((j << m) + i).imag();
        float z = -Dz.load((j << m) + i).real();
 
        float dx_du = l * (Dx.load((j << m) + ip).imag() - x);
        float dy_du = l * (Dy.load((j << m) + ip).imag() - y);
 
        float dx_dv = l * (Dx.load((jp << m) + i).imag() - x);
        float dy_dv = l * (Dy.load((jp << m) + i).imag() - y);
 
        float zdu = dZdu.load((j << m) + i).imag();
        float zdv = dZdv.load((j << m) + i).imag();
 
        glm::vec3 u = glm::vec3(s * dx_du + 1.f, s * dy_du + 0.f, h * zdu);
        glm::vec3 v = glm::vec3(s * dx_dv + 0.f, s * dy_dv + 1.f, h * zdv);
        glm::vec3 n = glm::normalize(glm::cross(u, v));
 
        pos[(j << m) + i] = glm::vec4(x, y, z, 0);
        nrm[(j << m) + i] = glm::vec4(n.x, n.y, n.z, 0);
      }
    }
  }
 
 
  void encodeTextureDataQuux(glm::u8vec4* pos,
    glm::u8vec4* nrm,
    float& magnitude0Out,
    float& magnitude1Out,
    const Cogs::Core::ComplexArray& dZdu,
    const Cogs::Core::ComplexArray& dZdv,
    const Cogs::Core::ComplexArray& Dx,
    const Cogs::Core::ComplexArray& Dy,
    const Cogs::Core::ComplexArray& Dz,
    const float magnitude0In,
    const float magnitude1In,
    const float L,
    const float significantWaveHeight,
    const size_t begin,
    const size_t end,
    const size_t m)
  {
    float magnitude0 = 0.f;
    float magnitude1 = 0.f;
    const size_t N = (size_t)1 << m;
    const auto M = N - 1; // assume N is a power of two.
    const auto s = 1.f / L;
    const auto l = L / N;
    const auto h = significantWaveHeight;
 
    const float scale0 = 0.5f * 255.f / magnitude0In;
    const float scale1 = 0.5f * 255.f / magnitude1In;
    for (size_t j = begin; j < end; j++) {
      size_t jp = (j + 1) & M;
 
      for (size_t i = 0; i < N; i++) {
        size_t ip = (i + 1) & M;
 
        glm::vec3 p = glm::vec3(Dx.load((j << m) + i).imag(),
          Dy.load((j << m) + i).imag(),
          -Dz.load((j << m) + i).real());
 
        float dx_du = l * (Dx.load((j << m) + ip).imag() - p.x);
        float dy_du = l * (Dy.load((j << m) + ip).imag() - p.y);
 
        float dx_dv = l * (Dx.load((jp << m) + i).imag() - p.x);
        float dy_dv = l * (Dy.load((jp << m) + i).imag() - p.y);
 
        float zdu = dZdu.load((j << m) + i).imag();
        float zdv = dZdv.load((j << m) + i).imag();
 
        glm::vec3 u = glm::vec3(s * dx_du + 1.f, s * dy_du + 0.f, h * zdu);
        glm::vec3 v = glm::vec3(s * dx_dv + 0.f, s * dy_dv + 1.f, h * zdv);
        glm::vec3 n = glm::normalize(glm::cross(u, v));
 
        magnitude0 = std::max(std::max(std::abs(p.x), std::abs(p.y)),
          std::max(std::abs(p.z), magnitude0));
 
        // Note: z-component of normal is usually quite close to 1,
        // so we encode difference from 1 instead to make range match n.x and n.y.
        magnitude1 = std::max(std::max(std::abs(n.x), std::abs(n.y)),
          std::max(std::abs(n.z - 1.f), magnitude1));
 
        pos[(j << m) + i] = glm::vec4(glm::clamp(scale0 * (p + glm::vec3(magnitude0In)), glm::vec3(0.f), glm::vec3(255.f)), 0);
        nrm[(j << m) + i] = glm::vec4(glm::clamp(scale1 * glm::vec3(n.x + magnitude1In,
          n.y + magnitude1In,
          n.z - 1.f), glm::vec3(0.f), glm::vec3(255.f)), 0);
      }
    }
    magnitude0Out = magnitude0;
    magnitude1Out = magnitude1;
 
  }
}
 
using glm::normalize;
using glm::clamp;
using glm::dvec2;
 
void Cogs::Core::BasicOceanSystem::initialize(Context * context)
{
  ComponentSystem::initialize(context);
 
  DisplacementTexH = context->textureManager->create();
  NormalTexH = context->textureManager->create();
  TangentsTexH = context->textureManager->create();
 
  DisplacementTexH->setName("BasicOcean.Displacement");
  NormalTexH->setName("BasicOcean.Normals");
  TangentsTexH->setName("BasicOcean.Tangents");
  
  oceanTaskGroup = context->taskManager->createGroup(TaskManager::GlobalQueue);
 
  setupMaterial();
  setupWaveSpectrum();
  
  Variable* animateVar = context->variables->get(animateKey);
  if (animateVar->isEmpty()) {
    animateVar->setBool(true);
  }
  Variable* eightbitVar = context->variables->get(eightbitKey);
  if (eightbitVar->isEmpty()) {
    eightbitVar->setBool(true);
  }
}
 
void Cogs::Core::BasicOceanSystem::cleanup(Context * context)
{
  if (oceanTaskGroup.isValid()) {
    context->taskManager->destroy(oceanTaskGroup);
  }
  DisplacementTexH = TextureHandle();
  NormalTexH = TextureHandle();
  TangentsTexH = TextureHandle();
  ReflectionTextureH = TextureHandle();
 
  oceanMaterial = MaterialHandle();
  oceanMaterial2 = MaterialHandle();
}
 
void Cogs::Core::BasicOceanSystem::setupMaterial()
{
  auto adaptiveMat = context->materialManager->loadMaterial("AdaptiveGridMaterial.material");
  oceanMaterial = context->materialManager->loadMaterial("BasicOceanMaterial.material");
  oceanMaterial2 = context->materialManager->loadMaterial("BasicSkyMaterial.material");
 
 
  //NOTE: This should all be moved to after loading has been performed.
  context->materialManager->processLoading();
 
  oceanMaterial->options.cullMode = CullMode::None;
 
  auto m = oceanMaterial.resolve();
 
  DisplacementKey = m->getTextureKey("Displacement");
  NormalKey = m->getTextureKey("Difference");
 
  TangentsKey = m->getTextureKey("Tangents");
 
  m->setTextureProperty(DisplacementKey, DisplacementTexH);
  m->setTextureProperty(NormalKey, NormalTexH);
  m->setTextureProperty(TangentsKey, TangentsTexH);
 
  ReflectionTextureH = context->textureManager->create();
  ReflectionTextureH->setName("BasicOcean.Reflection");
 
  PlanarReflectionKey = m->getTextureKey("PlanarReflection");
  m->setTextureProperty(PlanarReflectionKey, ReflectionTextureH);
  m->setTextureAddressMode(PlanarReflectionKey, SamplerState::Clamp);
 
  cameraYAxisKey = m->getVec4Key("cameraYAxis");
  waterColorKey = m->getVec4Key("waterColor");
  waveDirectionKey = m->getVec2Key("waveDirection");
  camPlaneDirKey = m->getVec2Key("camPlaneDir");
  significantWaveHeightKey = m->getFloatKey("significantWaveHeight");
  fftTileScaleKey = m->getFloatKey("fftTileScale");
  camAzimuthKey = m->getFloatKey("camAzimuth");
  seaLevelKey = m->getFloatKey("seaLevel");
  reflectionBrightnessKey = m->getFloatKey("reflectionBrightness");
  phaseShiftNoiseFrequencyKey = m->getFloatKey("phaseShiftNoiseFrequency");
  phaseShiftNoisePeriodKey = m->getFloatKey("phaseShiftNoisePeriod");
}
 
void Cogs::Core::BasicOceanSystem::updateTextureResolution(BasicOceanComponent* oceanComp)
{
  const auto textureResolution = static_cast<uint32_t>(std::max(1,context->variables->get("ocean.reflectionTextureResolution", 1024)));
  auto reflectionTex = context->textureManager->get(ReflectionTextureH);
 
  TextureFormat format = TextureFormat::R8G8B8A8_UNORM_SRGB;
  if (oceanComp != nullptr) {
    format = parseTextureFormat(oceanComp->reflectionTexFormat, format);
  }
 
  if ((reflectionTex->description.width != textureResolution)
      || (reflectionTex->description.format != format))
  {
    reflectionTex->description.width = textureResolution;
    reflectionTex->description.height = textureResolution;
    reflectionTex->description.format = format;
    reflectionTex->description.flags = TextureFlags::RenderTarget;
    reflectionTex->setChanged();
  }
}
 
void Cogs::Core::BasicOceanSystem::setupWaveSpectrum()
{
  size_t N = size_t(1) << fftTileResolutionLog2;
  size_t N_times_N = N*N;
 
  constexpr float oneOvertwoPi = 1.f / (2.0f * std::numbers::pi_v<float>);
  const float dominantWaveLength = oneOvertwoPi*(g*dominantWavePeriod*dominantWavePeriod);
  
  fftTileExtent = dominantWaveLength;
 
  P.resize(N_times_N);
  frqH0.resize(N_times_N);
  frqDx.resize(N_times_N);
  frqDy.resize(N_times_N);
  frqDz.resize(N_times_N);
  frqdDzdu.resize(N_times_N);
  frqdDzdv.resize(N_times_N);
 
  spcDx.resize(N_times_N);
  spcDy.resize(N_times_N);
  spcDz.resize(N_times_N);
  spcdDzdu.resize(N_times_N);
  spcdDzdv.resize(N_times_N);
 
  fftScratch.resize(fastGenericFourierTransform2DScratchsize(N_times_N));
 
  const float omega_p = (0.855f*g) / windSpeed;
 
  createDirectionalWaveSpectrum(P, N, fftTileExtent, omega_p, windSpeed, 0.f /*windDirection*/, 1.f, dominantWavePeriod);
  createRandomizedWaveSpectrumInstance(frqH0, P, 42, N);
}
 
void Cogs::Core::BasicOceanSystem::updateTextures(const float magnitudeIn0, const float magnitudeIn1)  // Max magnitude of data kind 0 found while encoding
{
  CpuInstrumentationScope(SCOPE_SYSTEMS, "OceanSystem::updateTextures");
 
  const uint16_t N = uint16_t(1) << fftTileResolutionLog2;
 
  auto realTex = context->textureManager->get(DisplacementTexH);
  auto imagTex = context->textureManager->get(NormalTexH);
  auto tangentsTex = context->textureManager->get(TangentsTexH);
 
  size_t M = 0;
  switch (waves)
  {
  case BasicOceanWaves::Default: M = static_cast<size_t>(N) * N; break;
  case BasicOceanWaves::Quux: M = static_cast<size_t>(N); break;
  default: assert(false && "Illegal wave type"); break;
  }
 
  size_t split = std::min(size_t(8), 1 + (context->engine->workParallel() ? context->taskManager->getQueueConcurrency(TaskManager::GlobalQueue) : size_t(0)));
  size_t incr = (M + split - 1) / split;
 
  TaskId gr = 1 < split ? context->taskManager->createGroup(TaskManager::GlobalQueue) : NoTask;
 
  std::vector<float> magnitudes_scratch0(split, 0.f);
  std::vector<float> magnitudes_scratch1(split, 0.f);
  if (eightBit) {
    MappedTexture<glm::u8vec4> real = realTex->map<glm::u8vec4>(N, N, false);
    switch (waves)
    {
    case BasicOceanWaves::Default:
    {
      MappedTexture<glm::u8vec4> tangents = tangentsTex->map<glm::u8vec4>(N, N, true);
      for (size_t i = 0; i < split; ++i) {
        TaskFunction f = [this, i, magnitudeIn0, magnitudeIn1, incr, M, &tangents, &real, &magnitudes_scratch0, &magnitudes_scratch1]() {
          CpuInstrumentationScope(SCOPE_SYSTEMS, "OceanSystem::encodeTextureData");
          encodeTextureDataTangents(tangents.data(), real.data(), magnitudes_scratch0[i], magnitudes_scratch1[i],
                            spcdDzdu, spcdDzdv, spcDx, spcDy, spcDz,
                            magnitudeIn0, magnitudeIn1,
                            i * incr, std::min(M, (i + 1) * incr));
        };
        if (1 < split) { context->taskManager->enqueueChild(gr, f); } else { f(); }
      }
 
      if (gr.isValid()) {
        context->taskManager->destroy(gr);
        gr = NoTask;
      }
      break;
    }
    case BasicOceanWaves::Quux:
      MappedTexture<glm::u8vec4> imag = imagTex->map<glm::u8vec4>(N, N, true);
      for (size_t i = 0; i < split; ++i) {
        TaskFunction f = [this, i, magnitudeIn0, magnitudeIn1, incr, M, &real, &imag, &magnitudes_scratch0, &magnitudes_scratch1]() {
          CpuInstrumentationScope(SCOPE_SYSTEMS, "OceanSystem::encodeTextureData");
          encodeTextureDataQuux(real.data(), imag.data(), magnitudes_scratch0[i], magnitudes_scratch1[i],
                                spcdDzdu, spcdDzdv, spcDx, spcDy, spcDz,
                                magnitudeIn0, magnitudeIn1,
                                fftTileExtent, significantWaveHeight,
                                i * incr, std::min(M, (i + 1) * incr),
                                fftTileResolutionLog2);
        };
        if (1 < split) { context->taskManager->enqueueChild(gr, f); } else { f(); }
      }
 
      if (gr.isValid()) {
        context->taskManager->destroy(gr);
        gr = NoTask;
      }
      break;
    }
  }
  else {
    MappedTexture<glm::vec4> real = realTex->map<glm::vec4>(N, N, true);
    MappedTexture<glm::vec4> imag = imagTex->map<glm::vec4>(N, N, true);
    switch (waves)
    {
    case BasicOceanWaves::Default:
    {
      MappedTexture<glm::vec4> tangents = tangentsTex->map<glm::vec4>(N, N, true);
      for (size_t i = 0; i < split; ++i) {
        TaskFunction f = [this, i, incr, M, &tangents, &real, &imag]() {
          CpuInstrumentationScope(SCOPE_SYSTEMS, "OceanSystem::encodeTextureData");
          encodeTextureData(tangents.data(), real.data(), imag.data(), 
                            spcdDzdu, spcdDzdv, spcDx, spcDy, spcDz,
                            i * incr, std::min(M, (i + 1) * incr));
        };
        if (1 < split) { context->taskManager->enqueueChild(gr, f); }
        else { f(); }
      }
 
      if (gr.isValid()) {
        context->taskManager->destroy(gr);
        gr = NoTask;
      }
      break;
    }
    case BasicOceanWaves::Quux:
      for (size_t i = 0; i < split; ++i) {
        TaskFunction f = [this, i, incr, M, &real, &imag]() {
          CpuInstrumentationScope(SCOPE_SYSTEMS, "OceanSystem::encodeTextureData");
          encodeTextureDataQuux(real.data(), imag.data(),
                                spcdDzdu, spcdDzdv, spcDx, spcDy, spcDz,
                                fftTileExtent, significantWaveHeight,
                                i * incr, std::min(M, (i + 1) * incr),
                                fftTileResolutionLog2);
        };
        if (1 < split) { context->taskManager->enqueueChild(gr, f); }
        else { f(); }
      }
 
      if (gr.isValid()) {
        context->taskManager->destroy(gr);
        gr = NoTask;
      }
      break;
    }
  }
  if (gr.isValid()) context->taskManager->destroy(gr);
 
  if (eightBit) {
    magnitudes_tmp.channel0 = *std::max_element(magnitudes_scratch0.begin(), magnitudes_scratch0.end());
    magnitudes_tmp.channel1 = *std::max_element(magnitudes_scratch1.begin(), magnitudes_scratch1.end());
  }
}
 
void Cogs::Core::BasicOceanSystem::updateTileMaterialInstances(const BasicOceanData & oceanData,
                                                               AdaptivePlanarGridComponent * /*gridComp*/,
                                                               AdaptivePlanarGridData & gridData,
                                                               const glm::mat4 & viewToWorld,
                                                               const glm::vec2 & viewPortSize)
{
  const glm::vec2 camPlaneDir = glm::normalize(glm::vec2(viewToWorld[2]));
  float camAzimuth = (0.5f / glm::pi<float>()) * acos(camPlaneDir.x);
  
  if (camPlaneDir.y > 0.f) {
    camAzimuth = -camAzimuth;
  }
 
  camAzimuth = camAzimuth + 0.75f;
 
  const vec4 camYAxis = vec4(normalize(vec3(viewToWorld[1])), 0.5f * viewPortSize.y);
  const vec2 waveDirection(cos(windDirection), sin(windDirection));
  const vec4 rgba(vec3(oceanData.color), clamp(1.f - oceanData.transparency, 0.f, 1.f));
 
  auto m = oceanMaterial.resolve();
  m->setVec4Property(cameraYAxisKey, camYAxis);
  m->setVec4Property(waterColorKey, rgba);
  m->setVec2Property(waveDirectionKey, vec2(waveDirection));
  m->setVec2Property(camPlaneDirKey, camPlaneDir);
  m->setFloatProperty(significantWaveHeightKey, significantWaveHeight);
  m->setFloatProperty(fftTileScaleKey, 1.f / fftTileExtent);
  m->setFloatProperty(camAzimuthKey, camAzimuth);
  m->setFloatProperty(seaLevelKey, oceanData.seaLevel);
  m->setFloatProperty(reflectionBrightnessKey, oceanData.reflectionBrightness);
  m->setFloatProperty(phaseShiftNoiseFrequencyKey, float(phaseShiftNoiseFrequency));
  m->setFloatProperty(phaseShiftNoisePeriodKey, float(phaseShiftNoiseFrequency*tilePeriod));
 
  m->setVec2Property(m->getVec2Key("viewportScale"), 2.f*viewPortSize);
  m->setVec2Property(m->getVec2Key("valueScale"), 2.f * glm::vec2(magnitudes.channel0, magnitudes.channel1));
 
 
  for (auto & tile : gridData.tiles) {
    if (oceanData.transparent) {
      tile.materialInstance->setTransparent();
    }
    else {
      tile.materialInstance->setOpaque();
    }
  }
}
 
void Cogs::Core::BasicOceanSystem::destroyDefaultReflectionCameraIfExists(Context* context, BasicOceanData& oceanData)
{
  if (!oceanData.defaultReflectionCamera) return;
 
  if (TransformComponent* transform = oceanData.defaultReflectionCamera->getComponent<TransformComponent>(); transform && transform->parent) {
    if (Entity* parent = transform->parent.resolve()->getContainer(); parent) {
      context->store->removeChild(parent, oceanData.defaultReflectionCamera.get());
    }
  }
 
  context->store->destroyEntity(oceanData.defaultReflectionCamera->getId());
  oceanData.defaultReflectionCamera.reset();
  LOG_DEBUG(logger, "Destroyed default reflection camera");
}
 
void Cogs::Core::BasicOceanSystem::destroyComponent(ComponentHandle component)
{
  BasicOceanData& oceanData = getData(component.resolveComponent<BasicOceanComponent>());
  destroyDefaultReflectionCameraIfExists(context, oceanData);
  base::destroyComponent(component);
}
 
void Cogs::Core::BasicOceanSystem::update(Context * context)
{
  if (!pool.size()) return;
 
  updateTextureResolution(&(*pool.begin()));
 
  bool encodingChanged = false;
  {
    // Use 8-bit encoding if float is not supported or it is requested
    bool eightBit = ((context->device->getCapabilities()->getDeviceCapabilities().FloatTextures == false) ||
                     (context->variables->get("basic-ocean.8bit", true)));
    encodingChanged = this->eightBit != eightBit;
    this->eightBit = eightBit;
  }
 
  bool visibleInstances = false;
  for (BasicOceanComponent& oceanComp : pool) {
    BasicOceanData& oceanData = getData(&oceanComp);
    EntityPtr reflectionCamera;
    if (oceanComp.reflectionCamera) {
      destroyDefaultReflectionCameraIfExists(context, oceanData);
      reflectionCamera = oceanComp.reflectionCamera;
    }
    else {
      if (!oceanData.defaultReflectionCamera) {
        if (CameraComponent* mainCameraComp = context->cameraSystem->getMainCamera(); mainCameraComp) {
          oceanData.defaultReflectionCamera = context->store->createChildEntity("ReflectionCamera", mainCameraComp->getContainer(), "BasicOcean Default Reflection Camera");
          LOG_DEBUG(logger, "Created default reflection camera");
        }
        else {
          return;
        }
      }
      reflectionCamera = oceanData.defaultReflectionCamera;
    }
 
    auto * gridComp = oceanComp.getComponent<AdaptivePlanarGridComponent>();
    auto * sceneComp = oceanComp.getComponent<SceneComponent>();
 
    visibleInstances = visibleInstances || sceneComp->visible;
 
    if (!oceanData.initialized) {
      context->adaptivePlanarGridSystem->registerMaterial(gridComp, oceanMaterial, initMaterialInstanceCallback, &oceanData);
      gridComp->layer = RenderLayers::Ocean;
      gridComp->setChanged();
      oceanData.initialized = true;
    }
 
    if (oceanData.reflectionCamera.lock() != reflectionCamera) {
      oceanData.reflectionCamera = reflectionCamera;
      if (ReflectionComponent* reflectionComponent = reflectionCamera->getComponent<ReflectionComponent>(); reflectionComponent) {
        reflectionComponent->texture = ReflectionTextureH;
        sceneComp->setChanged();  // Trigger EnableRender code below.
        LOG_DEBUG(logger, "Detected new reflection camera, setting up reflection component");
      }
      else {
        LOG_DEBUG(logger, "Reflection camera lacks ReflectionComponent");
      }
    }
 
    auto reflectionTex = context->textureManager->get(ReflectionTextureH);
 
    CameraComponent* refCamComp = reflectionCamera->getComponent<CameraComponent>();
    refCamComp->viewportSize = glm::vec2(reflectionTex->description.width, reflectionTex->description.height);
 
    if (sceneComp->hasChanged() || oceanComp.hasChanged()) {
      if (sceneComp->visible) {
        switch (oceanComp.reflection) {
        case BasicOceanReflection::Planar:
          refCamComp->flags |= CameraFlags::EnableRender;
          break;
        case BasicOceanReflection::EnvSkyBox:
        case BasicOceanReflection::EnvRadiance:
          refCamComp->flags &= ~CameraFlags::EnableRender;
          break;
        default:
          assert(false);
        }
      }
      else {
        refCamComp->flags &= ~CameraFlags::EnableRender;
      }
    }
 
    if (oceanComp.hasChanged() || encodingChanged) {
      specsChanged = true;
 
      const char * waveVariant = nullptr;
      switch (oceanComp.waves)
      {
      case BasicOceanWaves::Default: waveVariant = "Basic"; break;
      case BasicOceanWaves::Quux: waveVariant = "Quux"; break;
      default: assert(false);
      }
      
      const char* beliefSystem = nullptr;
      switch (oceanComp.beliefSystem) {
      case BasicOceanBeliefSystem::FlatEarth: beliefSystem = "FlatEarth"; break;
      case BasicOceanBeliefSystem::CurvedEarth: beliefSystem = "CurvedEarth"; break;
      default: assert(false);
      }
 
      const char * reflectionVariant = nullptr;
      switch (oceanComp.reflection) {
      case BasicOceanReflection::Planar: reflectionVariant = "Planar"; break;
      case BasicOceanReflection::EnvSkyBox: reflectionVariant = "EnvSkyBox"; break;
      case BasicOceanReflection::EnvRadiance: reflectionVariant = "EnvRadiance"; break;
      default: assert(false);
      }
 
      const char* lightModelVariant = nullptr;
      lightModelVariant = oceanComp.pbr ? "PBR" : "Phong";
 
      const char* encodingVariant = eightBit ? "Scaled" : "AsIs";
 
      if ((oceanData.encoding != encodingVariant)  || 
          (oceanData.waveVariant != waveVariant) ||
          (oceanData.reflectionVariant != reflectionVariant) ||
          (oceanData.lightModelVariant != lightModelVariant) ||
          (oceanData.beliefSystem != beliefSystem))
      {
        oceanData.encoding = encodingVariant;
        oceanData.waveVariant = waveVariant;
        oceanData.beliefSystem = beliefSystem;
        oceanData.reflectionVariant = reflectionVariant;
        oceanData.lightModelVariant = lightModelVariant;
 
        for (auto & instance : context->adaptivePlanarGridSystem->getData(gridComp).materialPool) {
          instance->setVariant("Encoding", encodingVariant);
          instance->setVariant("Waves", waveVariant);
          instance->setVariant("BeliefSystem", beliefSystem);
          instance->setVariant("Reflection", reflectionVariant);
          instance->setVariant("LightModel", lightModelVariant);
        }
      }
 
 
      oceanData.transparency = oceanComp.transparency;
      oceanData.color = oceanComp.color;
      oceanData.transparent = oceanComp.transparency > 0.0f;
      oceanData.seaLevel = oceanComp.seaLevel;
      oceanData.reflectionBrightness = oceanComp.reflectionBrightness;
 
      fftTileResolutionLog2 = std::max(1, oceanComp.fftTileResolutionLog2);
      significantWaveHeight = oceanComp.significantWaveHeight;
      dominantWavePeriod = oceanComp.dominantWavePeriod;
      windSpeed = oceanComp.windSpeed;
      windDirection = oceanComp.windDirection;
      waves = oceanComp.waves;
 
      auto adjustedDisplacement = std::max(significantWaveHeight, 0.01f);
      gridComp->displaceMin = glm::vec3(-0.75f * adjustedDisplacement) + glm::vec3(0.f, 0.f, oceanComp.seaLevel);
      gridComp->displaceMax = glm::vec3(0.75f * adjustedDisplacement)  + glm::vec3(0.f, 0.f, oceanComp.seaLevel);
      gridComp->setChanged();
    }
  }
 
  if (specsChanged) {
    setupWaveSpectrum();
 
    for (auto & oceanComp : pool) {
      auto gridComp = oceanComp.getComponent<AdaptivePlanarGridComponent>();
 
      auto c = glm::cos(windDirection);
      auto s = glm::sin(windDirection);
 
      context->adaptivePlanarGridSystem->setTexCoordTransform(gridComp, glm::mat2(c, s, -s, c), glm::vec2(fftTileExtent*tilePeriod));
    }
  }
 
  if (!visibleInstances) return;  // No visible instances, no point in updating shared data
 
  // We want to get one extra frame when we switch from animate == true to false
  // to set up the waves at a fixed time.
  bool doAnimate = context->variables->get("basic-ocean.animate", true);
  if (!doAnimate && !specsChanged) {
    if (animate) {
      specsChanged = true;
      animate = false;
    }
    else {
      return;
    }
  }
  else {
    animate = true;
  }
 
  context->engine->setDirty(); // Trigger new frame as ocean is animated
 
 
  const int N = 1 << fftTileResolutionLog2;
 
  const bool workParallel = context->engine->workParallel();
 
  extraStep = eightBit && ((doAnimate == false) && (this->animate || specsChanged));
 
  //TODO: Remove force-true.
  const float time = doAnimate ? static_cast<float>(context->time->getAnimationTime()) : 7.f;
  if (animate || specsChanged) {
    if (workParallel || true) {
      const size_t numTasks = 1 + context->taskManager->getQueueConcurrency(TaskManager::GlobalQueue);
      const size_t subRange = N / numTasks;
 
      auto gr = context->taskManager->createGroup(TaskManager::GlobalQueue);
 
      for (size_t i = 0; i < numTasks; ++i) {
        context->taskManager->enqueueChild(gr, [&, i, time, N, subRange]()
                                           {
                                             //CpuInstrumentationScope("Systems", "OceanSystem::disperseWaves");
 
                                             disperseWaves(frqdDzdu, frqdDzdv, frqDx, frqDy, frqDz, frqH0, time, fftTileExtent, N, i * subRange, (i + 1) * subRange);
                                           });
      }
 
      context->taskManager->enqueueChild(oceanTaskGroup, [this, context, gr]()
                                         {
                                           context->taskManager->wait(gr);
 
                                           context->taskManager->enqueueChild(gr, [&, context]() { fastInverseFourierTransform2D(context, fftScratch1, spcdDzdu, frqdDzdu, 1.f, fftTileResolutionLog2); });
                                           context->taskManager->enqueueChild(gr, [&, context]() { fastInverseFourierTransform2D(context, fftScratch2, spcdDzdv, frqdDzdv, 1.f, fftTileResolutionLog2); });
                                           context->taskManager->enqueueChild(gr, [&, context]() { fastInverseFourierTransform2D(context, fftScratch3, spcDx, frqDx, 1.f, fftTileResolutionLog2); });
                                           context->taskManager->enqueueChild(gr, [&, context]() { fastInverseFourierTransform2D(context, fftScratch4, spcDy, frqDy, 1.f, fftTileResolutionLog2); });
                                           context->taskManager->enqueueChild(gr, [&, context]() { fastInverseFourierTransform2D(context, fftScratch5, spcDz, frqDz, 1.f, fftTileResolutionLog2); });
 
                                           context->taskManager->destroy(gr);
 
                                           updateTextures(magnitudes.channel0, magnitudes.channel1);
                                           if (this->extraStep) {
                                             LOG_DEBUG(logger, "Extra texture encoder run to get magnitudes right");
                                             updateTextures(magnitudes_tmp.channel0, magnitudes_tmp.channel1);
                                           }
                                         });
 
    }
    else {
      disperseWaves(frqdDzdu, frqdDzdv, frqDx, frqDy, frqDz, frqH0, time, fftTileExtent, N, 0, N);
      fastInverseFourierTransform2D(context, fftScratch, spcdDzdu, frqdDzdu, 1.f, fftTileResolutionLog2);
      fastInverseFourierTransform2D(context, fftScratch, spcdDzdv, frqdDzdv, 1.f, fftTileResolutionLog2);
      fastInverseFourierTransform2D(context, fftScratch, spcDx, frqDx, 1.f, fftTileResolutionLog2);
      fastInverseFourierTransform2D(context, fftScratch, spcDy, frqDy, 1.f, fftTileResolutionLog2);
      fastInverseFourierTransform2D(context, fftScratch, spcDz, frqDz, 1.f, fftTileResolutionLog2);
    }
  }
 
}
 
void Cogs::Core::BasicOceanSystem::postUpdate(Context * context)
{
  ComponentSystem::postUpdate(context);
 
  if (!pool.size()) return;
 
  // Update materials after new material instances has been instantiated.
  for (auto & oceanComp : pool) {
    auto & oceanData = getData(&oceanComp);
    auto gridComp = oceanComp.getComponent<AdaptivePlanarGridComponent>();
    auto & gridData = context->adaptivePlanarGridSystem->getData(gridComp);
 
    const TransformComponent* lodRefComp = nullptr;
    if (auto e = gridComp->lodReference.lock(); e) {
      lodRefComp = e->getComponent<TransformComponent>();
    }
    if (!lodRefComp) {
      lodRefComp = context->cameraSystem->getMainCamera()->getComponent<TransformComponent>();
    }
 
    bool any = false;
    glm::vec2 viewportSize;
    for (auto & weak : gridComp->cameras) {
      if (auto entity = weak.lock(); entity) {
        if (auto * comp = entity->getComponent<CameraComponent>(); comp) {
          auto & camData = context->cameraSystem->getData(comp);
          viewportSize = glm::max(viewportSize, camData.viewportSize);
          any = true;
        }
      }
    }
    if (!any) {
      viewportSize = context->cameraSystem->getMainCameraData().viewportSize;
    }
 
    updateTileMaterialInstances(oceanData,
                                gridComp,
                                gridData,
                                context->transformSystem->getLocalToWorld(lodRefComp),
                                viewportSize);
  }
 
  if (animate || specsChanged) {
    context->taskManager->wait(oceanTaskGroup);
    // Update 
    if (eightBit) {
      if (!std::isfinite(magnitudes.channel0) || !std::isfinite(magnitudes.channel1) || specsChanged || this->extraStep) {
        magnitudes.channel0 = magnitudes_tmp.channel0;
        magnitudes.channel1 = magnitudes_tmp.channel1;
      }
      else {
        magnitudes.channel0 = 0.9f * magnitudes.channel0 + 0.1f * magnitudes_tmp.channel0;
        magnitudes.channel1 = 0.9f * magnitudes.channel1 + 0.1f * magnitudes_tmp.channel1;
      }
    }
  }
  specsChanged = false;
 
}