/*{
  "DESCRIPTION": "THE TIDAL CATHEDRAL — an original procedural fjord at first light. A monumental sea-carved arch, glacial peaks, conifer groves, a distant cascade, rippled reflections, wind-driven clouds and circling birds. All geometry, materials and light are synthesized here. The first pass generates a reusable noise lattice; no imported textures or external assets.",
  "CREDIT": "OpenAI / Codex",
  "ISFVSN": "2",
  "CATEGORIES": ["Landscape", "Raymarching"],
  "PASSES": [
    { "TARGET": "geologyNoise", "PERSISTENT": true },
    {}
  ],
  "INPUTS": [
    { "NAME": "cameraDrift", "LABEL": "Camera drift", "TYPE": "float", "DEFAULT": 0.65, "MIN": 0.0, "MAX": 1.0 },
    { "NAME": "sunHeight", "LABEL": "Sun elevation", "TYPE": "float", "DEFAULT": 0.27, "MIN": 0.12, "MAX": 0.85 },
    { "NAME": "atmosphere", "LABEL": "Valley mist", "TYPE": "float", "DEFAULT": 0.5, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "waterMotion", "LABEL": "Water motion", "TYPE": "float", "DEFAULT": 0.6, "MIN": 0.0, "MAX": 1.5 },
    { "NAME": "exposure", "LABEL": "Exposure", "TYPE": "float", "DEFAULT": 0.95, "MIN": 0.6, "MAX": 1.6 }
  ]
}*/

#define FAR 290.0
const mat2 TURN = mat2(0.80, -0.60, 0.60, 0.80);

float hash12(vec2 p) {
    vec3 q = fract(vec3(p.xyx) * 0.1031);
    q += dot(q, q.yzx + 33.33);
    return fract((q.x + q.y) * q.z);
}

float hash13(vec3 p) {
    p = fract(p * 0.1031);
    p += dot(p, p.zyx + 31.32);
    return fract((p.x + p.y) * p.z);
}

float noise2(vec2 p) {
    vec2 i = floor(p), f = fract(p);
    f = f * f * (3.0 - 2.0 * f);
    vec2 at = mod(i + 32.0, 64.0) + 0.5;
    float a = IMG_PIXEL(geologyNoise, at).g;
    float b = IMG_PIXEL(geologyNoise, at + vec2(1,0)).g;
    float c = IMG_PIXEL(geologyNoise, at + vec2(0,1)).g;
    float d = IMG_PIXEL(geologyNoise, at + vec2(1,1)).g;
    return mix(mix(a,b,f.x),mix(c,d,f.x),f.y);
}

float noise3(vec3 p) {
    vec3 i = floor(p), f = fract(p);
    f = f * f * (3.0 - 2.0 * f);
    i = mod(i, 16.0);
    vec2 tile0 = vec2(mod(i.z, 4.0), floor(i.z / 4.0));
    vec2 at = tile0 * 18.0 + i.xy + 1.5;
    vec2 a = IMG_PIXEL(geologyNoise, at).rb;
    vec2 b = IMG_PIXEL(geologyNoise, at + vec2(1,0)).rb;
    vec2 c = IMG_PIXEL(geologyNoise, at + vec2(0,1)).rb;
    vec2 d = IMG_PIXEL(geologyNoise, at + vec2(1,1)).rb;
    vec2 pair = mix(mix(a,b,f.x),mix(c,d,f.x),f.y);
    return mix(pair.x, pair.y, f.z);
}

float fbm2(vec2 p) {
    float v = 0.0, a = 0.5;
    for (int i = 0; i < 5; i++) {
        v += a * noise2(p);
        p = TURN * p * 2.07 + vec2(17.4, 9.2);
        a *= 0.48;
    }
    return v;
}

float fbm3(vec3 p) {
    float v = 0.0, a = 0.5;
    for (int i = 0; i < 4; i++) {
        v += a * noise3(p);
        p = p.yzx * 2.03 + vec3(11.7, 4.3, 8.1);
        a *= 0.48;
    }
    return v;
}

float ridge(vec2 p) {
    float v = 0.0, a = 0.55, weight = 1.0;
    for (int i = 0; i < 4; i++) {
        float r = 1.0 - abs(2.0 * noise2(p) - 1.0);
        r *= r;
        v += a * r * weight;
        weight = clamp(r * 1.8, 0.15, 1.0);
        p = TURN * p * 2.11 + 3.7;
        a *= 0.49;
    }
    return v;
}

float riverPath(float z) {
    return 4.0 + 5.0 * sin(z * 0.045) - 12.0 * smoothstep(42.0, 80.0, z);
}

float terrain(vec2 p) {
    float channel = riverPath(p.y);
    float bank = smoothstep(9.0, 25.0, abs(p.x - channel));
    float back = smoothstep(42.0, 80.0, p.y);
    float r = ridge(p * 0.035 + vec2(2.8, 7.3));
    float h = -2.7 + bank * (3.8 + 12.0 * fbm2(p * 0.085));
    float peaks = max(0.0, 1.0 - length((p - vec2(-39.0, 105.0)) * vec2(0.020, 0.019))) * 48.0;
    peaks = max(peaks, max(0.0, 1.0 - length((p - vec2(43.0, 123.0)) * vec2(0.021, 0.016))) * 57.0);
    peaks = max(peaks, max(0.0, 1.0 - length((p - vec2(-5.0, 178.0)) * vec2(0.014, 0.015))) * 72.0);
    peaks = max(peaks, max(0.0, 1.0 - length((p - vec2(-89.0, 152.0)) * vec2(0.020, 0.016))) * 61.0);
    h += back * (4.0 + peaks * (0.72 + 0.48 * r) + 5.0 * r);
    h *= 1.0 - 0.55 * exp(-abs(p.x - channel) * 0.3) * back;
    vec2 nearRock = (p - vec2(-11.0, -13.0)) / vec2(9.5, 11.0);
    h += 7.5 * exp(-dot(nearRock, nearRock) * 1.5);
    h += 0.12 * noise2(p * 2.8);
    return h;
}

// Elliptical stone vault, leaned and twisted once in world space.
// The enclosing box avoids evaluating erosion noise outside the arch's region.
float arch(vec3 p) {
    vec3 q = p - vec3(-7.0, 2.0, 27.0);
    vec3 b = abs(q) - vec3(21.0, 24.0, 9.0);
    float bound = length(max(b, 0.0));
    if (bound > 2.0) return bound + 0.2;
    q.x += q.y * 0.15 + sin(q.y * 0.24) * 0.65;
    q.z += q.x * 0.15;
    q.y += sin(q.x * 0.19) * 1.1;
    float thick = 2.15 + 1.1 * (1.0 - smoothstep(-12.0, 0.0, q.x)) + 0.65 * noise3(p * 0.17);
    float d = length(vec2(length(q.xy * vec2(1.0, 0.76)) - 12.0, q.z * 1.3)) - thick;
    d += (noise3(p * 0.46) - 0.5) * 1.65;
    d += (noise3(p * 1.7) - 0.5) * 0.28;
    return d * 0.40;
}

vec2 scene(vec3 p) {
    float ground = terrain(p.xz);
    float land = (p.y - ground) * 0.23;
    float vault = arch(p);
    vec2 res = land < vault ? vec2(land, 1.0) : vec2(vault, 2.0);
    if (p.y - ground < 4.2 && p.y < 21.0 && p.z > -5.0 && p.z < 94.0) {
        vec2 cell = floor(p.xz / 2.2);
        float seed = hash12(cell + 13.5);
        vec2 center = (cell + 0.5) * 2.2;
        center += vec2(seed - 0.5, hash12(cell + 8.1) - 0.5) * 0.6;
        float base = terrain(center);
        if (base > 1.2 && base < 10.5 && seed > 0.32 && noise2(center * 0.15) > 0.37) {
            float ht = 0.85 + seed * 1.4;
            vec3 q = p - vec3(center.x, base - 0.2, center.y);
            q.x += 0.04 * sin(q.y * 3.0 + seed * 8.0);
            float radius = 0.59 * max(1.0 - q.y / ht, 0.0);
            radius *= 0.84 + 0.12 * sin(q.y * 16.0 + seed * 8.0) + 0.12 * noise3(q * 13.0);
            float tree = max(length(q.xz) - radius, max(-q.y, q.y - ht)) * 0.45;
            if (tree < res.x) res = vec2(tree, 3.0);
        }
    }
    return res;
}

vec2 trace(vec3 ro, vec3 rd, float limit, bool reflected) {
    float t = 0.06;
    float mat = -1.0;
    for (int i = 0; i < 220; i++) {
        vec2 h = scene(ro + rd * t);
        float eps = max(0.004, t * (reflected ? 0.0014 : 0.00055));
        if (h.x < eps) { mat = h.y; break; }
        t += max(h.x * 0.91, 0.008);
        if (t > limit || (reflected && i > 145)) break;
    }
    return vec2(min(t, limit), mat);
}

vec3 normalAt(vec3 p, float t) {
    float e = max(0.008, t * 0.00035);
    vec3 n = vec3(0.0);
    for (int i = 0; i < 4; i++) {
        vec3 k = i == 0 ? vec3(1,-1,-1) : (i == 1 ? vec3(-1,-1,1) : (i == 2 ? vec3(-1,1,-1) : vec3(1,1,1)));
        n += k * scene(p + k * e).x;
    }
    return normalize(n);
}

vec3 sunDirection() { return normalize(vec3(-0.65, sunHeight, -0.50)); }

float shadowAt(vec3 p, vec3 lightDir) {
    float sh = 1.0, t = 0.13;
    for (int i = 0; i < 36; i++) {
        float h = scene(p + lightDir * t).x;
        sh = min(sh, 12.0 * max(h, 0.0) / t);
        t += clamp(h, 0.12, 4.0);
        if (sh < 0.015 || t > 65.0) break;
    }
    return clamp(sh, 0.0, 1.0);
}

vec3 sky(vec3 rd, bool clouds) {
    vec3 sd = sunDirection();
    float y = max(rd.y, 0.0);
    float sun = max(dot(rd, sd), 0.0);
    vec3 horizon = mix(vec3(0.48,0.60,0.67), vec3(0.82,0.56,0.34), smoothstep(0.15,-0.8,rd.x) * 0.8);
    vec3 col = mix(horizon, vec3(0.018, 0.085, 0.19), pow(y, 0.42));
    col += vec3(1.0, 0.49, 0.19) * pow(sun, 8.0) * 0.32;
    col += vec3(1.0, 0.65, 0.28) * pow(sun, 70.0) * 0.4;
    col += vec3(9.0, 6.4, 3.6) * smoothstep(0.9995, 0.99985, sun);
    if (clouds && rd.y > 0.015) {
        float trans = 1.0;
        vec3 cloudCol = vec3(0.0);
        float start = 70.0 / rd.y;
        float stride = 22.0 / (rd.y * 9.0);
        for (int i = 0; i < 9; i++) {
            vec3 p = rd * (start + (float(i) + 0.5) * stride);
            p.x += TIME * 0.7;
            vec3 cp = p * vec3(0.010, 0.035, 0.010);
            float shape = fbm3(cp * 1.5 + vec3(1.1, 0, 3.0));
            float edge = smoothstep(70.0, 76.0, p.y) * (1.0 - smoothstep(84.0, 92.0, p.y));
            float den = smoothstep(0.46, 0.66, shape) * edge;
            float a = (1.0 - exp(-den * stride * 0.16)) * smoothstep(0.03, 0.19, rd.y);
            float rim = pow(sun, 12.0);
            vec3 lighting = mix(vec3(0.25, 0.33, 0.42), vec3(0.95, 0.87, 0.76), clamp(0.45 + (0.59 - shape) * 3.0 + rim * 0.5, 0.0, 1.0));
            cloudCol += trans * a * lighting;
            trans *= 1.0 - a;
        }
        col = col * trans + cloudCol;
        // High, stretched cirrus adds a second scale to the sky.
        vec2 cp = rd.xz / max(rd.y, 0.05) * vec2(0.9, 4.5) + vec2(TIME * 0.006, 0);
        float veil = smoothstep(0.55, 0.76, fbm2(cp));
        col = mix(col, vec3(0.91, 0.9, 0.85), veil * 0.22);
    }
    return col;
}

vec3 haze(vec3 col, vec3 ro, vec3 rd, float t) {
    float dy = rd.y * 0.10;
    float integral = abs(dy) < 0.0001 ? t : (1.0 - exp(-dy * t)) / dy;
    float mist = integral * exp(-max(ro.y, 0.0) * 0.10);
    float optical = t * 0.0014 + atmosphere * 0.006 * max(mist, 0.0);
    float lit = pow(max(dot(rd, sunDirection()), 0.0), 6.0);
    vec3 air = mix(vec3(0.32, 0.50, 0.59), vec3(1.1, 0.82, 0.50), lit * 0.72);
    return mix(col, air, 1.0 - exp(-optical));
}

vec3 shadeRock(vec3 p, vec3 rd, float t, float material, bool reflected) {
    vec3 n = normalAt(p, t);
    vec3 sd = sunDirection();
    if (material > 2.5) {
        float sh = shadowAt(p + n * 0.12, sd);
        float needle = noise3(p * 18.0);
        vec3 green = mix(vec3(0.018, 0.038, 0.021), vec3(0.080, 0.115, 0.041), needle);
        vec3 light = vec3(0.32, 0.47, 0.58) + vec3(1.5, 1.05, 0.54) * max(dot(n, sd), 0.0) * sh;
        light += vec3(0.28, 0.36, 0.12) * pow(max(dot(rd, sd), 0.0), 3.0);
        return green * light;
    }
    // Microrelief is tangential, preserving the vault's geometric silhouette.
    float be = 0.035;
    vec3 bump = vec3(noise3((p + vec3(be,0,0)) * 7.0) - noise3((p - vec3(be,0,0)) * 7.0),
                     noise3((p + vec3(0,be,0)) * 7.0) - noise3((p - vec3(0,be,0)) * 7.0),
                     noise3((p + vec3(0,0,be)) * 7.0) - noise3((p - vec3(0,0,be)) * 7.0));
    n = normalize(n - (bump - n * dot(bump, n)) * 0.65 / (1.0 + t * 0.009));
    float grain = fbm3(p * vec3(1.3, 1.8, 1.3));
    float layer = sin(p.y * 5.2 + noise3(p * 0.42) * 8.0 + p.x * 0.3);
    float fine = noise3(p * 13.0);
    vec3 stone = mix(vec3(0.10, 0.095, 0.08), vec3(0.46, 0.33, 0.21), grain);
    stone *= 0.83 + 0.26 * smoothstep(-0.85, 0.15, layer) + 0.15 * fine;
    float cracks = smoothstep(0.025, 0.12, abs(sin(p.y * 2.9 + noise3(p * 0.58) * 5.0)));
    stone *= 0.72 + 0.28 * cracks;
    float moss = smoothstep(0.38, 0.82, n.y + (grain - 0.5) * 0.65);
    moss *= 1.0 - smoothstep(19.0, 29.0, p.y);
    vec3 vegetation = mix(vec3(0.035, 0.07, 0.014), vec3(0.18, 0.24, 0.045), grain);
    stone = mix(stone, vegetation, moss * 0.95);
    float snowLine = p.y - 25.0 + 5.0 * fbm2(p.xz * 0.16);
    float snow = smoothstep(0.0, 4.0, snowLine) * smoothstep(0.38, 0.83, n.y + (grain - 0.5) * 0.35);
    if (material > 1.5) snow = 0.0;
    stone = mix(stone, vec3(0.62, 0.72, 0.79), snow);
    float wet = 1.0 - smoothstep(0.15, 1.7, p.y);
    stone *= 1.0 - wet * 0.5;
    float ao = 1.0;
    if (!reflected) {
        float occ = 0.0;
        for (int i = 1; i <= 3; i++) {
            float s = float(i) * 0.65;
            occ += max(0.0, s - scene(p + n * s).x / 0.4) / exp2(float(i));
        }
        ao = clamp(1.0 - occ * 0.45, 0.3, 1.0);
    }
    float sh = shadowAt(p + n * 0.09, sd);
    float diffuse = max(dot(n, sd), 0.0) * sh;
    vec3 ambient = vec3(0.25, 0.38, 0.53) * (0.38 + 0.62 * max(n.y, 0.0)) * ao;
    ambient += vec3(0.13, 0.17, 0.12) * max(0.0, -n.y + 0.25);
    vec3 col = stone * (ambient + vec3(2.2, 1.30, 0.65) * diffuse);
    float stream = (1.0 - smoothstep(0.28, 0.72, abs(p.x - riverPath(p.z)))) * smoothstep(52.0, 58.0, p.z) * (1.0 - smoothstep(88.0, 100.0, p.z));
    if (material < 1.5 && p.y < 22.0) {
        float streak = noise2(vec2(p.x * 15.0, p.y * 2.0 + TIME * 2.5));
        vec3 cascade = mix(vec3(0.17, 0.38, 0.42), vec3(0.68, 0.79, 0.77), streak * 0.8 + 0.2);
        col = mix(col, cascade * (0.65 + 0.35 * sh), stream);
    }
    vec3 halfDir = normalize(sd - rd);
    float spec = pow(max(dot(n, halfDir), 0.0), mix(28.0, 85.0, wet));
    col += vec3(1.1, 0.9, 0.64) * spec * sh * (0.04 + 0.4 * wet + 0.12 * snow);
    return col;
}

float wave(vec2 p) {
    float time = TIME * waterMotion;
    float v = sin(p.x * 1.9 + p.y * 0.7 - time * 1.4) * 0.35;
    v += sin(p.x * -0.8 + p.y * 2.6 + time * 1.7) * 0.20;
    v += sin(p.x * 4.8 + p.y * 3.1 - time * 2.1) * 0.04;
    v += noise2(p * 6.0 + vec2(time * 0.5, 0.0)) * 0.08;
    return v;
}

vec3 waterNormal(vec3 p, float t) {
    float scale = mix(0.025, 0.006, smoothstep(10.0, 110.0, t));
    float w = wave(p.xz);
    return normalize(vec3((w - wave(p.xz + vec2(0.06, 0.0))) * scale / 0.06,
                           1.0, (w - wave(p.xz + vec2(0.0, 0.06))) * scale / 0.06));
}

vec3 water(vec3 p, vec3 ro, vec3 rd, float t, vec3 n, vec3 reflection) {
    float w = wave(p.xz);
    float depth = max(0.0, -terrain(p.xz));
    vec3 submerged = mix(vec3(0.12, 0.23, 0.16), vec3(0.007, 0.085, 0.095), 1.0 - exp(-depth * 0.65));
    float caustic = pow(max(0.0, 0.5 + 0.5 * sin(p.x * 3.5 + w * 4.0) * sin(p.z * 4.1 - w * 5.0)), 5.0);
    submerged += vec3(0.10, 0.18, 0.12) * caustic * exp(-depth * 0.8);
    float fresnel = 0.035 + 0.965 * pow(1.0 - max(dot(-rd, n), 0.0), 5.0);
    vec3 col = mix(submerged, reflection, fresnel);
    vec3 sd = sunDirection();
    float glitter = pow(max(dot(reflect(-sd, n), -rd), 0.0), 350.0);
    col += vec3(5.0, 3.7, 1.9) * glitter;
    float shore = min(depth, max(arch(p) * 2.0, 0.0));
    float foam = (1.0 - smoothstep(0.08, 0.42, shore)) * smoothstep(0.2, 0.68, noise2(p.xz * 9.0 + w));
    col = mix(col, vec3(0.55, 0.68, 0.65), foam * 0.6);
    return haze(col, ro, rd, t);
}

vec3 tonemap(vec3 x) {
    x *= exposure;
    return clamp((x * (2.51 * x + 0.03)) / (x * (2.43 * x + 0.59) + 0.14), 0.0, 1.0);
}

float birds(vec3 ro, vec3 rd, float distanceLimit) {
    float flock = 0.0;
    for(int i=0;i<5;i++) {
        float fi = float(i);
        vec3 center = vec3(8.0 + fi * 1.6 + 3.0 * sin(TIME * 0.12),
                           25.0 + sin(fi * 1.9) * 1.3 + sin(TIME * 0.22) * 0.6,
                           54.0 + fi * 2.0);
        float t = (center.z - ro.z) / max(rd.z,0.01);
        if(t > 0.0 && t < distanceLimit) {
            vec2 q = (ro + rd * t - center).xy;
            float wing = sin(TIME * 3.5 + fi * 1.7) * 0.38;
            float curve = abs(q.y - abs(q.x) * wing - 0.06 * pow(abs(q.x) / 0.38,2.0));
            float pixel = t / RENDERSIZE.y * 0.55;
            float silhouette = (1.0-smoothstep(0.015,0.015+pixel,curve)) * (1.0-smoothstep(0.28,0.4,abs(q.x)));
            flock = max(flock,silhouette);
        }
    }
    return flock;
}

void main() {
    if (PASSINDEX == 0) {
        vec4 prior = IMG_THIS_PIXEL(geologyNoise);
        if (prior.a > 0.5) { gl_FragColor = prior; return; }
        // A padded, 16-cubed noise lattice generated entirely by this shader.
        // Explicit interpolation keeps terrain derivatives smooth across GPUs.
        vec2 pixel = floor(gl_FragCoord.xy);
        vec2 tile = floor(pixel / 18.0);
        vec3 lattice = vec3(mod(mod(pixel,18.0) - 1.0,16.0), mod(tile.x + tile.y * 4.0,16.0));
        float groundNoise = hash12(mod(pixel,64.0) - 32.0);
        gl_FragColor = vec4(hash13(lattice),groundNoise,hash13(mod(lattice+vec3(0,0,1),16.0)),1.0);
        return;
    }
    vec2 uv = (gl_FragCoord.xy - 0.5 * RENDERSIZE) / RENDERSIZE.y;
    float clockTime = TIME * 0.075;
    vec3 ro = vec3(1.2 + sin(clockTime) * 1.8 * cameraDrift,
                   5.0 + sin(clockTime * 0.7) * 0.3 * cameraDrift, -23.0);
    vec3 target = vec3(-1.0, 10.3, 35.0);
    vec3 forward = normalize(target - ro);
    vec3 right = normalize(cross(vec3(0,1,0), forward));
    vec3 up = cross(forward, right);
    vec3 rd = normalize(forward * 1.08 + right * uv.x + up * uv.y);
    float waterT = rd.y < -0.001 ? -ro.y / rd.y : FAR;
    vec3 rayOrigin = ro, rayDir = rd;
    vec3 waterPoint = vec3(0.0), waterN = vec3(0,1,0);
    vec2 hit = vec2(FAR, -1.0);
    bool isWater = false;
    float limit = min(waterT, FAR);
    // One tracing call site and one material call site serve both ray segments.
    for (int bounce = 0; bounce < 2; bounce++) {
        hit = trace(rayOrigin, rayDir, limit, bounce == 1);
        if (bounce == 0 && hit.y < 0.0 && waterT < FAR) {
            isWater = true;
            waterPoint = ro + rd * waterT;
            waterN = waterNormal(waterPoint, waterT);
            rayOrigin = waterPoint + vec3(0.0, 0.10, 0.0);
            rayDir = reflect(rd, waterN);
            rayDir.y = max(rayDir.y, 0.015);
            rayDir = normalize(rayDir);
            limit = 230.0;
        } else break;
    }
    vec3 col = vec3(0.0);
    if (hit.y > 0.0) {
        vec3 p = rayOrigin + rayDir * hit.x;
        col = shadeRock(p, rayDir, hit.x, hit.y, isWater);
        col = haze(col, rayOrigin, rayDir, hit.x);
    } else col = sky(rayDir, true);
    if (isWater) col = water(waterPoint, ro, rd, waterT, waterN, col);
    float occlusionDistance = isWater ? waterT : (hit.y > 0.0 ? hit.x : FAR);
    col = mix(col, vec3(0.055,0.073,0.085), birds(ro,rd,occlusionDistance) * 0.8);
    // Restrained optical bloom around the actual sun direction.
    float glow = pow(max(dot(rd, sunDirection()), 0.0), 24.0);
    col += vec3(0.14, 0.073, 0.025) * glow;
    col = pow(tonemap(col), vec3(1.0 / 2.2));
    vec2 screen = isf_FragNormCoord;
    float vignette = pow(max(16.0 * screen.x * screen.y * (1.0-screen.x) * (1.0-screen.y), 0.0), 0.10);
    col *= 0.82 + 0.18 * vignette;
    col += (hash12(gl_FragCoord.xy) - 0.5) / 255.0;
    gl_FragColor = vec4(col, 1.0);
}