Procedural Water with Caustics & Sky

const float EPSILON = 0.001;
const float WATER_HEIGHT = 0.5;
const int MAX_STEPS = 256;
const float MAX_DIST = 50.0;
const float MIN_DIST = 0.0001;
const vec3 LIGHT_POS = vec3(0.0, 0.15, 20.0);
const vec3 SUN_COLOR = vec3(1.0, 0.6, 0.3);
const vec3 SKY_COLOR = vec3(0.2, 0.4, 0.8);
const vec3 HORIZON_COLOR = vec3(0.9, 0.5, 0.3);
const float WATER_REFLECTION_DISTANCE = 12.0;

float hash21(vec2 p) {
    p = fract(p * vec2(233.34, 851.73));
    p += dot(p, p + 23.45);
    return fract(p.x * p.y);
}

float noise(vec2 p) {
    vec2 i = floor(p);
    vec2 f = fract(p);
    float a = hash21(i);
    float b = hash21(i + vec2(1.0, 0.0));
    float c = hash21(i + vec2(0.0, 1.0));
    float d = hash21(i + vec2(1.0, 1.0));
    vec2 u = f * f * (3.0 - 2.0 * f);
    return mix(a, b, u.x) + (c - a) * u.y * (1.0 - u.x) + (d - b) * u.x * u.y;
}

float fbm(vec2 p) {
    float value = 0.0;
    float amplitude = 0.5;
    float frequency = 1.0;
    float sum = 0.0;
    mat2 rotation = mat2(cos(0.5), sin(0.5), -sin(0.5), cos(0.5));
    for(int i = 0; i < 8; i++) {
        value += amplitude * noise(p * frequency);
        p = rotation * p * 1.1;
        sum += amplitude;
        amplitude *= 0.5;
        frequency *= 2.1;
    }
    return value / sum;
}

float waterHeight(vec2 pos, float time) {
    vec2 movement = vec2(time * 0.2, time * 0.1);
    float dist = length(pos);
    float largeWaveAtten = 1.0 / (1.0 + dist * 0.1);
    float smallWaveAtten = 1.0 / (1.0 + dist * 0.3);
    float largeWaves = fbm(pos * 1.0 + movement) * largeWaveAtten;
    float smallWaves = fbm(pos * 3.0 - movement * 1.3) * smallWaveAtten;
    float waves = largeWaves * 0.15 + smallWaves * 0.05;
    waves += fbm(pos * 8.0 - movement * 0.7) * smallWaveAtten * 0.02;
    float minDist = WATER_REFLECTION_DISTANCE * 0.8;
    float maxDist = WATER_REFLECTION_DISTANCE * 1.2;
    float distFade = 1.0 - smoothstep(minDist, maxDist, dist);
    waves *= 0.8 + 0.2 * smoothstep(-1.0, 1.0, waves);
    return waves * distFade;
}

vec3 waterNormal(vec2 pos, float time) {
    float dist = length(pos);
    float baseEps = 0.05;
    float distFactor = 0.01;
    float adaptiveEps = baseEps + dist * distFactor;
    vec2 epsX = vec2(adaptiveEps, 0.0);
    vec2 epsY = vec2(0.0, adaptiveEps);
    float r1 = waterHeight(pos + epsX, time);
    float r2 = waterHeight(pos + epsX * 0.5, time);
    float l1 = waterHeight(pos - epsX, time);
    float l2 = waterHeight(pos - epsX * 0.5, time);
    float t1 = waterHeight(pos + epsY, time);
    float t2 = waterHeight(pos + epsY * 0.5, time);
    float b1 = waterHeight(pos - epsY, time);
    float b2 = waterHeight(pos - epsY * 0.5, time);
    float dX = ((r1 - l1) * 0.3 + (r2 - l2) * 0.7) / (2.0 * adaptiveEps);
    float dY = ((t1 - b1) * 0.3 + (t2 - b2) * 0.7) / (2.0 * adaptiveEps);
    float normalStrength = 1.0 / (1.0 + dist * 0.1);
    return normalize(vec3(-dX * normalStrength, 1.0, -dY * normalStrength));
}

float calculateCaustics(vec3 pos, float time) {
    vec2 waterPos = pos.xz;
    vec3 normal = waterNormal(waterPos, time);
    float ior = 1.33;
    vec3 lightDir = normalize(LIGHT_POS - pos);
    vec3 refracted = refract(lightDir, normal, 1.0 / ior);
    float depth = (pos.y - (-1.0)) / refracted.y;
    vec2 causticPos = pos.xz + refracted.xz * depth;
    float causticIntensity = fbm(causticPos * 5.0 + time * 0.1);
    return pow(max(causticIntensity, 0.0), 2.0);
}

float sceneSDF(vec3 p) {
    float ground = p.y + 1.0;
    float water = p.y - WATER_HEIGHT - waterHeight(p.xz, iTime);
    return min(ground, water);
}

vec2 rayMarch(vec3 ro, vec3 rd) {
    float depth = 0.0;
    float id = 0.0;
    for(int i = 0; i < MAX_STEPS; i++) {
        vec3 p = ro + rd * depth;
        float dist = sceneSDF(p);
        float minDist = MIN_DIST * (1.0 + depth * 0.1);
        if(dist < minDist) {
            if(p.y > -0.99) id = 1.0;
            break;
        }
        if(depth > MAX_DIST) break;
        float stepSize = max(dist * 0.5, minDist);
        depth += stepSize;
    }
    return vec2(depth, id);
}

float rayleighPhase(float cosTheta) {
    return 3.0 * (1.0 + cosTheta * cosTheta) / (16.0 * 3.14159);
}

float miePhase(float cosTheta) {
    const float g = 0.76;
    float g2 = g * g;
    return (1.0 - g2) / pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5);
}

vec3 getSkyColor(vec3 rd) {
    vec3 sunDir = normalize(LIGHT_POS);
    float cosTheta = dot(rd, sunDir);
    float sunDot = max(cosTheta, 0.0);
    float rayleigh = rayleighPhase(cosTheta) * 0.15;
    float mie = miePhase(cosTheta) * 0.05;
    float sunIntensity = 15.0;
    vec3 sun = SUN_COLOR * pow(sunDot, 1500.0) * sunIntensity;
    vec3 corona = vec3(pow(sunDot, 6.0), pow(sunDot, 8.0), pow(sunDot, 10.0)) * 0.8;
    vec3 scattering = mix(vec3(rayleigh) * vec3(0.5, 0.7, 1.0), vec3(mie) * SUN_COLOR, 0.2);
    float horizon = pow(1.0 - abs(rd.y), 3.0);
    vec3 sky = mix(SKY_COLOR, HORIZON_COLOR, horizon);
    return sky + sun + corona + scattering;
}

vec3 toneMap(vec3 color) {
    color *= 0.8;
    float luminance = dot(color, vec3(0.2126, 0.7152, 0.0722));
    float mappedLuminance = luminance * (1.0 + luminance / (5.0 * 5.0)) / (1.0 + luminance);
    color = color * (mappedLuminance / (luminance + 0.001));
    float saturationAdjust = 1.0 / (1.0 + luminance * 0.3);
    color = mix(vec3(luminance), color, saturationAdjust);
    color = mix(color, vec3(luminance), pow(luminance, 2.0) * 0.2);
    return color;
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    vec2 uv = (fragCoord - 0.5 * iResolution.xy) / iResolution.y;
    vec3 ro = vec3(0.0, 2.5, -6.0);
    vec3 lookAt = vec3(0.0, 0.4, 10.0);
    vec3 forward = normalize(lookAt - ro);
    vec3 right = normalize(cross(forward, vec3(0.0, 1.0, 0.0)));
    vec3 up = cross(right, forward);
    float fov = 1.0;
    vec3 rd = normalize(forward + right * uv.x * fov + up * uv.y * fov);
    vec2 march = rayMarch(ro, rd);
    vec3 pos = ro + rd * march.x;
    vec3 col;
    if(march.x < MAX_DIST) {
        if(march.y > 0.5) {
            vec3 normal = waterNormal(pos.xz, iTime);
            vec3 reflectDir = reflect(rd, normal);
            float distanceToCamera = length(pos - ro);
            float baseFresnel = pow(1.0 - max(dot(normal, -rd), 0.0), 5.0);
            float distanceFresnel = smoothstep(0.0, WATER_REFLECTION_DISTANCE, distanceToCamera);
            float fresnel = mix(baseFresnel, 1.0, distanceFresnel * 0.8);
            float sunReflection = pow(max(dot(reflectDir, normalize(LIGHT_POS)), 0.0), 128.0);
            sunReflection *= 1.0 + distanceFresnel * 2.0;
            vec3 specular = SUN_COLOR * sunReflection * 3.0;
            vec3 closeWaterColor = vec3(0.1, 0.2, 0.3);
            vec3 farWaterColor = mix(HORIZON_COLOR, SUN_COLOR, 0.3);
            vec3 baseWaterColor = mix(closeWaterColor, farWaterColor, distanceFresnel);
            vec3 waterColor = mix(baseWaterColor, getSkyColor(reflectDir), fresnel);
            col = waterColor + specular;
        } else {
            float caustics = calculateCaustics(pos, iTime);
            float sunShadow = max(dot(normalize(pos - LIGHT_POS), vec3(0.0, 1.0, 0.0)), 0.0);
            col = vec3(0.2) + caustics * vec3(0.4, 0.6, 0.8) + sunShadow * SUN_COLOR * 0.2;
        }
        float fogAmount = 1.0 - exp(-march.x * 0.05);
        col = mix(col, getSkyColor(rd), fogAmount);
    } else {
        col = getSkyColor(rd);
    }
    col = toneMap(col);
    col = pow(col, vec3(0.4545));
    col = col * (1.0 + col) / (1.0 + col * col);
    fragColor = vec4(col, 1.0);
}