/*
    Procedural ocean ("Seascape" method): domain-warped choppy octaves for
    the wave field, plus fully analytic caustics derived from the same wave
    trains. Requires lib/common.glsl included first.
*/

// One sea octave: choppy waves on a noise-warped domain (curved, irregular
// wave fronts; no diamond artifacts)
float seaOctave(vec2 uv, float choppy) {
    uv += vnoise(uv) * 2.0 - 1.0;
    vec2 wv = 1.0 - abs(sin(uv));
    vec2 swv = abs(cos(uv));
    wv = mix(wv, swv, wv);
    return pow(1.0 - pow(wv.x * wv.y, 0.65), choppy);
}

// ============================ WAVES_V2 ============================
// Coherent directional spectrum: dominant ANALYTIC Gerstner-lite swell
// (long wavelength, fixed wind direction) + wind-aligned anisotropic chop
// (waveChopLOD). Total height = chop + swell; normal = finite differences
// on the CHOP only + ANALYTIC swell slope, barely attenuated at distance
// (the far ocean keeps rolling). SELF-CONTAINED: needs only common.glsl +
// settings.glsl (gbuffers_water.vsh includes water.glsl WITHOUT
// lighting.glsl - no uniforms here, time t is passed as an argument).
#if WAVES_V2 == 1
// Swell: height AND slope share one evaluator (declared BEFORE any use).
// Stokes-lite trains (cos + in-phase 2nd harmonic: pinched crests, wide
// troughs, ZERO mean -> shifts neither displacement nor the whitecap
// threshold); wave groups from very low-frequency vnoise (group gradient
// neglected: ~0.01 cycle/block, invisible in slope).
// WAVE_ITERATIONS swell trains: the 1st at full strength (lambda ~26), each
// next one disperses (+-30 deg around the wind, sin(i*2.4)*0.52), shortens
// (x0.58), weakens (x0.55) and speeds up (x1.35) -> rich crossed sea.
// Non-parallel fronts via domain warp (~9 blocks on ~50-block noise: crests
// undulate instead of being straight lines). The warp gradient is neglected
// in the slope (slow) - same approximation as the groups. The analytic
// caustics below derive from the SAME train set and follow these waves by
// construction: KEEP BOTH SUMS SYNCHRONIZED.
void swellEval(vec2 p, float t, out float height, out vec2 slope) {
    float ang = WIND_DIRECTION * 0.017453293;
    float amp = 0.38 * SWELL_AMOUNT; // swell IS the wave; chop is only detail on top
    float lambda = 26.0;
    float spd = 0.60 * WAVE_SPEED;
    float grp = 0.62 + 0.38 * vnoise(p * 0.013 + vec2(t * 0.009, -t * 0.007));
    vec2 wp2 = p + (vec2(vnoise(p * 0.020 + 3.7), vnoise(p * 0.020 + 9.1)) - 0.5) * 6.5;
    float q = 0.30; // Gerstner crest pinch (lower = softer crests)
    height = 0.0;
    slope = vec2(0.0);
    for (int i = 0; i < WAVE_ITERATIONS; i++) {
        // Directional spectrum: +-30 deg spread around the WIND. i=0 is
        // dead downwind; the sea rolls one way with variety, like a real
        // sea state (wide crossing angles made it unreadable).
        float aI = ang + sin(float(i) * 2.4) * 0.52;
        vec2 d = vec2(cos(aI), sin(aI));
        float k = 6.2831853 / lambda;
        float ph = dot(wp2, d) * k - t * spd;
        float a = amp * grp;
        height += a * (cos(ph) + q * cos(2.0 * ph));
        slope  += d * (a * k * (-sin(ph) - 2.0 * q * sin(2.0 * ph)));
        amp *= 0.60;
        lambda *= 0.58;
        spd *= 1.35;
    }
}
float swellHeight(vec2 p, float t) {
    float h; vec2 s;
    swellEval(p, t, h, s);
    return h;
}
vec2 swellSlope(vec2 p, float t) {
    float h; vec2 s;
    swellEval(p, t, h, s);
    return s;
}
#endif

// CHOP height, variable detail level (oct = 3 far ... 5 near).
// WAVES_V2=1: chop is ANISOTROPIC in wind space (field varies faster downwind
// than across = crests stretched perpendicular to the wind), 2 trains both
// drifting DOWNWIND, reduced octave rotation (~16 deg vs ~37) so the alignment
// survives across octaves. WAVES_V2=0 keeps the legacy isotropic body.
float waveChopLOD(vec2 p, float t, int oct) {
#if WAVES_V2 == 1
    float ang = WIND_DIRECTION * 0.017453293;
    vec2 wx = vec2(cos(ang), sin(ang));
    vec2 pw = vec2(dot(p, wx), dot(p, vec2(-wx.y, wx.x)));
    vec2 uv = vec2(pw.x * 1.06, pw.y * 0.78); // wind-aligned anisotropy
    // Tuning knobs: freq = base shape size (0.115 ~ 9 blocks), amp/choppy =
    // crest sharpness, decay 0.22 = how fast fine octaves die out. Raising
    // freq or the decay brings back a rough "gravel" look.
    float freq = 0.115;
    float amp = 0.40;
    float choppy = 2.4;
    float tt = 1.0 + t * 0.8 * WAVE_SPEED;

    float h = 0.0;
    for (int i = 0; i < 5; i++) {
        if (i >= oct) break;
        float d = seaOctave((uv - vec2(tt, tt * 0.14)) * freq, choppy)
                + seaOctave((uv - vec2(tt * 0.55, -tt * 0.09)) * freq * 1.13, choppy);
        h += d * amp;
        uv = mat2(1.92, 0.55, -0.55, 1.92) * uv;
        freq *= 1.9;
        amp *= 0.22;
        choppy = mix(choppy, 1.0, 0.3);
    }
    return h;
#else
    float freq = 0.14;
    float amp = 0.62;
    float choppy = 4.0;
    vec2 uv = p;
    uv.x *= 0.75;
    float tt = 1.0 + t * 0.8 * WAVE_SPEED;

    float h = 0.0;
    for (int i = 0; i < 5; i++) {
        if (i >= oct) break;
        float d = seaOctave((uv + tt) * freq, choppy)
                + seaOctave((uv - tt * 0.8) * freq, choppy);
        h += d * amp;
        uv = mat2(1.6, 1.2, -1.2, 1.6) * uv;
        freq *= 1.9;
        amp *= 0.28;
        choppy = mix(choppy, 1.0, 0.2);
    }
    return h;
#endif
}

// Total sea height = chop (+ swell when WAVES_V2). KEEP THIS SIGNATURE:
// consumers are gbuffers_water.vsh (displacement), composite1 (whitecaps),
// gbuffers_water.fsh (normals).
float waveHeightLOD(vec2 p, float t, int oct) {
#if WAVES_V2 == 1
    return waveChopLOD(p, t, oct) + swellHeight(p, t);
#else
    return waveChopLOD(p, t, oct);
#endif
}

float waveHeight(vec2 p, float t) {
    return waveHeightLOD(p, t, 5);
}

#if WAVE_NORMALS_V3 == 1
// Height field used for NORMALS ONLY (two-texture approach): layer A broad
// (~3.3 blocks) drifting downwind, layer B medium (~1.4 blocks) at ~22 deg
// with slight lateral counter-drift -> soft interference, round lobes, no
// choppy crests. Displacement (waveHeightLOD), whitecaps and caustics keep
// using their own fields.
float waterTexV3(vec2 p, float t) {
    float ang = WIND_DIRECTION * 0.017453293;
    vec2 wdir = vec2(cos(ang), sin(ang));
    vec2 pw = vec2(dot(p, wdir), dot(p, vec2(-wdir.y, wdir.x)) * 0.82);
    vec2 uvA = pw * 0.30 - vec2(t * 0.16 * WAVE_SPEED, 0.0);
    float hA = vnoise(uvA) * 0.72 + vnoise(uvA * 2.15 + 7.3) * 0.28;
    vec2 pb = mat2(0.927, 0.375, -0.375, 0.927) * pw;
    vec2 uvB = pb * 0.72 - vec2(t * 0.34 * WAVE_SPEED, -t * 0.05 * WAVE_SPEED);
    float hB = vnoise(uvB) * 0.70 + vnoise(uvB * 2.3 + 3.1) * 0.30;
    // Layer C: ripple built from seaOctave (self-warped domain -> organic
    // curved fronts, no perlin grid, soft choppy profile). 2 scales (~1.3 and
    // ~0.65 block) drifting downwind = the dense living surface texture.
    vec2 uvC = vec2(pw.x * 2.2, pw.y * 2.9) - vec2(t * 0.9 * WAVE_SPEED, t * 0.13 * WAVE_SPEED);
    float hC = seaOctave(uvC, 1.8) * 0.62
             + seaOctave(uvC * 2.05 + 3.7, 1.6) * 0.38;
    return hA * 0.50 + hB * 0.24 + hC * 0.26;
}
#endif

// World-space normal by finite differences, detail fading with distance
vec3 waveNormalWorld(vec2 p, float t, float strength, float chopFade) {
#if WAVE_NORMALS_V3 == 1
    // Soft 2-layer normals; the analytic WAVES_V2 swell (~26-block rollers) is
    // added on top -> the ocean ROLLS while the surface stays SMOOTH (clean
    // reflections, tidy sunset glint).
    // epsilon 0.22 resolves the ~0.65-block seaOctave ripple; slope gain 0.62
    float e = 0.22;
    float h0 = waterTexV3(p, t);
    float hx = waterTexV3(p + vec2(e, 0.0), t);
    float hz = waterTexV3(p + vec2(0.0, e), t);
    float slopeAtt3 = mix(0.45, 1.0, clamp(chopFade, 0.0, 1.0));
    vec2 slope = vec2(h0 - hx, h0 - hz) / e * strength * slopeAtt3 * 0.62;
    #if WAVES_V2 == 1
        float swellAtt3 = mix(0.75, 1.0, clamp(chopFade, 0.0, 1.0));
        slope -= swellSlope(p, t) * strength * swellAtt3;
    #endif
    #ifdef WATER_MICRONORMAL
        // Fine grain kept but reduced (x0.10) so the surface stays smooth
        float mf3 = smoothstep(0.45, 0.95, chopFade);
        if (mf3 > 0.001) {
            float md3 = 0.06;
            vec2 q3 = p * 3.7 + vec2(t * 0.45, -t * 0.33);
            float n03 = vnoise(q3);
            float nx3 = vnoise(q3 + vec2(md3, 0.0));
            float nz3 = vnoise(q3 + vec2(0.0, md3));
            slope += vec2((n03 - nx3) / md3, (n03 - nz3) / md3) * strength * mf3 * 0.10;
        }
    #endif
    return normalize(vec3(slope.x, 1.0, slope.y));
#else
    // Legacy normals. Octave count 4/3/2 by distance, widened epsilon and a
    // chop slope gain capped at 0.85 keep the micro-chop from resolving into
    // a gravel-like grain; the swell (x1.25 below) draws the surface shape.
    int oct = (chopFade > 0.5) ? 4 : ((chopFade > 0.15) ? 3 : 2);
    float e = mix(0.16, 0.11, clamp(chopFade, 0.0, 1.0));
    // Finite differences on the CHOP only: the long swell is added
    // ANALYTICALLY below (a short epsilon cannot resolve a 26-block wave).
    float h0 = waveChopLOD(p, t, oct);
    float hx = waveChopLOD(p + vec2(e, 0.0), t, oct);
    float hz = waveChopLOD(p + vec2(0.0, e), t, oct);
    // Slope attenuation with distance: kills far specular shimmer
    float slopeAtt = mix(0.30, 0.85, clamp(chopFade, 0.0, 1.0));
    vec2 slope = vec2((h0 - hx) / e, (h0 - hz) / e) * strength * slopeAtt;
    #if WAVES_V2 == 1
        // ANALYTIC swell slope, barely attenuated at distance so the far ocean
        // keeps rolling instead of flattening. Sign convention: slope = -dh/dxz
        // (final normal is (slope.x, 1.0, slope.y)).
        float swellAtt = mix(0.75, 1.0, clamp(chopFade, 0.0, 1.0));
        slope -= swellSlope(p, t) * strength * swellAtt * 1.25;
    #endif
    #ifdef WATER_MICRONORMAL
        // High-frequency grain, near-field only (faded by chopFade)
        float mf = smoothstep(0.45, 0.95, chopFade);
        if (mf > 0.001) {
            float md = 0.06;
            vec2 q = p * 3.7 + vec2(t * 0.45, -t * 0.33);
            float n0 = vnoise(q);
            float nx = vnoise(q + vec2(md, 0.0));
            float nz = vnoise(q + vec2(0.0, md));
            slope += vec2((n0 - nx) / md, (n0 - nz) / md) * strength * mf * 0.07; // faint surface shiver; above ~0.12 it reads as rock
        }
    #endif
    return normalize(vec3(slope.x, 1.0, slope.y));
#endif
}

// Animated Voronoi edges (cell points dancing).
// NOT used by the caustics anymore (the cellular net was decorative: it did
// not derive from the wave field). Kept as a generic building block.
float voroEdge(vec2 p, float t) {
    vec2 ip = floor(p);
    vec2 fp = fract(p);
    float f1 = 8.0;
    float f2 = 8.0;
    for (int x = -1; x <= 1; x++) {
        for (int y = -1; y <= 1; y++) {
            vec2 g = vec2(float(x), float(y));
            vec2 o = vec2(hash12(ip + g), hash12(ip + g + 19.7));
            o = 0.5 + 0.42 * sin(t * 0.7 + o * 6.2831853);
            float d = length(g + o - fp);
            if (d < f1) { f2 = f1; f1 = d; }
            else if (d < f2) { f2 = d; }
        }
    }
    return f2 - f1;
}

/*
    ==================== FULLY ANALYTIC CAUSTICS ====================
    A caustic is the convergence of rays refracted by a SUM OF COSINE WAVE
    TRAINS. The Hessian of a sum of cosines is exact and smooth by
    construction, so there is no finite differencing and no noise. Two train
    families:
      - SWELL: the WAVE_ITERATIONS trains, constants STRICTLY in sync with
        swellEval (amp 0.38*SWELL_AMOUNT, lambda 26, decay 0.60/0.58/1.35,
        q=0.30, spread ang+sin(i*2.4)*0.52, spd 0.60*WAVE_SPEED, warp wp2,
        groups grp) -> the broad pattern rolls with the RENDERED swell.
      - 6 SHORT trains (lambda 5.30/3.70/2.60/1.90/1.35/0.95, NON harmonic
        ratios): the fine net. Directions spread around the wind (alternating
        offsets, +-~52 deg, no two neighbours aligned), distinct initial
        phases, steepness a*k = 0.16..0.10 as a decreasing spectrum, speeds
        from GRAVITY DISPERSION omega ~ sqrt(k) (calibrated on the reference
        train: lambda 26 at 0.60*WAVE_SPEED) -> ripples oscillate faster than
        the swell, as in reality.
    ANTI-GRID: well-spread directions + non-harmonic lambdas + offset phases
    + shared slow warp wp2 (curved fronts) + very-low-frequency PER-TRAIN
    amplitude modulation (gradient neglected, like the swell groups) + the
    swell on top -> rich living interference, no visible lattice.
    EXACT SPECTRAL BLUR: gaussian-blurring the field with radius sigma is
    exactly weighting each train by exp(-0.5*(k*sigma)^2) (Fourier transform
    of the kernel; the filtering happens IN the spectrum). Zero extra taps:
    with depth/distance the fine filaments MELT into wide soft blobs instead
    of switching off.
    PERF: one Hessian per call = WAVE_ITERATIONS+6 trains of a few cos/sin
    each. CAUSTIC_SCALE is INERT (caustic scale IS the wave scale). Pure ALU,
    no uniforms, no 430-only functions: safe in both v120 (gbuffers_water.vsh)
    and v430 (composite1).
*/

// SWELL trains -> Hessian accumulation. h = a*(cos(ph) + q*cos(2ph)) along d
// gives d2h/dd2 = -a*k^2*(cos + 4q*cos2); in world space the contribution is
// (d outer d) * c2. SAME constants as swellEval - KEEP IN SYNC (at
// CAUSTIC_SPEED=1 the phases match the rendered surface exactly).
// sigma = spectral blur radius (see header). Warp and groups are passed in
// (shared with the short trains; their gradients are neglected, they are slow).
void causticSwellH(vec2 wp2, float grp, float t, float sigma,
                   inout float hxx, inout float hzz, inout float hxz) {
    float ang = WIND_DIRECTION * 0.017453293;
    float amp = 0.38 * SWELL_AMOUNT; // sync swellEval
    float lambda = 26.0;
    float spd = 0.60 * WAVE_SPEED;   // sync swellEval - holds at any WAVE_SPEED
    float q = 0.30;                  // sync swellEval (pinched crests)
    for (int i = 0; i < WAVE_ITERATIONS; i++) {
        float aI = ang + sin(float(i) * 2.4) * 0.52; // sync directional spectrum
        vec2 d = vec2(cos(aI), sin(aI));
        float k = 6.2831853 / lambda;
        float ph = dot(wp2, d) * k - t * spd;
        float w1 = exp(-0.5 * k * k * sigma * sigma); // exact spectral blur (fundamental, wavenumber k)
        float w2 = w1 * w1 * w1 * w1;                 // 2nd Stokes harmonic lives at 2k -> exp(-0.5*(2k*sigma)^2) = w1^4
        float c2 = -(amp * grp) * k * k * (cos(ph) * w1 + 4.0 * q * cos(2.0 * ph) * w2);
        hxx += d.x * d.x * c2;
        hzz += d.y * d.y * c2;
        hxz += d.x * d.y * c2;
        amp *= 0.60;
        lambda *= 0.58;
        spd *= 1.35;
    }
}

// One SHORT train: pure cosine (short ripples are nearly sinusoidal, and it
// is cheaper). steep = a*k (steepness) -> the Hessian contribution is
// -steep*k*cos(ph), the amplitude alone is never needed. The amplitude
// modulation (vnoise ~0.02-0.04/block + slow drift, gradient neglected) makes
// the net live region by region instead of looking statistically uniform.
// omg is multiplied by WAVE_SPEED here.
void causticTrainH(vec2 wp2, vec2 p, float t, float ang, float sigma,
                   float lam, float omg, float dOff, float ph0, float steep,
                   float mFrq, vec2 mOff,
                   inout float hxx, inout float hzz, inout float hxz) {
    float k = 6.2831853 / lam;
    float aT = ang + dOff;
    vec2 d = vec2(cos(aT), sin(aT));
    float ph = dot(wp2, d) * k - t * (omg * WAVE_SPEED) + ph0;
    float m = 0.55 + 0.90 * vnoise(p * mFrq + mOff + vec2(t * 0.011, -t * 0.008));
    float w = exp(-0.5 * k * k * sigma * sigma); // exact spectral blur
    float c2 = -steep * k * cos(ph) * m * w;
    hxx += d.x * d.x * c2;
    hzz += d.y * d.y * c2;
    hxz += d.x * d.y * c2;
}

// FULL caustic field: slow warp + groups (sync swellEval), swell, then the 6
// short trains. Columns: lambda | omega (= 0.60*sqrt(26/lambda), gravity
// dispersion) | direction offset (rad, alternating, no two neighbouring
// lambdas aligned) | initial phase | steepness a*k | frequency and offset of
// the amplitude modulation. All constant -> the compiler folds the direction
// cos/sin and the weights.
void causticFieldH(vec2 p, float t, float sigma,
                   out float hxx, out float hzz, out float hxz) {
    float ang = WIND_DIRECTION * 0.017453293;
    vec2 wp2 = p + (vec2(vnoise(p * 0.020 + 3.7), vnoise(p * 0.020 + 9.1)) - 0.5) * 6.5; // sync swellEval
    float grp = 0.62 + 0.38 * vnoise(p * 0.013 + vec2(t * 0.009, -t * 0.007));           // sync swellEval
    hxx = 0.0;
    hzz = 0.0;
    hxz = 0.0;
    causticSwellH(wp2, grp, t, sigma, hxx, hzz, hxz);
    //                                   lam    omega    dir    ph0   a*k    mFrq   mOff
    causticTrainH(wp2, p, t, ang, sigma, 5.30, 1.3289, -0.82, 1.7, 0.160, 0.021, vec2(11.3,  3.1), hxx, hzz, hxz);
    causticTrainH(wp2, p, t, ang, sigma, 3.70, 1.5905,  0.55, 4.2, 0.150, 0.026, vec2( 7.9, 21.4), hxx, hzz, hxz);
    causticTrainH(wp2, p, t, ang, sigma, 2.60, 1.8974, -0.28, 0.9, 0.140, 0.031, vec2( 2.6, 15.8), hxx, hzz, hxz);
    causticTrainH(wp2, p, t, ang, sigma, 1.90, 2.2195,  0.91, 5.5, 0.130, 0.036, vec2(17.2,  8.5), hxx, hzz, hxz);
    causticTrainH(wp2, p, t, ang, sigma, 1.35, 2.6331,  0.16, 2.6, 0.115, 0.024, vec2( 5.4, 27.7), hxx, hzz, hxz);
    causticTrainH(wp2, p, t, ang, sigma, 0.95, 3.1389, -0.61, 3.9, 0.100, 0.033, vec2(23.1, 12.9), hxx, hzz, hxz);
}

// Optics (paraxial refraction, eta = 1/1.33): a vertical refracted ray hits
// the bottom displaced by ~cf*grad(h), with cf ~ (1-eta)*column ~ 0.25*depth.
// Jacobian of the surface->bottom map: J = I + cf*H; received irradiance is
// 1/|det J|. det->0: rays CONVERGE (filament); det>1: divergence (penumbra
// between filaments). Floor 0.10 is the anti-firefly clamp (I caps at 10;
// the soft remap does the rest, no hard clamp here).
float causticIllum(float hxx, float hzz, float hxz, float cf) {
    float a11 = 1.0 + cf * hxx;
    float a22 = 1.0 + cf * hzz;
    float det = a11 * a22 - cf * cf * hxz * hxz;
    return 1.0 / max(abs(det), 0.10);
}

// SOFT anti-firefly remap, CALIBRATED against the fixed thresholds in
// composite1 (site 1: clamp(X*cscale - 0.22..0.30, 0, 1) with cscale
// 0.70..1.30; chroma sites: clamp(X - 0.20..0.25, 0, 1) * exp(-depth*0.13)).
// Reference points (keep them if you retune):
//   flat water det=1 -> I=1 -> base 0.36: clears the thresholds as a DISCREET
//   glow (site 1: 0.03..0.25 before gain) without ever blowing out sand;
//   divergence I<1 -> linear ramp 0.36*I down to 0 (DARK between filaments -
//   light is redistributed, not added); never negative;
//   convergence I ~ 1.6..2.5 (typical) -> 0.87..1.42: crisp filament that
//   saturates downstream at peak cscale;
//   above that: C1 exponential knee (slope continuous at x=1.05) toward the
//   1.65 asymptote - progressive compression, no hard plateau (a flat clamp
//   produces white blotches).
float causticTone(float I) {
    float x = 0.36 * min(I, 1.0) + 0.85 * max(I - 1.0, 0.0);
    float over = max(x - 1.05, 0.0);
    return min(x, 1.05) + 0.60 * (1.0 - exp(-over / 0.60));
}

// LARGE octave (god-rays / modulation): the SAME sum reduced to the SWELL
// trains, wide sigma (0.9) -> big smooth cells that roll with the swell.
// Normalization (floor 0.20, -0.9, saturate, then ^2) must stay: consumers
// (UW_GODRAY_CAUSTIC_LOCK shafts, the 0.55+0.90*large modulation in
// causticCombine) expect a [0..1] value.
float causticLarge(vec2 p, float t) {
    float tt = t * CAUSTIC_SPEED;
    vec2 wp2 = p + (vec2(vnoise(p * 0.020 + 3.7), vnoise(p * 0.020 + 9.1)) - 0.5) * 6.5;
    float grp = 0.62 + 0.38 * vnoise(p * 0.013 + vec2(tt * 0.009, -tt * 0.007));
    float hxx = 0.0;
    float hzz = 0.0;
    float hxz = 0.0;
    causticSwellH(wp2, grp, tt, 0.9, hxx, hzz, hxz);
    // Long focal length of the large octave, driven by CAUSTIC_LARGE_SCALE
    // (default 0.28 -> cf 1.35), floored at 0.3.
    float cf = 4.5 * max(CAUSTIC_LARGE_SCALE, 0.3);
    float a11 = 1.0 + cf * hxx;
    float a22 = 1.0 + cf * hzz;
    float det = a11 * a22 - cf * cf * hxz * hxz;
    float w = clamp(1.0 / max(abs(det), 0.20) - 0.9, 0.0, 1.0);
    return w * w;
}

// GOD-RAY APERTURE MASK - dedicated to the volumetric shafts (composite1,
// UW_GODRAY_CAUSTIC_LOCK==1). DO NOT wire it into the caustic render sites:
// causticLarge (determinant 1/|det J| -> convergence folds are thin LINES by
// nature) stays the source of waveCausticFull/.y. Here we want the opposite
// geometry: well-separated 6-15 block PATCHES with 5-20 blocks of penumbra.
// Principle: curvature (Laplacian hxx+hzz) of a VERY smooth swell field -
// spectral blur sigma grows with the water column (real defocus): only the
// lambda 26/15 trains survive, so the field varies over ~8-16 blocks. Under a
// crest DOME (negative Laplacian = converging lens) the whole sheet focuses
// -> a WIDE light column; under troughs, penumbra. Soft smoothstep threshold,
// no determinant.
// /SWELL_AMOUNT normalization keeps the column geometry independent of the
// amplitude slider (0 = flat sea -> null mask, consistent).
// Calibration (normalized, WAVE_ITERATIONS=3, mean grp 0.81): per-train
// curvature amplitudes ~0.030/0.034/0.024 -> RMS ~0.029; low threshold 0.015
// (~0.5 sigma, ~30% of the area lightly touched), full 0.055 (~1.9 sigma, ~3%
// = dome cores). At WATER_QUALITY 0 (2 trains) RMS ~0.024, slightly rarer
// columns. Animated by the SAME phases as the rendered swell (warp, groups,
// speeds in sync) so columns drift ~2.5 blocks/s and breathe with the wave
// groups. COST: identical to causticLarge - it is a replacement, not an add.
float causticShaftMask(vec2 p, float t, float depthM) {
    float tt = t * CAUSTIC_SPEED;
    vec2 wp2 = p + (vec2(vnoise(p * 0.020 + 3.7), vnoise(p * 0.020 + 9.1)) - 0.5) * 6.5; // sync swellEval
    float grp = 0.62 + 0.38 * vnoise(p * 0.013 + vec2(tt * 0.009, -tt * 0.007));         // sync swellEval
    float hxx = 0.0;
    float hzz = 0.0;
    float hxz = 0.0;
    // Base sigma 1.15 keeps the lambda~8.75 train alive (weight ~0.71) -> fine
    // 3-6 block BLADES elongated along the swell, spaced 4-10 blocks, rather
    // than round 8-16 block blobs. Raise it for wider, rounder patches.
    float sigma = (1.15 + min(depthM * 0.019, 0.45)) * UW_SHAFT_SCALE; // UW_SHAFT_SCALE = spacing between blades; depth term = soft defocus with the water column
    causticSwellH(wp2, grp, tt, sigma, hxx, hzz, hxz);
    float lap = -(hxx + hzz) / max(SWELL_AMOUNT, 0.25); // >0 under crests (converging lens)
    // Thresholds scaled by UW_SHAFT_FINESSE (default 1.3) -> 2-4 block shafts
    // separated by 10-20 blocks of penumbra.
    return smoothstep(0.022 * UW_SHAFT_FINESSE, 0.080 * UW_SHAFT_FINESSE, lap);
}

// SMALL caustic (high frequency). sharp in [0..1]: 1=crisp, 0=soft.
// KEEP THIS SIGNATURE (composite1 god-rays; sharp decreases with toSurf).
// No depth is passed, so a proxy column is used: crisp ~ 1.6 blocks (short
// ripple focal length), soft ~ 6 blocks. The blur lives in sigma (spectral).
float causticSmall(vec2 p, float t, float sharp) {
    // CAUSTIC_SPEED 1.0 = EXACT sync with the surface (the swell trains land
    // on swellEval's phases); != 1.0 deliberately desyncs the caustic from the
    // waves (stylized).
    float tt = t * CAUSTIC_SPEED;
    float s = clamp(sharp, 0.0, 1.0);
    float hxx, hzz, hxz;
    causticFieldH(p, tt, mix(0.85, 0.12, s), hxx, hzz, hxz); // sigma capped at 0.85 (see causticCombine)
    float cf = 0.25 * clamp(mix(6.0, 1.6, s), 1.2, 8.0) * CAUSTIC_FOCUS; // 1.2 floor matches causticCombine (inert here, mix >= 1.6)
    return causticTone(causticIllum(hxx, hzz, hxz, cf));
}
// Compat entry point, behaviour IDENTICAL to causticSmall(sharp=1).
float waveCaustic(vec2 p, float t) {
    return causticSmall(p, t, 1.0);
}

// Shared worker for Full/Soft: ONE field blurred by blurTotal (sigma =
// mix(0.12, 0.85, blurTotal)), focal length from the REAL depth (cf = 0.25 *
// column * CAUSTIC_FOCUS), large octave used as MODULATION (never additive).
// .x = combined net (render), .y = large octave alone (god-rays).
// The 0.85 sigma cap matters: at sigma 1.8 the short-train spectral weights
// collapse (0.10 at lambda 5.3, 0.009 at 3.7, ~0 below) leaving only the
// lambda-26 swell -> flat 0.36 output, under composite1's 0.20..0.30
// thresholds, i.e. nothing on screen. At 0.85: lambda 5.3 keeps w ~ 0.60, 3.7
// ~ 0.35, 2.6 ~ 0.12 -> credible defocus (wide soft dancing blobs) instead of
// extinction. Typical values at defaults (2-3 block column): filaments
// causticTone ~ 1.1..1.5, flat zone ~ 0.36.
vec2 causticCombine(vec2 p, float t, float depthM, float blurTotal) {
    float tt = t * CAUSTIC_SPEED; // 1.0 = exact sync (see causticSmall)
    float hxx, hzz, hxz;
    causticFieldH(p, tt, mix(0.12, 0.85, blurTotal), hxx, hzz, hxz);
    // Focal floor 1.2: the ripple focus sits at ~1 block of column, where real
    // filaments are already sharp; a lower floor makes them timid.
    float cf = 0.25 * clamp(depthM, 1.2, 8.0) * CAUSTIC_FOCUS;
    float fine = causticTone(causticIllum(hxx, hzz, hxz, cf));
    float large = causticLarge(p, t);
    // The large octave MODULATES the fine net around 1.0 (never an additive
    // cyan patch). lm saturates at 1.0 once CAUSTIC_LARGE_STRENGTH >= 2.0.
    float lm = clamp(CAUSTIC_LARGE_STRENGTH * 0.5, 0.0, 1.0);
    return vec2(fine * mix(1.0, 0.55 + 0.90 * large, lm), large);
}
// Depth defocus window. Real caustics are SHARPEST at 1-4 blocks of column
// (ripple focal length) and only defocus beyond ~8-15 m. So: full sharpness
// up to 2.5 blocks, then a fade over a window of 4x the slider (min 8 blocks)
// - at the 3.0 default, full blur around ~14.5 blocks, where the water already
// hides the bottom anyway. A window narrower than this puts full sigma inside
// the physically sharp zone (invisible caustics, milky bottom).
// Slider semantics: higher = later/softer defocus, lower = earlier/stronger.
float causticDepthBlur(float depthM) {
    return clamp((depthM - 2.5) / max(CAUSTIC_DEPTH_BLUR * 4.0, 8.0), 0.0, 1.0);
}
// Small + large + depth blur. depthM = water column in blocks.
// The depth blur removes the short trains SPECTRALLY (fine filaments melt
// into wide swell-borne blobs): the caustic widens and softens physically
// instead of fading out.
vec2 waveCausticFull(vec2 p, float t, float depthM) {
    float blurTotal = causticDepthBlur(depthM);
    return causticCombine(p, t, depthM, blurTotal);
}
// Distance fade-to-blur: like waveCausticFull but with blurX [0..1] driven by
// EYE DISTANCE. The 1-(1-a)(1-b) combination stays bounded by 1, so far away
// the caustic is WIDE and SOFT but still present. blurX=0 == Full exactly.
vec2 waveCausticSoft(vec2 p, float t, float depthM, float blurX) {
    float dB = causticDepthBlur(depthM);
    float blurTotal = 1.0 - (1.0 - dB) * (1.0 - clamp(blurX, 0.0, 1.0));
    return causticCombine(p, t, depthM, blurTotal);
}
// TRUE chromatic dispersion, nearly free: the index of water varies with
// wavelength -> the FOCAL LENGTH varies per channel. ONE field, 3 determinants
// with a scaled cf (x0.975 R / 1.0 G / x1.035 B ~ n(lambda)): near the focus
// (det->0) the per-channel gap blows up, giving a subtle spectral fringe on
// the EDGE of the filaments, exactly where it shows in reality.
// chromOff is IGNORED (legacy spatial offset of the voronoi net); dispersion
// now lives in the focal length. Signature kept for composite1.
vec3 waveCausticChroma(vec2 p, float t, float depthM, float blurX, vec2 chromOff) {
    float dB = causticDepthBlur(depthM);
    float blurTotal = 1.0 - (1.0 - dB) * (1.0 - clamp(blurX, 0.0, 1.0));
    float tt = t * CAUSTIC_SPEED;
    float hxx, hzz, hxz;
    causticFieldH(p, tt, mix(0.12, 0.85, blurTotal), hxx, hzz, hxz); // sigma capped at 0.85
    float cf = 0.25 * clamp(depthM, 1.2, 8.0) * CAUSTIC_FOCUS;      // focal floor 1.2
    float large = causticLarge(p, t);
    float lm = clamp(CAUSTIC_LARGE_STRENGTH * 0.5, 0.0, 1.0);
    float lgMod = mix(1.0, 0.55 + 0.90 * large, lm);
    vec3 disp = vec3(0.975, 1.0, 1.035); // n(red) < n(green) < n(blue)
    vec3 outC;
    for (int c = 0; c < 3; c++) {
        outC[c] = causticTone(causticIllum(hxx, hzz, hxz, cf * disp[c])) * lgMod;
    }
    return outC;
}
// Anti-repetition foam pattern: 2-pass domain warp + octave rotation.
float foamPattern(vec2 p, float t, float s) {
    p *= s;
    vec2 w1 = vec2(fbm(p * 0.35 + vec2(t * 0.05, -t * 0.04)),
                   fbm(p * 0.35 + vec2(5.2 - t * 0.045, 1.7 + t * 0.05)));
    p += (w1 - 0.5) * 2.4;
    p = mat2(0.80, 0.60, -0.60, 0.80) * p;
    vec2 w2 = vec2(fbm(p * 0.9 - vec2(t * 0.07, t * 0.05)),
                   fbm(p * 0.9 + vec2(t * 0.06, 3.1 - t * 0.055)));
    p += (w2 - 0.5) * 1.1;
    return fbm(p * 1.6 + t * 0.10);
}

// ================= CACHED CAUSTIC PATH (manual CSE) =================
// Same math as waveCausticFull/Soft, but the work the original path pays
// TWICE per pixel (once in causticFieldH, once in causticLarge, both with
// the SAME p and tt) is evaluated ONCE:
//   - wp2 slow warp (2 vnoise) and grp groups (1 vnoise): hoisted, 3 vnoise
//     saved per pixel;
//   - the WAVE_ITERATIONS swell trains: per train d, k, ph, cos(ph) and
//     cos(2*ph) are sigma-independent; only the spectral weights w1/w2
//     differ between the fine sigma and the large sigma (0.9). One loop
//     accumulates BOTH Hessians -> WAVE_ITERATIONS train phase/trig
//     evaluations saved per pixel.
// Every arithmetic expression matches the originals term for term (same
// parse tree), so the output is bit-identical -- no gate needed. Consumed
// by composite1 water site 1 only; the UW caustic sites and the VL keep
// the original entry points.

// Dual-sigma swell Hessian: one phase evaluation, two spectral weightings.
// Constants sync swellEval / causticSwellH - KEEP ALL THREE IN SYNC.
void causticSwellH2_cached(vec2 wp2, float grp, float t, float sigmaA, float sigmaB,
                           inout float hxxA, inout float hzzA, inout float hxzA,
                           inout float hxxB, inout float hzzB, inout float hxzB) {
    float ang = WIND_DIRECTION * 0.017453293;
    float amp = 0.38 * SWELL_AMOUNT; // sync swellEval
    float lambda = 26.0;
    float spd = 0.60 * WAVE_SPEED;   // sync swellEval
    float q = 0.30;                  // sync swellEval (pinched crests)
    for (int i = 0; i < WAVE_ITERATIONS; i++) {
        float aI = ang + sin(float(i) * 2.4) * 0.52; // sync directional spectrum
        vec2 d = vec2(cos(aI), sin(aI));
        float k = 6.2831853 / lambda;
        float ph = dot(wp2, d) * k - t * spd;
        // Shared per-train terms (sigma-independent)
        float cp  = cos(ph);
        float cp2 = cos(2.0 * ph);
        float base = -(amp * grp) * k * k;
        // Per-sigma spectral weights (2nd harmonic at 2k -> w^4, see causticSwellH)
        float w1A = exp(-0.5 * k * k * sigmaA * sigmaA);
        float w2A = w1A * w1A * w1A * w1A;
        float w1B = exp(-0.5 * k * k * sigmaB * sigmaB);
        float w2B = w1B * w1B * w1B * w1B;
        float cA = base * (cp * w1A + 4.0 * q * cp2 * w2A);
        float cB = base * (cp * w1B + 4.0 * q * cp2 * w2B);
        hxxA += d.x * d.x * cA;
        hzzA += d.y * d.y * cA;
        hxzA += d.x * d.y * cA;
        hxxB += d.x * d.x * cB;
        hzzB += d.y * d.y * cB;
        hxzB += d.x * d.y * cB;
        amp *= 0.60;
        lambda *= 0.58;
        spd *= 1.35;
    }
}

// Large octave from a PRECOMPUTED swell Hessian (math sync causticLarge).
float causticLargeH_cached(float hxx, float hzz, float hxz) {
    float cf = 4.5 * max(CAUSTIC_LARGE_SCALE, 0.3); // sync causticLarge
    float a11 = 1.0 + cf * hxx;
    float a22 = 1.0 + cf * hzz;
    float det = a11 * a22 - cf * cf * hxz * hxz;
    float w = clamp(1.0 / max(abs(det), 0.20) - 0.9, 0.0, 1.0);
    return w * w;
}

// Cached worker: math sync causticCombine (fine sigma mix(0.12,0.85), large
// sigma 0.9, focal floor 1.2, large-octave modulation) and causticFieldH
// (short train table VERBATIM). wp2/grp computed once, swell trains walked
// once for both Hessians.
vec2 causticCombine_cached(vec2 p, float t, float depthM, float blurTotal) {
    float tt = t * CAUSTIC_SPEED; // 1.0 = exact sync (see causticSmall)
    float ang = WIND_DIRECTION * 0.017453293;
    vec2 wp2 = p + (vec2(vnoise(p * 0.020 + 3.7), vnoise(p * 0.020 + 9.1)) - 0.5) * 6.5; // sync swellEval
    float grp = 0.62 + 0.38 * vnoise(p * 0.013 + vec2(tt * 0.009, -tt * 0.007));         // sync swellEval
    float sigma = mix(0.12, 0.85, blurTotal); // sigma capped at 0.85 (see causticCombine)
    float hxx = 0.0; float hzz = 0.0; float hxz = 0.0; // fine field
    float lxx = 0.0; float lzz = 0.0; float lxz = 0.0; // large octave
    causticSwellH2_cached(wp2, grp, tt, sigma, 0.9, hxx, hzz, hxz, lxx, lzz, lxz);
    //                                    lam    omega    dir    ph0   a*k    mFrq   mOff   (VERBATIM causticFieldH)
    causticTrainH(wp2, p, tt, ang, sigma, 5.30, 1.3289, -0.82, 1.7, 0.160, 0.021, vec2(11.3,  3.1), hxx, hzz, hxz);
    causticTrainH(wp2, p, tt, ang, sigma, 3.70, 1.5905,  0.55, 4.2, 0.150, 0.026, vec2( 7.9, 21.4), hxx, hzz, hxz);
    causticTrainH(wp2, p, tt, ang, sigma, 2.60, 1.8974, -0.28, 0.9, 0.140, 0.031, vec2( 2.6, 15.8), hxx, hzz, hxz);
    causticTrainH(wp2, p, tt, ang, sigma, 1.90, 2.2195,  0.91, 5.5, 0.130, 0.036, vec2(17.2,  8.5), hxx, hzz, hxz);
    causticTrainH(wp2, p, tt, ang, sigma, 1.35, 2.6331,  0.16, 2.6, 0.115, 0.024, vec2( 5.4, 27.7), hxx, hzz, hxz);
    causticTrainH(wp2, p, tt, ang, sigma, 0.95, 3.1389, -0.61, 3.9, 0.100, 0.033, vec2(23.1, 12.9), hxx, hzz, hxz);
    float cf = 0.25 * clamp(depthM, 1.2, 8.0) * CAUSTIC_FOCUS; // focal floor 1.2 (sync causticCombine)
    float fine = causticTone(causticIllum(hxx, hzz, hxz, cf));
    float large = causticLargeH_cached(lxx, lzz, lxz);
    float lm = clamp(CAUSTIC_LARGE_STRENGTH * 0.5, 0.0, 1.0);
    return vec2(fine * mix(1.0, 0.55 + 0.90 * large, lm), large);
}

// Drop-in cached entry points -- output identical to waveCausticFull/Soft.
vec2 waveCausticFull_cached(vec2 p, float t, float depthM) {
    float blurTotal = causticDepthBlur(depthM);
    return causticCombine_cached(p, t, depthM, blurTotal);
}
vec2 waveCausticSoft_cached(vec2 p, float t, float depthM, float blurX) {
    float dB = causticDepthBlur(depthM);
    float blurTotal = 1.0 - (1.0 - dB) * (1.0 - clamp(blurX, 0.0, 1.0));
    return causticCombine_cached(p, t, depthM, blurTotal);
}

// ================= FAST WATER NORMALS (WAVE_FAST_NORMALS) =================
// ADDITIVE fast path for the legacy Seascape v2 normals (WAVES_V2=1,
// WAVE_NORMALS_V3=0). Two savings over waveNormalWorld:
//   (a) manual CSE: the wind-space transform of waveChopLOD (rotation +
//       anisotropy) is LINEAR in p, so it is applied ONCE and the two
//       finite-difference offsets become constant vectors in wind space
//       (the per-octave warps/rotations depend on the SHIFTED uv and stay
//       per-call - they cannot be fused);
//   (b) tighter octave LOD, measured at 1080p (fov 70 -> focal ~771 px, eye
//       ~1.8 bl above water; ground footprint per pixel along view ~
//       d^2/(771*1.8)). Per octave the uv domain scales x~2.0 AND freq x1.9,
//       so wavelengths are ~8.7 / 2.3 / 0.61 / 0.16 bl. Nyquist (lambda/2 per
//       pixel): octave 3 is sub-pixel past ~11 bl, octave 2 past ~21-45 bl.
//       New tiers: 4 octaves < ~17 bl (chopFade > 0.60), 3 < ~58 (> 0.30),
//       2 < ~288 (> 0.08), 1 beyond (slopeAtt floor 0.30 hides the residue).
//       Legacy tiers were 4 < 25 bl and 3 < 142 bl.
// NOTE a 2-call chop gradient is NOT offered: the field is scalar and a 2D
// gradient needs 3 samples minimum (forward differences are already in use;
// "central on one axis" would ADD a 4th sample, not remove one).
// Epsilon, slope attenuation, swell and micro-grain match waveNormalWorld
// term for term - at equal octave count the output is bit-identical.
#if WAVE_FAST_NORMALS == 1 && WAVES_V2 == 1
// Chop with PRE-TRANSFORMED wind-space uv. Body sync waveChopLOD (WAVES_V2
// branch) minus the transform - KEEP IN SYNC.
float waveChopUVFast(vec2 uv, float t, int oct) {
    float freq = 0.115;
    float amp = 0.40;
    float choppy = 2.4;
    float tt = 1.0 + t * 0.8 * WAVE_SPEED;

    float h = 0.0;
    for (int i = 0; i < 5; i++) {
        if (i >= oct) break;
        float d = seaOctave((uv - vec2(tt, tt * 0.14)) * freq, choppy)
                + seaOctave((uv - vec2(tt * 0.55, -tt * 0.09)) * freq * 1.13, choppy);
        h += d * amp;
        uv = mat2(1.92, 0.55, -0.55, 1.92) * uv;
        freq *= 1.9;
        amp *= 0.22;
        choppy = mix(choppy, 1.0, 0.3);
    }
    return h;
}

vec3 waveNormalWorldFast(vec2 p, float t, float strength, float chopFade) {
    float ang = WIND_DIRECTION * 0.017453293;
    vec2 wx = vec2(cos(ang), sin(ang)); // compile-time constants (WIND_DIRECTION is a define)
    // uv(p) = (dot(p,wx)*1.06, dot(p,(-wx.y,wx.x))*0.78) - sync waveChopLOD
    vec2 uv0 = vec2(dot(p, wx) * 1.06, dot(p, vec2(-wx.y, wx.x)) * 0.78);
    // Tight LOD tiers (see header). chopFade 0.60/0.30/0.08 = ~17/58/288 bl.
    int oct = (chopFade > 0.60) ? 4 : ((chopFade > 0.30) ? 3 : ((chopFade > 0.08) ? 2 : 1));
    float e = mix(0.16, 0.11, clamp(chopFade, 0.0, 1.0)); // sync waveNormalWorld
    // World-space finite-difference offsets mapped through the linear transform
    vec2 dux = vec2(wx.x * 1.06, -wx.y * 0.78) * e; // +x world step in uv space
    vec2 duz = vec2(wx.y * 1.06,  wx.x * 0.78) * e; // +z world step in uv space
    float h0 = waveChopUVFast(uv0, t, oct);
    float hx = waveChopUVFast(uv0 + dux, t, oct);
    float hz = waveChopUVFast(uv0 + duz, t, oct);
    float slopeAtt = mix(0.30, 0.85, clamp(chopFade, 0.0, 1.0)); // sync waveNormalWorld
    vec2 slope = vec2((h0 - hx) / e, (h0 - hz) / e) * strength * slopeAtt;
    // ANALYTIC swell slope, sync waveNormalWorld (WAVES_V2 guaranteed here)
    float swellAtt = mix(0.75, 1.0, clamp(chopFade, 0.0, 1.0));
    slope -= swellSlope(p, t) * strength * swellAtt * 1.25;
    #ifdef WATER_MICRONORMAL
        // High-frequency grain, near-field only (sync waveNormalWorld)
        float mf = smoothstep(0.45, 0.95, chopFade);
        if (mf > 0.001) {
            float md = 0.06;
            vec2 q = p * 3.7 + vec2(t * 0.45, -t * 0.33);
            float n0 = vnoise(q);
            float nx = vnoise(q + vec2(md, 0.0));
            float nz = vnoise(q + vec2(0.0, md));
            slope += vec2((n0 - nx) / md, (n0 - nz) / md) * strength * mf * 0.07;
        }
    #endif
    return normalize(vec3(slope.x, 1.0, slope.y));
}
#endif
