/*{
  "DESCRIPTION": "A stone sundial that knows what time it is. The sun is computed from the DATE uniform - your machine's real clock and calendar: solar declination from the day of the year, hour angle from the seconds since local midnight. The gnomon shadow therefore falls on the true hour line, the light swings north and south with the season, and 3am is genuinely 3am. At night the same date drives the moon phase and position (it casts its own shadow) and the rotation of the star field. Set your latitude; nudge clockOffset for longitude-within-zone and DST. Turn useRealTime off to scrub the day and the year by hand.",
  "CREDIT": "meridian",
  "CATEGORIES": ["daylight", "astronomy"],
  "INPUTS": [
    { "NAME": "latitude",    "TYPE": "float", "DEFAULT": 45.0, "MIN": -58.0, "MAX": 66.0 },
    { "NAME": "useRealTime", "TYPE": "bool",  "DEFAULT": true },
    { "NAME": "solarHour",   "TYPE": "float", "DEFAULT": 15.4, "MIN": 0.0,  "MAX": 24.0 },
    { "NAME": "dayOfYear",   "TYPE": "float", "DEFAULT": 172.0,"MIN": 1.0,  "MAX": 365.0 },
    { "NAME": "clockOffset", "TYPE": "float", "DEFAULT": 0.0,  "MIN": -120.0, "MAX": 120.0 },
    { "NAME": "viewAngle",   "TYPE": "float", "DEFAULT": 22.0, "MIN": -75.0, "MAX": 75.0 },
    { "NAME": "stoneWear",   "TYPE": "float", "DEFAULT": 0.55, "MIN": 0.0,  "MAX": 1.0 },
    { "NAME": "starGain",    "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.0,  "MAX": 2.0 },
    { "NAME": "exposure",    "TYPE": "float", "DEFAULT": 1.0,  "MIN": 0.3,  "MAX": 2.5 }
  ]
}*/

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#endif

#define PI   3.14159265359
#define TAU  6.28318530718
#define D2R  0.01745329252

#define PLATE_R 1.05
#define PLATE_T 0.13
#define GN_HX   0.023
#define GROUND_Y (-1.55)
#define FOCAL   0.87

// ---------------------------------------------------------------- globals
float gSP, gCP, gSGN;      // sin/cos/sign of latitude
float gGL, gGH;            // gnomon base length, apex height
vec3  gSun, gMoon;         // unit vectors, world (x=east, y=up, z=south)
float gSunI, gMoonI, gMoonK;
vec3  gSunCol, gMoonCol;
vec3  gZen, gHor;          // sky palette keys
vec2  gSunH;               // sun azimuth on the horizon
float gNight;              // 0 day .. 1 astronomical night
float gStarF;              // star visibility, lags the night factor
float gStarRot;
vec3  gAmb;                // hemispheric skylight
float gWarm, gFrost;       // stone thermal state
float gAdapt;              // scene-uniform eye adaptation
float gGlowI;
vec3  gGlowCol;            // colour of the twilight horizon
vec3  gFill;               // warm directional fill from the glowing horizon
vec3  gFillDir;

// ---------------------------------------------------------------- noise
float hash21(vec2 p){
  p = fract(p*vec2(123.34, 456.21));
  p += dot(p, p + 45.32);
  return fract(p.x*p.y);
}
float vnoise(vec2 p){
  vec2 i = floor(p), f = fract(p);
  f = f*f*(3.0 - 2.0*f);
  float a = hash21(i), b = hash21(i + vec2(1.0,0.0));
  float c = hash21(i + vec2(0.0,1.0)), d = hash21(i + vec2(1.0,1.0));
  return mix(mix(a,b,f.x), mix(c,d,f.x), f.y);
}
float fbm(vec2 p){
  float s = 0.0, a = 0.5;
  for (int i = 0; i < 4; i++){ s += a*vnoise(p); p *= 2.03; a *= 0.5; }
  return s;
}
// gradient of two octaves of value noise - surface relief without dFdx
vec2 bumpGrad(vec2 tc){
  float e = 0.015;
  float a0 = vnoise(tc), a1 = vnoise(tc*3.7 + 5.0);
  vec2 g;
  g.x = (vnoise(tc + vec2(e,0.0)) - a0) + 0.45*(vnoise(tc*3.7 + 5.0 + vec2(e,0.0)) - a1);
  g.y = (vnoise(tc + vec2(0.0,e)) - a0) + 0.45*(vnoise(tc*3.7 + 5.0 + vec2(0.0,e)) - a1);
  return g/e;
}

// ---------------------------------------------------------------- geometry
// offset vector from the closest point on segment a-b to p
vec2 segOff(vec2 p, vec2 a, vec2 b){
  vec2 pa = p - a, ba = b - a;
  float h = clamp(dot(pa,ba)/max(dot(ba,ba), 1e-6), 0.0, 1.0);
  return pa - ba*h;
}

float sdPlate(vec3 p){
  float r = 0.016;
  vec2 d = vec2(length(p.xz) - (PLATE_R - r), abs(p.y + PLATE_T*0.5) - (PLATE_T*0.5 - r));
  return min(max(d.x,d.y), 0.0) + length(max(d, 0.0)) - r;
}
float sdGnomon(vec3 p){
  float zz = p.z*gSGN;                 // triangle occupies zz in [-gGL, 0]
  float t  = gGH/max(gGL, 1e-4);       // tan|phi|
  float c  = inversesqrt(1.0 + t*t);
  float d  = -p.y;
  d = max(d, -(zz + gGL));
  d = max(d, (p.y + zz*t)*c);
  d = max(d, abs(p.x) - GN_HX);
  return d - 0.006;
}
float sdColumn(vec3 p){
  float rr = 0.022;
  vec2 d1 = vec2(length(p.xz) - (0.60 - rr), abs(p.y + 0.215) - (0.085 - rr));
  float collar = min(max(d1.x,d1.y),0.0) + length(max(d1,0.0)) - rr;
  vec2 d2 = vec2(length(p.xz) - 0.40, abs(p.y + 1.02) - 0.74);
  float shaft = min(max(d2.x,d2.y),0.0) + length(max(d2,0.0));
  return min(collar, shaft);
}
float map(vec3 p){
  return min(min(sdPlate(p), sdGnomon(p)), sdColumn(p));
}
vec3 calcNormal(vec3 p){
  vec2 e = vec2(1.0, -1.0)*0.0012;
  return normalize(e.xyy*map(p + e.xyy) + e.yyx*map(p + e.yyx) +
                   e.yxy*map(p + e.yxy) + e.xxx*map(p + e.xxx));
}
// ray vs bounding sphere of the whole monument
bool boundHit(vec3 ro, vec3 rd, out float t0, out float t1){
  vec3 oc = ro - vec3(0.0, -0.55, 0.0);
  float b = dot(oc, rd);
  float c = dot(oc,oc) - 2.45*2.45;
  float h = b*b - c;
  if (h < 0.0) return false;
  h = sqrt(h);
  t0 = max(-b - h, 0.0); t1 = -b + h;
  return t1 > 0.0;
}
float softShadow(vec3 ro, vec3 rd, float k){
  float b0, b1;
  if (!boundHit(ro, rd, b0, b1)) return 1.0;
  float res = 1.0, t = max(b0, 0.012);
  for (int i = 0; i < 44; i++){
    float h = map(ro + rd*t);
    res = min(res, k*h/t);
    t += clamp(h, 0.012, 0.30);
    if (res < 0.004 || t > b1) break;
  }
  return clamp(res, 0.0, 1.0);
}
float ao(vec3 p, vec3 n){
  float o = 0.0, s = 1.0;
  for (int i = 0; i < 4; i++){
    float d = 0.022 + 0.06*float(i);
    o += (d - map(p + n*d))*s;
    s *= 0.72;
  }
  return clamp(1.0 - 2.4*o, 0.0, 1.0);
}

// ---------------------------------------------------------------- sky
vec3 starField(vec3 d){
  vec3 pole = vec3(0.0, abs(gSP), -gSGN*gCP);
  vec3 e1 = vec3(1.0, 0.0, 0.0);
  vec3 e2 = cross(pole, e1);
  float a = dot(d,e1), b = dot(d,e2), c = dot(d,pole);
  float cs = cos(gStarRot), sn = sin(gStarRot);
  vec3 sc = vec3(a*cs - b*sn, a*sn + b*cs, c);

  vec3 q = sc/(abs(sc.x) + abs(sc.y) + abs(sc.z));
  vec2 oc = (q.z >= 0.0) ? q.xy : (vec2(1.0) - abs(q.yx))*sign(q.xy);

  // milky way: a great circle band in the fixed sky frame
  float bandv = dot(normalize(sc), vec3(-0.868, -0.196, 0.456));
  float mw = exp(-bandv*bandv*26.0);
  mw *= 0.30 + 1.10*fbm(oc*7.0 + 11.0);
  mw *= 0.45 + 0.55*smoothstep(0.02, 0.16, abs(bandv));   // dark rift
  vec3 col = vec3(0.52, 0.53, 0.62)*mw*0.055;

  for (int L = 0; L < 2; L++){
    float dens = (L == 0) ? 44.0 : 88.0;
    vec2 g = oc*dens;
    vec2 ic = floor(g), fc = fract(g);
    float h = hash21(ic + float(L)*57.13);
    vec2 ctr = vec2(fract(h*197.31), fract(h*71.77));
    float mag = pow(fract(h*3771.3), 5.5)*((L == 0) ? 1.0 : 0.38);
    vec2 dv = fc - ctr;
    float d2 = dot(dv,dv);
    float tw = 0.78 + 0.22*sin(TIME*(1.4 + 3.1*fract(h*17.3)) + h*TAU);
    vec3 sCol = mix(vec3(0.70,0.79,1.0), vec3(1.0,0.80,0.60), fract(h*13.71));
    col += sCol*mag*tw*(0.0016/(d2 + 0.00022));
  }
  return col*smoothstep(-0.03, 0.20, d.y);
}

vec3 skyColor(vec3 d, bool withBodies){
  float h = clamp(d.y, 0.0, 1.0);
  vec3 col = mix(gHor, gZen, pow(h, 0.44));

  // horizon glow at the sun azimuth (persists while the sun is below)
  vec2 dh = d.xz;
  float ts = dot(normalize(dh + vec2(1e-4)), gSunH);
  float ev = abs(d.y);                       // decay both sides of the horizon
  float glow = pow(max(ts, 0.0), 2.6)*exp(-ev*6.2);
  col += gGlowCol*glow*gGlowI;

  // Facing away from the sunset you see the earth's own shadow rise out of the
  // horizon, capped by the pink Belt of Venus. Its top sits at the antisolar
  // point, i.e. as far above the horizon as the sun is below it.
  float anti = max(-ts, 0.0)*max(-ts, 0.0);
  float st   = clamp(-gSun.y, 0.0, 0.40);
  float inSh = smoothstep(st, st - 0.075, d.y)*step(0.0, d.y);
  float bu   = (d.y - st)/0.070;
  float belt = exp(-bu*bu);
  float tw   = smoothstep(0.030, -0.010, gSun.y)*smoothstep(-0.330, -0.150, gSun.y);
  col = mix(col, col*vec3(0.66, 0.74, 0.98), inSh*anti*tw*0.85);
  col += vec3(0.34, 0.135, 0.165)*belt*anti*tw*1.30;

  col += starField(d)*gStarF*starGain;

  // high cirrus, drifting; lit from below when the sun is low
  if (d.y > 0.002){
    vec2 cp = d.xz/(d.y + 0.055)*0.42 + vec2(TIME*0.0035, TIME*0.0012);
    float cv = fbm(vec2(cp.x*0.42, cp.y*1.9) + 6.0);
    cv = smoothstep(0.47, 0.80, cv)*smoothstep(0.004, 0.075, d.y);
    vec3 cCol = gHor*0.85 + gGlowCol*gGlowI*1.15*pow(max(ts, 0.0), 1.8)
              + gSunCol*gSunI*0.052;
    col = mix(col, cCol, cv*0.80);
  }

  if (withBodies){
    float sdot = dot(d, gSun);
    float sunUp = smoothstep(-0.030, 0.004, gSun.y);
    col += gSunCol*(pow(max(sdot,0.0), 2200.0)*22.0 + pow(max(sdot,0.0), 26.0)*0.45)*sunUp;
    col += vec3(1.0,0.94,0.85)*smoothstep(0.99990, 0.99996, sdot)*26.0*sunUp;

    // moon: real limb, real terminator (its normal dotted with the sun)
    float md = dot(d, gMoon);
    float moonUp = smoothstep(-0.020, 0.020, gMoon.y);
    float rho = 0.0125;
    if (md > 0.9995 && moonUp > 0.0){
      vec3 e = (d - gMoon*md)/rho;
      float r = min(length(e), 1.0);
      vec3 nm = normalize(gMoon*sqrt(max(1.0 - r*r, 0.0)) + normalize(e + vec3(1e-5))*r);
      float lit = max(dot(nm, gSun), 0.0);
      float mott = 0.80 + 0.32*fbm(e.xy*3.6 + 7.0);
      float disc = smoothstep(rho*1.02, rho*0.90, length(d - gMoon*md));
      col += vec3(1.0,0.97,0.90)*disc*(lit*mott*4.2 + 0.055)*moonUp;
    }
    col += vec3(0.55,0.62,0.85)*pow(max(md,0.0), 900.0)*0.55*gMoonK*moonUp;
  }
  return col;
}

// ---------------------------------------------------------------- dial face
// direction of the hour line for hour angle H (radians), in plate (x,z)
vec2 hourDir(float H){
  return normalize(vec2(abs(gSP)*sin(H), -gSGN*cos(H)) + vec2(1e-5, 0.0));
}

// engrave the dial face. returns vec4(signedDist, dirx, dirz, halfWidth)
vec4 dialField(vec2 q){
  float bestD = 1e6, bestW = 0.01;
  vec2  bestDir = vec2(0.0, 1.0);
  float r = length(q);

  // concentric rings bounding the numeral band, plus an inner ring
  for (int i = 0; i < 3; i++){
    float rr = (i == 0) ? 0.720 : ((i == 1) ? 0.955 : 0.300);
    float w  = (i == 1) ? 0.0115 : 0.0075;
    float dd = abs(r - rr) - w;
    if (dd < bestD){ bestD = dd; bestW = w; bestDir = q/max(r,1e-4)*sign(r - rr); }
  }

  // hour lines + half-hour ticks
  for (int k = 0; k < 17; k++){
    float hh = float(k) + 4.0;                       // 04:00 .. 20:00
    float Hk = (hh - 12.0)*15.0*D2R;
    vec2  hd = hourDir(Hk);
    bool noon = (k == 8);
    float w  = noon ? 0.0150 : 0.0092;
    float r1 = noon ? 0.955 : 0.720;
    vec2  o  = segOff(q, hd*0.150, hd*r1);
    float dl = length(o);
    float dd = dl - w;
    if (dd < bestD){ bestD = dd; bestW = w; bestDir = o/max(dl,1e-4); }

    vec2 td = hourDir(Hk + 7.5*D2R);                 // half-hour tick
    vec2 o2 = segOff(q, td*0.620, td*0.720);
    float dl2 = length(o2);
    float dd2 = dl2 - 0.0062;
    if (dd2 < bestD){ bestD = dd2; bestW = 0.0062; bestDir = o2/max(dl2,1e-4); }
  }

  // roman numeral for the nearest hour line (numeral band only)
  if (r > 0.700 && r < 0.980){
    vec2 u = q/max(r, 1e-4);
    float Hn = atan(u.x/max(abs(gSP), 0.06), -gSGN*u.y);
    float hf = 12.0 + degrees(Hn)/15.0;
    float hi = floor(hf + 0.5);
    if (hi >= 4.0 && hi <= 20.0){
      float rn = mod(hi - 1.0, 12.0) + 1.0;          // 1..12
      int   hn = int(rn + 0.5);
      int lead = 0; int nI = 0;
      if (hn >= 10){ lead = 3; nI = hn - 10; }
      else if (hn == 9){ lead = 0; nI = 0; }
      else if (hn >= 5){ lead = 2; nI = hn - 5; }
      else { lead = 0; nI = hn; }
      int nchar = (hn == 9) ? 2 : (((lead > 0) ? 1 : 0) + nI);

      vec2 hd = hourDir((hi - 12.0)*15.0*D2R);
      // numerals sit radially but stand upright along the noon axis, so they
      // read from the standing position on the equator-facing side of the dial
      vec2 Rw = vec2(gSGN, 0.0), Uw = vec2(0.0, -gSGN);
      float sz = 0.125;
      vec2  v  = q - hd*0.840;
      vec2  g  = vec2(dot(v, Rw), dot(v, Uw))/sz;
      float adv = 0.40;
      g.x += adv*(float(nchar) - 1.0)*0.5;           // centre the string

      for (int s = 0; s < 4; s++){
        if (s >= nchar) break;
        int c;
        if (hn == 9)       c = (s == 0) ? 1 : 3;
        else if (lead > 0) c = (s == 0) ? lead : 1;
        else               c = 1;
        vec2 gl = vec2(g.x - adv*float(s), g.y);
        for (int j = 0; j < 2; j++){
          vec2 a, b; bool on = false;
          if (c == 1){ if (j == 0){ a = vec2(0.0,-0.46); b = vec2(0.0,0.46); on = true; } }
          else if (c == 2){ float sx = (j == 0) ? -0.26 : 0.26; a = vec2(sx,0.46); b = vec2(0.0,-0.46); on = true; }
          else            { float sy = (j == 0) ? -0.46 : 0.46; a = vec2(-0.26,sy); b = vec2(0.26,-sy); on = true; }
          if (on){
            vec2 o3 = segOff(gl, a, b);
            float dl3 = length(o3)*sz;
            float dd3 = dl3 - 0.0078;
            if (dd3 < bestD){
              bestD = dd3; bestW = 0.0082;
              bestDir = (o3.x*Rw + o3.y*Uw)/max(length(o3), 1e-4);
            }
          }
        }
      }
    }
  }
  return vec4(bestD, bestDir, bestW);
}

// ---------------------------------------------------------------- solar
float solarAlt(float hourVal, float sdec, float cdec){
  float Hv = (hourVal - 12.0)*15.0*D2R;
  return sdec*gSP + cdec*gCP*cos(Hv);
}

void main(){
  vec2 uv = (gl_FragCoord.xy - 0.5*RENDERSIZE.xy)/RENDERSIZE.y;

  // ---- clock ----------------------------------------------------
  float yr, N, hh;
  if (useRealTime){
    yr = DATE.x;
    float M = DATE.y, D = DATE.z;
    float leap = ((mod(yr,4.0) < 0.5) && ((mod(yr,100.0) > 0.5) || (mod(yr,400.0) < 0.5))) ? 1.0 : 0.0;
    float K = 2.0 - leap;
    N  = floor(275.0*M/9.0) - K*floor((M + 9.0)/12.0) + D - 30.0;
    hh = DATE.w/3600.0;
  } else {
    yr = 2026.0; N = dayOfYear; hh = solarHour;
  }
  float lam = TAU*(N - 81.0)/365.0;                 // sun ecliptic longitude
  float Bq  = TAU*(N - 81.0)/364.0;
  float eot = 9.87*sin(2.0*Bq) - 7.53*cos(Bq) - 1.5*sin(Bq);   // minutes
  if (useRealTime) hh += clockOffset/60.0 + eot/60.0;          // clock -> apparent solar

  float phi = latitude*D2R;
  gSP = sin(phi); gCP = cos(phi); gSGN = (latitude >= 0.0) ? 1.0 : -1.0;

  float decl = 23.44*D2R*sin(lam);
  float sdec = sin(decl), cdec = cos(decl);
  float Hs   = (hh - 12.0)*15.0*D2R;
  gSun = normalize(vec3(-cdec*sin(Hs),
                         sdec*gSP + cdec*gCP*cos(Hs),
                       -(sdec*gCP - cdec*gSP*cos(Hs))));

  // ---- moon (mean synodic model; ignores the 5 deg orbital tilt) ----
  float y2 = yr - 2000.0;
  float leaps = floor((y2 + 3.0)/4.0) - floor((y2 + 99.0)/100.0) + floor((y2 + 399.0)/400.0);
  float days = 365.0*y2 + leaps + N - 1.0 + hh/24.0;
  float age  = fract((days - 5.76)/29.530588);
  float lamM = lam + TAU*age;
  float dM   = 23.44*D2R*sin(lamM);
  float HM   = Hs + TAU*age;
  float sdm = sin(dM), cdm = cos(dM);
  gMoon = normalize(vec3(-cdm*sin(HM),
                          sdm*gSP + cdm*gCP*cos(HM),
                        -(sdm*gCP - cdm*gSP*cos(HM))));
  gMoonK = clamp((1.0 - dot(gSun, gMoon))*0.5, 0.0, 1.0);

  // ---- gnomon geometry (style always parallel to the earth axis) ----
  float tphi = abs(gSP)/max(gCP, 0.04);
  gGL = min(0.575, 0.700/max(tphi, 0.02));
  gGH = gGL*tphi;

  // ---- lighting -------------------------------------------------
  float sa = gSun.y;
  gSunH = normalize(vec2(gSun.x, gSun.z) + vec2(1e-4));
  gStarRot = -Hs - lam;

  float am = 1.0/max(sa + 0.075, 0.075);
  // reddening from airmass, renormalised so a low sun goes red rather than merely dim
  vec3 ext = exp(-vec3(0.13, 0.36, 0.80)*(am - 1.0)*0.45);
  gSunCol = vec3(1.0, 0.92, 0.82)*ext/pow(max(max(ext.r, ext.g), ext.b), 0.95);
  gSunI   = smoothstep(-0.022, 0.030, sa)*4.60;
  gMoonCol = vec3(0.66, 0.76, 1.0);
  gMoonI   = smoothstep(-0.020, 0.045, gMoon.y)*pow(gMoonK, 1.35)*0.300;

  // sky palette keys, blended across the twilight bands
  vec3 NZ = vec3(0.0045,0.0075,0.020), NH = vec3(0.013,0.018,0.038);
  vec3 TZ = vec3(0.022,0.052,0.190),   TH = vec3(0.125,0.068,0.090);
  vec3 GZ = vec3(0.100,0.260,0.720),   GH = vec3(1.150,0.560,0.270);
  vec3 DZ = vec3(0.075,0.300,1.050),   DH = vec3(0.620,0.860,1.250);
  float ta = smoothstep(-0.310, -0.035, sa);
  float tb = smoothstep(-0.030,  0.075, sa);
  float tc = smoothstep( 0.090,  0.420, sa);
  gZen = mix(mix(mix(NZ, TZ, ta), GZ, tb), DZ, tc);
  gHor = mix(mix(mix(NH, TH, ta), GH, tb), DH, tc);
  gGlowI = smoothstep(-0.200, -0.010, sa)*(1.0 - smoothstep(0.020, 0.300, sa))*1.05
         + smoothstep(-0.010, 0.150, sa)*0.26;
  gGlowCol = mix(vec3(0.62,0.15,0.12), vec3(1.0,0.52,0.20), smoothstep(-0.09, 0.06, sa));
  gFillDir = normalize(vec3(gSunH.x, 0.13, gSunH.y));
  gFill    = gGlowCol*gGlowI*0.80;
  gNight = smoothstep(0.02, -0.16, sa);
  gStarF = smoothstep(-0.035, -0.235, sa);
  gAmb   = mix(gHor, gZen, 0.72)*0.88 + gMoonCol*gMoonI*0.30;

  // stone thermal state: lagged insolation + the season mean altitude
  float lagged = max(solarAlt(hh - 2.6, sdec, cdec), 0.0);
  float season = sdec*gSP;
  gWarm  = clamp(2.0*lagged + 1.6*season + 0.20, 0.0, 1.0);
  gFrost = clamp((0.30 - gWarm)*2.6, 0.0, 1.0)*gNight;

  float key = gSunI*max(sa, 0.0)*0.90 + gMoonI*1.60 + dot(gAmb, vec3(0.33))
            + dot(gFill, vec3(0.33))*1.30;
  gAdapt = clamp(0.48/pow(key + 0.030, 0.72), 0.12, 1.15);

  // ---- camera ---------------------------------------------------
  float az = viewAngle*D2R;
  vec3 ro = vec3(sin(az)*2.02, 0.74, gSGN*cos(az)*2.02);   // stand on the sunward side
  float pitch = 17.0*D2R;
  vec3 head = -vec3(sin(az), 0.0, gSGN*cos(az));
  vec3 fw = normalize(head*cos(pitch) - vec3(0.0,1.0,0.0)*sin(pitch));
  vec3 rt = normalize(cross(fw, vec3(0.0,1.0,0.0)));
  vec3 upv = cross(rt, fw);
  vec3 rd = normalize(uv.x*rt + uv.y*upv + FOCAL*fw);

  // ---- trace ----------------------------------------------------
  float tHit = 1e9; int kind = 0;   // 1 = stone, 2 = ground
  float b0, b1;
  if (boundHit(ro, rd, b0, b1)){
    float t = max(b0, 0.001);
    for (int i = 0; i < 96; i++){
      float d = map(ro + rd*t);
      if (d < 0.0008*t){ tHit = t; kind = 1; break; }
      t += d;
      if (t > b1) break;
    }
  }
  if (rd.y < -0.0005){
    float tg = (GROUND_Y - ro.y)/rd.y;
    if (tg > 0.0 && tg < tHit){ tHit = tg; kind = 2; }
  }

  vec3 col;
  if (kind == 0){
    col = skyColor(rd, true);
    // distant ridge: a silhouette for the sun and moon to rise and set behind
    float azr = atan(rd.x, -rd.z);
    float rn1 = fbm(vec2(azr*1.70 + 4.0, 2.0));
    float rn2 = fbm(vec2(azr*6.20 + 9.0, 6.0));
    float rh  = 0.006 + 0.060*rn1*rn1 + 0.013*rn2*rn2;
    float aaw = 1.4/RENDERSIZE.y;
    float rm  = smoothstep(rh + aaw, rh - aaw, rd.y);
    vec3 hz   = mix(gHor, gZen, 0.10);
    col = mix(col, mix(vec3(0.016,0.019,0.018), hz, 0.52), rm*0.95);
  } else {
    vec3 p = ro + rd*tHit;
    vec3 n, alb;
    float occ = 1.0, sheen = 0.0, grooveOcc = 1.0;

    if (kind == 1){
      n = calcNormal(p);
      float pxw = tHit/(FOCAL*RENDERSIZE.y)/max(abs(dot(n,rd)), 0.18);
      occ = ao(p, n);

      // limestone body: triplanar-ish coord, top face vs cylinder wall
      bool face = abs(n.y) > 0.5;
      vec2 tc = face ? p.xz*3.1 : vec2(atan(p.z, p.x)*1.55, p.y*3.4);
      float grain = fbm(tc*5.0);
      float mott  = fbm(tc*1.35 + 4.0);
      alb = mix(vec3(0.505,0.462,0.392), vec3(0.660,0.615,0.535), mott);
      alb *= 0.74 + 0.48*grain;
      float lich = smoothstep(0.54, 0.80, fbm(tc*2.3 + 19.0))*stoneWear;
      alb = mix(alb, vec3(0.255,0.290,0.195), lich*0.62*(0.30 + 0.70*max(n.y,0.0)));
      // surface relief: pitting and tooth
      vec2 bg = bumpGrad(tc*7.0)*(0.020 + 0.030*stoneWear);
      vec3 tanU, tanV;
      if (face){ tanU = vec3(1.0,0.0,0.0); tanV = vec3(0.0,0.0,1.0); }
      else     { tanU = normalize(vec3(-p.z, 0.0, p.x) + vec3(1e-4)); tanV = vec3(0.0,1.0,0.0); }
      n = normalize(n - (tanU*bg.x + tanV*bg.y));

      // engraved face
      bool top = (p.y > -0.004) && (length(p.xz) < PLATE_R - 0.012) && (n.y > 0.5);
      if (top){
        vec4 f = dialField(p.xz);
        float inside = smoothstep(pxw*0.7, -pxw*0.7, f.x);
        float u = clamp(f.x/f.w + 1.0, 0.0, 1.0);
        float depth = (1.0 - u*u)*inside*mix(1.0, 0.60, stoneWear);
        vec2 grad = 2.0*u*f.yz/f.w*0.0072*inside*mix(1.0, 0.6, stoneWear);
        n = normalize(n + vec3(-grad.x, 0.0, -grad.y));
        // a cut groove holds dirt, sees less sky, and shadows itself
        alb = mix(alb, alb*vec3(0.52,0.50,0.47), depth);
        occ *= 1.0 - 0.62*depth;
        grooveOcc = 1.0 - 0.52*depth*depth;
        sheen = depth*0.35;
      }
    } else {
      // ground: dry heath, receding into haze
      vec2 gp = p.xz;
      float patch = fbm(gp*0.55 + 2.0);
      float tuft  = fbm(gp*3.1);
      float fine  = fbm(gp*13.0 + 5.0);
      alb = mix(vec3(0.030,0.038,0.020), vec3(0.115,0.104,0.052), tuft*0.65 + patch*0.55);
      alb *= 0.50 + 0.85*fine;
      vec2 gg = bumpGrad(gp*9.0)*0.10;
      n = normalize(vec3((tuft - 0.5)*0.45 - gg.x, 1.0, (patch - 0.5)*0.45 - gg.y));
      occ = 0.68 + 0.32*tuft;
    }

    // frost on upward faces, cold nights
    float up = max(n.y, 0.0);
    float sparkleH = hash21(floor(p.xz*310.0) + floor(p.y*310.0));
    float frostM = gFrost*up*smoothstep(0.35, 0.75, fbm(p.xz*7.0 + 31.0))
                 * ((kind == 1) ? 1.0 : 0.40);
    // sparkle only survives at close range, or it aliases into static
    float spark = step(0.9986, sparkleH)*exp(-tHit*0.55);
    alb = mix(alb, vec3(0.62,0.68,0.78), frostM*0.60);

    // ---- direct light ----
    float ndlS = max(dot(n, gSun), 0.0);
    vec3 lit = vec3(0.0);
    if (ndlS > 0.001 && gSunI > 0.001){
      float sh = softShadow(p + n*0.006, gSun, 78.0);
      lit += gSunCol*gSunI*ndlS*sh;
      vec3 hv = normalize(gSun - rd);
      lit += gSunCol*gSunI*sh*pow(max(dot(n,hv),0.0), 34.0)*(0.05 + sheen)*ndlS;
    }
    float ndlM = max(dot(n, gMoon), 0.0);
    if (ndlM > 0.001 && gMoonI > 0.0006){
      float shm = softShadow(p + n*0.006, gMoon, 44.0);
      lit += gMoonCol*gMoonI*ndlM*shm;
      vec3 hv = normalize(gMoon - rd);
      lit += gMoonCol*gMoonI*shm*pow(max(dot(n,hv),0.0), 26.0)*(0.8 + 6.0*frostM)*ndlM;
    }

    // the glowing twilight horizon is a broad warm source of its own
    lit += gFill*max(dot(n, gFillDir), 0.0);
    lit *= grooveOcc;

    // ---- ambient: hemispheric sky + ground bounce ----
    vec3 amb = gAmb*(0.35 + 0.65*up)*occ;
    amb += vec3(0.10,0.11,0.07)*gAmb*0.9*(0.30 + 0.70*(1.0 - up))*occ*0.55;

    col = alb*(lit + amb);

    // stone still holding the day heat, only readable once it is dark
    col += vec3(0.075,0.031,0.014)*gWarm*gNight*occ*(0.35 + 0.65*up)*alb;
    col += vec3(0.85,0.92,1.0)*frostM*spark*(0.9 + 40.0*gMoonI);

    // aerial perspective
    float fog = 1.0 - exp(-tHit*((kind == 2) ? 0.012 : 0.008));
    vec3 fogc = skyColor(normalize(rd + vec3(0.0, 0.02, 0.0)), false);
    col = mix(col, fogc, fog*((kind == 2) ? 1.0 : 0.55));
  }

  // ---- output ---------------------------------------------------
  col *= gAdapt*exposure;
  // Purkinje shift: night vision loses colour and leans blue
  float lum = dot(col, vec3(0.2126, 0.7152, 0.0722));
  float scot = gNight*smoothstep(0.16, 0.008, lum);
  col = mix(col, vec3(lum)*vec3(0.72, 0.86, 1.22), scot*0.62);

  col = vec3(1.0) - exp(-col*1.55);
  col = pow(max(col, 0.0), vec3(0.4545));
  float l2 = dot(col, vec3(0.2126, 0.7152, 0.0722));
  col = clamp(mix(vec3(l2), col, 1.22), 0.0, 1.0);      // saturation
  col = mix(col, col*col*(3.0 - 2.0*col), 0.45);        // S-curve
  col *= 1.0 - 0.26*dot(uv, uv);                       // vignette
  col += (hash21(gl_FragCoord.xy + fract(TIME)*17.0) - 0.5)/220.0;   // dither
  gl_FragColor = vec4(clamp(col, 0.0, 1.0), 1.0);
}