/*{
  "DESCRIPTION": "THE TIDAL CATHEDRAL II — the original procedural fjord refined with human feedback: stable shoreline tracing, weathered mineral surfaces, smaller conifers and richer alpine detail. Procedural geometry, lighting and internal caches; no imported textures. A second study of the original composition.",
  "CREDIT": "OpenAI / Codex",
  "ISFVSN": "2",
  "CATEGORIES": ["Landscape", "Raymarching"],
  "PASSES": [
    { "TARGET": "geologyNoise", "PERSISTENT": true },
    { "TARGET": "forestData", "PERSISTENT": true, "FLOAT": 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 },
    { "NAME": "showBoat", "LABEL": "Sailboat", "TYPE": "bool", "DEFAULT": true }
  ]
}*/

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

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);
}

// Value and analytic gradient from the same lattice. The derivatives remain
// smooth without six additional noise samples or a screen-space approximation.
vec4 noiseGradient(vec3 p) {
    vec3 i=mod(floor(p),16.0), f=fract(p);
    vec3 df=6.0*f*(1.0-f);
    f=f*f*(3.0-2.0*f);
    vec2 at=vec2(mod(i.z,4.0),floor(i.z/4.0))*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 layers=mix(mix(a,b,f.x),mix(c,d,f.x),f.y);
    vec2 dx=mix(b-a,d-c,f.y)*df.x;
    vec2 dy=mix(c-a,d-b,f.x)*df.y;
    return vec4(mix(layers.x,layers.y,f.z),mix(dx.x,dx.y,f.z),mix(dy.x,dy.y,f.z),(layers.y-layers.x)*df.z);
}

float fracturedNoise(vec3 p) {
    vec3 i=mod(floor(p),16.0),f=fract(p);
    vec2 at=vec2(mod(i.z,4.0),floor(i.z/4.0))*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);
}

// Three-dimensional joint planes and their gradient. The two closest stone
// cells define a fissure; unlike painted contour lines it changes illumination.
vec4 stoneJoints(vec3 p) {
    vec3 cell=floor(p), f=fract(p);
    float d1=100.0,d2=100.0;
    vec3 v1=vec3(0),v2=vec3(0);
    for(int i=0;i<27;i++) {
        float fi=float(i);
        vec3 o=vec3(mod(fi,3.0),mod(floor(fi/3.0),3.0),floor(fi/9.0))-1.0;
        vec3 c=cell+o;
        vec3 point=0.18+0.64*vec3(hash13(c),hash13(c+19.7),hash13(c+41.3));
        vec3 v=o+point-f;
        float d=dot(v,v);
        if(d<d1){d2=d1;v2=v1;d1=d;v1=v;}else if(d<d2){d2=d;v2=v;}
    }
    float a=sqrt(max(d1,0.00001)), b=sqrt(max(d2,0.00001));
    return vec4(b-a,v1/a-v2/b);
}

vec4 rockRelief(vec3 p,float footprint) {
    vec3 gradient=vec3(0.0);
    float height=0.0,frequency=1.8,amplitude=0.24;
    for(int i=0;i<4;i++) {
        float filtering=1.0-smoothstep(0.35,1.15,footprint*frequency);
        vec4 v=noiseGradient(p*frequency+vec3(9.7,2.3,5.1)*float(i));
        float softness=sqrt((v.x-0.5)*(v.x-0.5)+0.0025);
        float signum=(v.x-0.5)/softness;
        height+=(0.5-softness)*amplitude*filtering;
        gradient-=v.yzw*signum*frequency*amplitude*filtering;
        frequency*=2.85;amplitude*=0.43;
    }
    return vec4(height,gradient);
}

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;
    if(d>0.5) return (d-0.195)*0.32;
    vec3 chips=vec3(p.x+0.31*p.z,p.y+0.13*p.x,p.z-0.27*p.y);
    d += (fracturedNoise(chips * 2.0) - 0.5) * 0.32;
    d += (fracturedNoise(chips * 5.9+7.1) - 0.5) * 0.07;
    return d * 0.32;
}

vec3 boatSpace(vec3 p) {
    p-=vec3(7.0+0.14*sin(TIME*0.10),0.07+0.009*sin(TIME*1.25),7.0);
    p.xz=mat2(0.94,-0.342,0.342,0.94)*p.xz;
    p.xy=mat2(cos(0.025*sin(TIME*1.1)),-sin(0.025*sin(TIME*1.1)),sin(0.025*sin(TIME*1.1)),cos(0.025*sin(TIME*1.1)))*p.xy;
    return p;
}

float boat(vec3 p) {
    float bound=length(p-vec3(7.0+0.14*sin(TIME*0.10),0.57+0.009*sin(TIME*1.25),7.0))-1.16;
    if(bound>0.4) return bound;
    vec3 q=boatSpace(p);
    float hull=max((length(q/vec3(0.49,0.11,0.145))-1.0)*0.1,q.y-0.035);
    vec3 pole=q-vec3(-0.08,0.57,0.0);
    float mast=max(length(pole.xz)-0.008,abs(pole.y)-0.57);
    vec2 sail=q.xy-vec2(-0.065,0.18);
    float mainSail=max(max(-sail.x,-sail.y),dot(sail,vec2(0.882,0.471))-0.405);
    mainSail=max(mainSail,abs(q.z-0.017-0.035*sin(max(sail.y,0.0)*3.4))-0.009);
    vec2 jib=vec2(-q.x-0.10,q.y-0.15);
    float frontSail=max(max(-jib.x,-jib.y),dot(jib,vec2(0.922,0.388))-0.32);
    frontSail=max(frontSail,abs(q.z+0.015)-0.006);
    return min(min(hull,mast),min(mainSail,frontSail))*0.65;
}

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(showBoat) { float vessel=boat(p);if(vessel<res.x)res=vec2(vessel,4.0); }
    if (p.y - ground < 3.0 && p.y < 18.0 && p.z > -8.0 && p.z < 106.0) {
        vec2 cell = floor(p.xz / FOREST_SPACING);
        vec2 lookup = cell + vec2(64.0,16.0);
        vec4 data = IMG_PIXEL(forestData, lookup + 0.5);
        float seed = data.b;
        vec2 center = (cell + 0.5) * FOREST_SPACING;
        center += vec2(seed - 0.5, hash12(cell + 8.1) - 0.5) * 0.72;
        float ht = data.g;
        if (ht > 0.1 && lookup.x >= 0.0 && lookup.x < 128.0 && lookup.y >= 0.0 && lookup.y < 96.0) {
            vec3 q = p - vec3(center.x, data.r - 0.05, center.y);
            q.x += 0.024 * sin(q.y * 3.0 + seed * 8.0);
            float crown = 0.28 * ht * max(1.0 - q.y / ht, 0.0);
            crown *= 0.80 + 0.16 * sin(q.y * 27.0 + seed * 11.0) + 0.18 * noise3(q * 21.0);
            float foliage = max(length(q.xz) - crown, max(0.13*ht-q.y, q.y - ht));
            float trunk = max(length(q.xz)-0.022*ht,max(-q.y,q.y-ht*0.82));
            float tree = min(foliage,trunk) * 0.48;
            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 < 768; 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 > 383)) break;
    }
    return vec2(min(t, limit), mat);
}

vec3 normalAt(vec3 p, float t) {
    float e = max(0.008, t * 0.45 / max(RENDERSIZE.y,360.0));
    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 geometricNormal = normalAt(p, t);
    vec3 n=geometricNormal;
    vec3 sd = sunDirection();
    if(material>3.5) {
        vec3 q=boatSpace(p);
        vec3 paint=q.y<0.045?vec3(0.035,0.075,0.105):vec3(0.72,0.60,0.40);
        if(abs(q.x+0.08)<0.015)paint=vec3(0.17,0.115,0.065);
        float seams=1.0-0.10*(1.0-smoothstep(0.015,0.05,abs(sin(q.y*24.0+q.x*4.0))));
        return paint*seams*(vec3(0.44,0.49,0.53)+vec3(1.5,1.2,0.80)*abs(dot(n,sd)));
    }
    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;
    }
    float footprint=t/(RENDERSIZE.y*1.08);
    vec4 relief=rockRelief(p,footprint*(reflected?1.8:1.0));
    vec3 jointScale=vec3(0.82,1.75,0.82);
    float jointWidth=max(0.014,footprint*0.25);
    float broken=smoothstep(0.42,0.67,noise3(p*0.55+21.0));
    vec4 joints=vec4(10.0,0.0,0.0,0.0);
    if(broken>0.01 && footprint<0.24)joints=stoneJoints(p*jointScale+vec3(0.4,1.8,2.3));
    float crevice=exp(-joints.x/jointWidth)*broken*(1.0-smoothstep(0.09,0.24,footprint));
    vec3 fractureGradient=0.012/jointWidth*crevice*joints.yzw*jointScale;
    float bedCoordinate=dot(p,vec3(0.27,2.25,0.13));
    float bedding=sin(bedCoordinate*6.2831+noise3(p*0.28)*2.0);
    float bedFilter=1.0-smoothstep(0.05,0.22,footprint);
    vec3 bedGradient=vec3(0.27,2.25,0.13)*0.035*cos(bedCoordinate*6.2831)*bedFilter;
    vec3 gradient=relief.yzw+fractureGradient+bedGradient;
    gradient-=geometricNormal*dot(gradient,geometricNormal);
    float grain=fbm3(p*vec3(1.7,0.22,1.7));
    float smallGrain=noise3(p*8.5);
    float moss=smoothstep(0.38,0.78,geometricNormal.y+(grain-0.45)*0.55);
    moss*=1.0-smoothstep(17.0,26.0,p.y);
    n=normalize(n-gradient*mix(1.15,0.40,moss));
    float iron=smoothstep(0.53,0.75,fbm3(p*0.30+11.0));
    float mineral=clamp((grain-0.43)*2.7+0.46,0.0,1.0);
    vec3 stone=mix(vec3(0.13,0.135,0.13),vec3(0.43,0.39,0.31),mineral);
    stone=mix(stone,vec3(0.31,0.19,0.10),iron*0.55);
    stone*=0.9+0.18*smallGrain+0.025*bedding;
    stone*=1.0-0.19*crevice*(1.0-moss);
    vec3 vegetation=mix(vec3(0.026,0.058,0.012),vec3(0.17,0.225,0.040),grain);
    vegetation*=0.78+0.42*smallGrain;
    stone=mix(stone,vegetation,moss);
    float snowLine = p.y - 25.0 + 5.0 * fbm2(p.xz * 0.16);
    float snow = smoothstep(0.0, 4.0, snowLine) * smoothstep(0.35, 0.82, geometricNormal.y + (grain - 0.5) * 0.35);
    if (material > 1.5) snow = 0.0;
    stone = mix(stone, vec3(0.62, 0.72, 0.79)*(0.92+0.08*smallGrain), snow);
    n=normalize(mix(n,geometricNormal,snow*0.70));
    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 + geometricNormal * s).x / 0.4) / exp2(float(i));
        }
        ao = clamp(1.0 - occ * 0.45, 0.3, 1.0);
    }
    ao*=mix(1.0-0.4*crevice,1.0,max(snow,moss));
    float sh = shadowAt(p + geometricNormal * 0.09, sd);
    float diffuse = max(dot(n, sd), 0.0) * sh * (1.0-0.22*crevice*(1.0-snow));
    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(18.0, 85.0, wet));
    col += vec3(1.1, 0.9, 0.64) * spec * sh * (0.009 + 0.4 * wet + 0.08 * 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);
    if(showBoat) {
        vec3 q=boatSpace(p);
        float behind=-q.x;
        float spread=0.12+0.12*behind;
        float wake=exp(-abs(abs(q.z)-spread)*65.0)*smoothstep(0.3,0.65,behind)*(1.0-smoothstep(1.0,2.7,behind));
        col=mix(col,vec3(0.37,0.57,0.54),wake*0.28);
    }
    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;
    }
    if (PASSINDEX == 1) {
        vec4 prior = IMG_THIS_PIXEL(forestData);
        if(prior.a > 0.5) { gl_FragColor=prior; return; }
        vec2 pixel=floor(gl_FragCoord.xy);
        if(pixel.x >= 128.0 || pixel.y >= 96.0) { gl_FragColor=vec4(0,0,0,1); return; }
        vec2 cell=pixel-vec2(64.0,16.0);
        float seed=hash12(cell+13.5);
        vec2 center=(cell+0.5)*FOREST_SPACING + vec2(seed-0.5,hash12(cell+8.1)-0.5)*0.72;
        float base=terrain(center);
        vec2 slope=vec2(terrain(center+vec2(0.4,0))-base,terrain(center+vec2(0,0.4))-base)/0.4;
        float grove=noise2(center*0.12);
        float ht=0.65+0.75*seed;
        float fertile=smoothstep(0.7,1.8,base)*(1.0-smoothstep(8.0,12.0,base));
        fertile*=1.0-smoothstep(0.7,1.6,length(slope));
        if(center.y < -5.0 || center.y > 102.0 || seed > fertile*0.95 || grove < 0.3) ht=0.0;
        gl_FragColor=vec4(base,ht,seed,1);
        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);
}