/********************************************************************************/
/*                                                                              */
/*    Gerstner Wave Water Animation                                            */
/*    Ported/adapted from Noble Shaders (https://github.com/BelmuTM/Noble)     */
/*    Copyright (C) 2026  Belmu - licensed under the GNU General Public        */
/*    License v3.0. See LUT_LICENSE.txt for the full license text.             */
/*                                                                              */
/*    Adapted to run inside Unbound's material pipeline (view-space normals,   */
/*    playerPos/cameraPosition world position convention, GERSTNER_WAVES       */
/*    option block in lib/common.glsl).                                       */
/*                                                                              */
/********************************************************************************/

#ifdef GERSTNER_WAVES

#ifndef PI
    #define PI 3.14159265358979323846
#endif
#define GW_TAU 6.28318530717958647692
#define GW_EPS 1e-5
const float gwG = 9.81; // Earth's gravitational constant

// ---- Minimal self-contained noise/FBM (from Noble's include/utility/rng.glsl) ----
float gwHash12(vec2 p) {
    return fract(sin(dot(p, vec2(127.1, 311.7))) * 43758.5453123);
}

float gwNoise(vec2 uv) {
    vec2 ip = floor(uv);
    vec2 u  = fract(uv);
    u = u * u * (3.0 - 2.0 * u);

    float res = mix(
        mix(gwHash12(ip), gwHash12(ip + vec2(1.0, 0.0)), u.x),
        mix(gwHash12(ip + vec2(0.0, 1.0)), gwHash12(ip + vec2(1.0, 1.0)), u.x), u.y);
    return res * res;
}

const float gwFbmLacunarity  = 2.0;
const float gwFbmPersistance = 0.5;

float gwFBM(vec2 uv, int octaves, float frequency) {
    float height    = 0.0;
    float amplitude = 1.0;

    for (int i = 0; i < octaves; i++) {
        height    += gwNoise(uv * frequency) * amplitude;
        frequency *= gwFbmLacunarity;
        amplitude *= gwFbmPersistance;
    }
    return height;
}

// ---- Fast sqrt + rotation helpers (from Noble's include/utility/math.glsl) ----
float gwSqrtNewton(float x, float guess) { return 0.5 * (guess + x / guess); }
float gwFastSqrt(float x)  { return uintBitsToFloat((floatBitsToUint(x) >> 1) + 0x1FC00000u); }
float gwFastSqrtN1(float x) { return gwSqrtNewton(x, gwFastSqrt(x)); }

vec3 gwRotate(vec3 vector, vec3 from, vec3 to) {
    float cosTheta = dot(from, to);
    if (abs(cosTheta) >= 0.9999) { return cosTheta < 0.0 ? -vector : vector; }
    vec3 axis = normalize(cross(from, to));

    vec2 sc = vec2(sqrt(1.0 - cosTheta * cosTheta), cosTheta);
    return sc.y * vector + sc.x * cross(axis, vector) + (1.0 - sc.y) * dot(axis, vector) * axis;
}

// ---- Gerstner waves (from Noble's include/fragment/water.glsl) ----
#define GW_SETUP()                                                              \
    float gwSpeed      = WAVE_SPEED;                                            \
    float gwSteepness  = WAVE_STEEPNESS;                                        \
    float gwAmplitude  = WAVE_AMPLITUDE;                                        \
    float gwWavelength = WAVE_LENGTH * 1.0;                                     \
    float gwTime        = frameTimeCounter * gwSpeed;                          \
                                                                                 \
    const float gwAngle    = radians(WAVE_ANGLE);                               \
    const mat2 gwRotation  = mat2(cos(gwAngle), -sin(gwAngle), sin(gwAngle), cos(gwAngle)); \
                                                                                 \
    vec2 gwDirection = vec2(0.1, 0.1);                                          \
                                                                                 \
    float gwFbmNoise = gwFBM(position * 4e-2, 3, 1.0);

#define GW_PARAMS_FACTOR()                       \
    gwSteepness  *= WAVE_STEEPNESS_MULTIPLIER;   \
    gwAmplitude  *= WAVE_AMPLITUDE_MULTIPLIER;   \
    gwWavelength *= WAVE_LENGTH_MULTIPLIER;      \
    gwTime       *= WAVE_TIME_MULTIPLIER;        \
    gwDirection  *= gwRotation;

#define GW_TIME_NOISE() \
    gwTime + gwFbmNoise * 5.0

float gwGerstnerWaves(vec2 coords, float time, float steepness, float amplitude, float lambda, vec2 direction) {
    direction = normalize(direction);

    float k = GW_TAU / lambda;
    float x = gwFastSqrtN1(gwG * k) * time - k * dot(direction, coords);

    return amplitude * pow(sin(x) * 0.5 + 0.5, steepness);
}

vec2 gwGerstnerWavesDerivative(vec2 coords, float time, float steepness, float amplitude, float lambda, vec2 direction) {
    direction = normalize(direction);

    float k = GW_TAU / lambda;
    float x = gwFastSqrtN1(gwG * k) * time - k * dot(direction, coords);

    float u    = sin(x) * 0.5 + 0.5;
    float dudx = -0.5 * cos(x) * k;

    return amplitude * steepness * pow(u, steepness - 1.0) * dudx * direction;
}

vec2 calculateWaveDerivativeGerstner(vec2 position, int octaves) {
    vec2 derivative = vec2(0.0);

    GW_SETUP();

    for (int i = 0; i < octaves; i++) {

        derivative += gwGerstnerWavesDerivative(
            position,
            GW_TIME_NOISE(),
            gwSteepness,
            gwAmplitude,
            gwWavelength,
            gwDirection
        );

        GW_PARAMS_FACTOR();
    }

    return derivative;
}

// Returns a TANGENT-SPACE pseudo normal-map sample (x, y, z) built from the
// Gerstner wave slope, meant as a drop-in replacement for the texture-sampled
// `normalMap` Unbound normally reads from gaux4 - the caller still applies
// its own `tbnMatrix` afterwards exactly like the original texture path does,
// so it stays consistent with the rest of Unbound's water shading/fresnel code.
vec3 gerstnerWaterNormalTangentSpace(vec2 worldPositionXZ, int octaves, float strength, out vec2 waveDerivativeOut) {
    vec2 waveDerivative = calculateWaveDerivativeGerstner(worldPositionXZ, octaves) * strength;

    // v8.0: soft-normalize the multi-octave slope sum so high-frequency
    // octaves can't push the summed slope past ~1 and saturate the normal
    // map (z clamps to 0). Saturation is what made the surface read as
    // dense, distorted and shimmering instead of a smooth swell - bounding
    // it keeps the waves full but stable at every viewing angle, and smooth
    // even in small water bodies like cauldrons / water seen through glass.
    waveDerivative *= 1.0 / (1.0 + length(waveDerivative));
    waveDerivativeOut = waveDerivative;

    vec3 normalMapM;
    normalMapM.xy = -waveDerivative;
    normalMapM.z  = sqrt(max0(1.0 - dot(normalMapM.xy, normalMapM.xy)));
    return normalMapM;
}

#endif // GERSTNER_WAVES
