/*{
  "DESCRIPTION": "Flying low over an endless high desert: Joshua trees, Mojave yuccas, creosote, cholla, brittlebush and granite boulders streaming past a couple of metres under the camera, under a sky that is solved rather than painted. There is no texture in here and no buffer. The world is an infinite field on a nine metre lattice, and every cell is one exact integer hash that decides what grows there, how big it is, how far it leans off centre and which way it faces -- so the ten thousandth tree costs exactly what the first one did. The hash is integer arithmetic carried in floats with every intermediate a whole number below 2^24, because a float hash quietly loses the bits that decide structure as its argument grows, and two GPUs will then disagree about which cells hold trees. Plants are billboards, but a billboard here is not a picture: it is a signed distance field in the plane facing the camera, so it has no resolution, its coverage is exact, it is antialiased against the real pixel footprint at its own distance, and it can be lit from behind. That last part is most of the shader -- at eleven degrees of sun elevation the light coming THROUGH a yucca rosette is brighter than the light bouncing off it, and a cholla's spines glow before its flesh does. Every radius is widened to at least one pixel first, because a twig thinner than its pixel does not go faint, it goes dashed, and two thousand plants doing that reads as the whole frame boiling. The ground is a sum of plane waves rather than fbm, which is a tenth the cost inside a march that runs it sixty times a pixel and lets every z frequency be forced to a whole multiple of 2*pi/720; on it are two soils, the pale alluvium that still moves and the dark varnished pavement that does not, braided dry washes cut through both, gravel, dry grass and windrows in the sand that fade out exactly when the pixel footprint stops resolving them. Shadows are walked along the sun direction rather than sampled from a neighbourhood, because a low sun throws a shadow five times an object's height and a 3x3 patch of cells would miss nearly all of them. The sky is single scatter by hand: Rayleigh, whose 1/lambda^4 cross section is why the zenith is blue and why the sun reddens as it drops, plus an aerosol haze with a hard forward lobe, which is the glare around the sun and the pale band along the whole horizon; on top of that a multiple-scattering floor carrying the anti-twilight arch, an altocumulus deck projected onto a plane so its rows converge to the horizon by themselves, cirrus stretched along the wind above it, and three ranges at four, thirteen and twenty-six kilometres. The clouds get their silver lining honestly, by sampling their own density one step toward the sun and subtracting it. Everything -- ground, plants, clouds, camera sway, bob, wind and grain -- is an exact 60 second loop at the default speed: 720 metres of travel, an eighty cell lattice, and every rate a whole number of cycles per minute. Sliders for speed, camera height, pitch, yaw and field of view, sun elevation, azimuth and power, haze, fog, cloud cover, cirrus, plant density, how much of the field is Joshua trees, plant scale, wind, terrain relief, sand ripples, gravel, washes, heat shimmer, mountains, shadows, backlight, saturation, exposure, glare, vignette and grain.",
  "CREDIT": "CC0",
  "ISFVSN": "2",
  "CATEGORIES": ["generator"],
  "INPUTS": [
    { "NAME": "uSpeed",     "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.0,    "MAX": 3.0,   "LABEL": "speed" },
    { "NAME": "uCamH",      "TYPE": "float", "DEFAULT": 1.70,  "MIN": 0.6,    "MAX": 12.0,  "LABEL": "camera height m" },
    { "NAME": "uPitch",     "TYPE": "float", "DEFAULT": 2.00,  "MIN": -15.0,  "MAX": 25.0,  "LABEL": "camera pitch" },
    { "NAME": "uYaw",       "TYPE": "float", "DEFAULT": 0.00,  "MIN": -45.0,  "MAX": 45.0,  "LABEL": "camera yaw" },
    { "NAME": "uFov",       "TYPE": "float", "DEFAULT": 62.0,  "MIN": 35.0,   "MAX": 95.0,  "LABEL": "field of view" },
    { "NAME": "uSunElev",   "TYPE": "float", "DEFAULT": 12.0,  "MIN": 1.0,    "MAX": 20.0,  "LABEL": "sun elevation" },
    { "NAME": "uSunAzim",   "TYPE": "float", "DEFAULT": -24.0, "MIN": -180.0, "MAX": 180.0, "LABEL": "sun azimuth" },
    { "NAME": "uSunPower",  "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.0,    "MAX": 3.0,   "LABEL": "sun power" },
    { "NAME": "uHaze",      "TYPE": "float", "DEFAULT": 0.58,  "MIN": 0.0,    "MAX": 2.0,   "LABEL": "haze" },
    { "NAME": "uFog",       "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.0,    "MAX": 3.0,   "LABEL": "aerial perspective" },
    { "NAME": "uClouds",    "TYPE": "float", "DEFAULT": 0.68,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "cloud cover" },
    { "NAME": "uCirrus",    "TYPE": "float", "DEFAULT": 0.38,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "cirrus" },
    { "NAME": "uDensity",   "TYPE": "float", "DEFAULT": 0.95,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "plant density" },
    { "NAME": "uJoshua",    "TYPE": "float", "DEFAULT": 0.50,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "joshua trees" },
    { "NAME": "uScale",     "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.5,    "MAX": 1.8,   "LABEL": "plant scale" },
    { "NAME": "uWind",      "TYPE": "float", "DEFAULT": 0.45,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "wind" },
    { "NAME": "uRelief",    "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.0,    "MAX": 2.5,   "LABEL": "terrain relief" },
    { "NAME": "uRipple",    "TYPE": "float", "DEFAULT": 0.50,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "sand ripples" },
    { "NAME": "uPebbles",   "TYPE": "float", "DEFAULT": 0.60,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "gravel" },
    { "NAME": "uWash",      "TYPE": "float", "DEFAULT": 0.80,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "dry washes" },
    { "NAME": "uShimmer",   "TYPE": "float", "DEFAULT": 0.40,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "heat shimmer" },
    { "NAME": "uMountains", "TYPE": "float", "DEFAULT": 0.80,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "mountains" },
    { "NAME": "uShadow",    "TYPE": "float", "DEFAULT": 0.85,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "shadows" },
    { "NAME": "uBacklight", "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.0,    "MAX": 2.0,   "LABEL": "backlight" },
    { "NAME": "uSat",       "TYPE": "float", "DEFAULT": 1.08,  "MIN": 0.0,    "MAX": 2.0,   "LABEL": "saturation" },
    { "NAME": "uExposure",  "TYPE": "float", "DEFAULT": 1.00,  "MIN": 0.0,    "MAX": 3.0,   "LABEL": "exposure" },
    { "NAME": "uBloom",     "TYPE": "float", "DEFAULT": 0.40,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "sun glare" },
    { "NAME": "uVignette",  "TYPE": "float", "DEFAULT": 0.00,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "vignette" },
    { "NAME": "uGrain",     "TYPE": "float", "DEFAULT": 0.25,  "MIN": 0.0,    "MAX": 1.0,   "LABEL": "grain" }
  ]
}*/
#define PI  3.14159265
#define TAU 6.28318531

// ---- globals the shading needs ---------------------------------------------
vec3  gSun;        // unit vector toward the sun
vec3  gSunCol;     // sunlight after its own trip through the atmosphere
vec3  gAmb;        // sky averaged over the hemisphere: the ambient term
vec3  gHorizon;    // colour the fog fades everything to, sampled per ray
float gPxK;        // world metres one pixel covers, per metre of distance

// ---- the two hashes --------------------------------------------------------
// Dave Hoskins' float hashes. Used only where the answer is a continuous
// quantity -- a grain of noise, a metre of jitter -- because fract() of a
// product quietly loses bits as its argument grows and two GPUs can then
// disagree about the last of them.
float h11(float p){
p = fract(p * 0.1031);
p *= p + 33.33;
p *= p + p;
return fract(p);
}
float h21(vec2 p){
vec3 p3 = fract(p.xyx * 0.1031);
p3 += dot(p3, p3.yzx + 33.33);
return fract((p3.x + p3.y) * p3.z);
}

// Anything that decides STRUCTURE -- whether a cell holds a Joshua tree or bare
// gravel -- uses this instead: integer arithmetic carried in floats where every
// intermediate is a whole number below 2^24, so float32 is exact and no two
// implementations can disagree. 4091 is prime and 3 does not divide 4090, so
// n -> n^3 is a permutation of the residues; the cube is split in two steps so
// the largest intermediate is 4090^2, just inside exact range.
#define IM 4091.0
float ihash(float n){
n = mod(n, IM);
n = mod(n * 3192.0 + 3096.0, IM); n = mod(mod(n * n, IM) * n, IM);
n = mod(n *  745.0 +  531.0, IM); n = mod(mod(n * n, IM) * n, IM);
n = mod(n * 1109.0 + 3358.0, IM); n = mod(mod(n * n, IM) * n, IM);
return n;                                    // a whole number in 0..4090
}

// ---- noise -----------------------------------------------------------------
// Value noise whose lattice wraps, so the pattern is exactly periodic. per is
// given in lattice cells; pass it scaled with the octave and the world-space
// period stays put as the frequency doubles.
float vnz(vec2 p, vec2 per){
vec2 i = floor(p), f = fract(p);
f = f * f * (3.0 - 2.0 * f);
vec2 a = mod(i, per), b = mod(i + 1.0, per);
return mix(mix(h21(vec2(a.x, a.y)), h21(vec2(b.x, a.y)), f.x),
   mix(h21(vec2(a.x, b.y)), h21(vec2(b.x, b.y)), f.x), f.y);
}
float stretch(float n, float k){ return clamp((n - 0.5) * k + 0.5, 0.0, 1.0); }
float vn1(float x, float per){
float i = floor(x), f = fract(x);
f = f * f * (3.0 - 2.0 * f);
return mix(h11(mod(i, per)), h11(mod(i + 1.0, per)), f);
}
// Ridged 1D fbm on a circle: the lattice closes on itself every per cells, so
// a horizon built from it has no seam behind the camera.
// How wide the plant can possibly be at height hh (in units of H), as a
// fraction of W. It must never be TIGHTER than the real silhouette -- that is
// the clipping bug again -- but every bit of slack in it is a candidate slot
// spent on empty air.
float profOf(float kd, float hh){
float t = clamp(hh, 0.0, 1.8);
float top = 1.0 - smoothstep(1.12, 1.46, t);         // nothing above any plant
if (kd < 2.5)                                         // Joshua: pole, then crown
 return (0.24 + 1.31 * smoothstep(0.14, 0.50, t)) * top;
if (kd < 6.5 || (kd > 9.5 && kd < 13.5))              // anything on a stem
 return (0.42 + 1.13 * smoothstep(0.04, 0.55, t)) * top;
return 1.55 * top;                                    // mounds, rocks, bone
}

float lodN(float cells){
float k = cells * gPxK * (1.0 / TAU);
return 1.0 / (1.0 + 4.0 * k * k);
}

// Ridged, but only partly: pure 1-|2v-1| on every octave gives a row of
// identical teeth, which is the one profile no mountain range has. Mixing the
// ridge back toward plain value noise, and letting the fine octaves fall away
// faster than they double, leaves a skyline with a few real summits on it.
float ridge1(float x, float per, float sd){
float a = 0.0, w = 0.5, s = 1.0, t = 0.0, rm = 0.52;
for (int i = 0; i < 7; i++){
 float lw = w * lodN(per * s);
 float v = vn1(x * s + sd, per * s);
 a += lw * mix(v, 1.0 - abs(v * 2.0 - 1.0), rm);
 t += lw;
 w *= 0.54; s *= 2.0; rm = min(rm * 1.28, 0.92);  // soft massif, sharp flanks
}
// A range is not a picket fence. Every octave ridged at full strength gives
// a row of identical spikes; an envelope at half the base frequency gives it
// a few dominant massifs with lower saddles between them, which is what
// stops every summit being the same height. per is even so the envelope's
// lattice still closes on itself round the circle.
return (a / max(t, 1e-5)) * (0.46 + 0.78 * vn1(x * 0.5 + sd + 7.0, per * 0.5));
}
float fbm2(vec2 p, vec2 per, float oct){
float a = 0.0, w = 0.5, s = 1.0;
for (int i = 0; i < 5; i++){
 if (float(i) >= oct) break;
 a += w * vnz(p * s, per * s);
 w *= 0.5; s *= 2.0;
}
return a;
}

// ---- the ground ------------------------------------------------------------
// A sum of plane waves rather than fbm, for two reasons: it is a tenth of the
// cost inside a march that evaluates it sixty times a pixel, and every z
// frequency can be forced to a whole multiple of 2*pi/720 so the terrain is
// exactly periodic over one minute of travel. The x frequencies are free.
// An alluvial fan is broad and gently graded anyway -- the interesting relief
// out here is in the plants, not the dirt.
float terr(vec2 p){
const float k = TAU / 720.0;
float h;
h  = 4.10 * sin(p.x * 0.00611 + p.y * (k *  2.0) + 0.7);
h += 2.75 * sin(p.x * 0.01471 - p.y * (k *  3.0) + 2.4);
h += 1.60 * sin(p.x * 0.00302 + p.y * (k *  5.0) - 1.1);
h += 0.90 * sin(p.x * 0.03310 + p.y * (k *  8.0) + 3.0);
h += 0.46 * sin(p.x * 0.02010 - p.y * (k * 13.0) + 0.3);
h += 0.22 * sin(p.x * 0.05170 + p.y * (k * 21.0) + 1.6);
return h * uRelief;
}

// ---- sky -------------------------------------------------------------------
// Single-scatter by hand rather than a lookup table. Two species of scatterer:
// Rayleigh, whose 1/lambda^4 cross section is what makes the zenith blue and
// what reddens the sun as it sets, and an aerosol haze with a strong forward
// lobe, which is the glare around the sun and the pale band that sits along
// the whole horizon on a dusty afternoon. Both are integrated with one number
// -- how much air the ray looks through -- which is all a flat-earth sky needs.
vec3 skyAt(vec3 d){
float mu = clamp(dot(d, gSun), -1.0, 1.0);
float y  = max(d.y, 0.0);
float vm = 1.0 / (y + 0.175);                 // view airmass, 1 at the zenith
vec3  be = vec3(0.44, 0.94, 1.98);            // Rayleigh, normalised to green

// Two beams, not one. The light that lights the SKY scatters high up, where
// the sun has crossed far less air than it has by the time it reaches the
// ground -- which is exactly why a desert zenith stays blue while the sun
// itself has already gone orange. A single airmass cannot produce both, and
// trying made the whole frame sepia.
vec3 hi   = exp(-be * (0.115 / max(gSun.y + 0.30, 0.10))) * uSunPower;
vec3 beam = mix(gSunCol, hi, smoothstep(0.0, 0.34, y));

vec3 sct = 1.0 - exp(-be * vm * 0.40);
float phR = 0.0597 * (1.0 + mu * mu);
vec3 col = sct * beam * phR * 19.5;

float g = 0.815;
float phM = (1.0 - g * g) / (12.566 * pow(max(1.0 + g * g - 2.0 * g * mu, 1e-3), 1.5));
col += (1.0 - exp(-vm * 0.026)) * gSunCol * phM * 7.5 * uHaze;

// Multiple scattering, which the single-scatter terms above cannot produce:
// it is why the sky opposite a low sun is not black and why the horizon
// holds a warm floor. The pink in it is the anti-twilight arch -- light that
// has already been reddened on its way past the terminator.
float lowSun = 1.0 - smoothstep(0.0, 0.42, gSun.y);
col += beam * vec3(0.32, 0.46, 0.72) * (0.030 + 0.016 * vm);
col += vec3(0.40, 0.17, 0.22) * lowSun * (1.0 - smoothstep(0.0, 0.26, y))
        * (0.5 - 0.5 * mu) * 0.50 * uHaze;

// Dust lifted off the floor. Only the last couple of degrees.
col = mix(col, mix(vec3(0.56, 0.46, 0.37), gSunCol, 0.50) * 0.60,
   (1.0 - smoothstep(0.0, 0.045, y)) * 0.30 * uHaze);
return max(col, vec3(0.0));
}

// The sun itself, plus the glare it throws into the lens. Kept out of skyAt so
// the fog can reuse the sky without ever fogging a scene toward a sun disc.
vec3 sunDisc(vec3 d){
float mu = clamp(dot(d, gSun), -1.0, 1.0);
float s  = length(d - gSun * mu);             // sin of the angle off centre
float R  = 0.0125;                            // wider than the real 0.0047, so it reads
float lim = sqrt(max(1.0 - (s / R) * (s / R), 0.0));
vec3 disc = gSunCol * 26.0 * (0.62 + 0.38 * lim) * smoothstep(R, R * 0.94, s);
vec3 glow = gSunCol * 1.45 * exp(-s * 46.0) + gSunCol * 0.42 * exp(-s * 8.5);
return (disc + glow) * step(0.0, d.y + 0.02);
}

// ---- clouds ----------------------------------------------------------------
// Projected onto a plane at a fixed altitude, which is what makes the rows
// converge to the horizon without anything being told to. Two decks: an
// altocumulus sheet that breaks into rows, and cirrus stretched across the
// wind above it. The lattice wraps every 24 cells in x and the wind carries it
// exactly 24 cells a minute, so the whole sky is on the same 60 second loop as
// the ground.
vec3 clouds(vec3 d, vec3 sky){
float cy = max(d.y, 0.010);
vec2  cp = d.xz / cy;
float t60 = TIME * (1.0 / 60.0);
vec2  toSun = normalize(gSun.xz + vec2(1e-4, 0.0));
vec2  per = vec2(24.0, 512.0);

// altocumulus
vec2 ap = cp * 0.42 + vec2(t60 * 24.0, 0.0);
float f = fbm2(ap, per, 5.0);
float rows = 0.5 + 0.5 * sin(cp.y * 0.55 + f * 4.0);      // breaks the sheet up
f = f * (0.62 + 0.38 * rows);
float cov = 1.0 - uClouds;
float dens = smoothstep(cov * 0.60 + 0.15, cov * 0.60 + 0.27, f);
// The lit edge: density a step toward the sun, subtracted. Where the cloud
// thins in that direction the sun is getting through, which is the whole of
// what a silver lining is.
float f2 = fbm2(ap + toSun * 0.85, per, 3.0);
float lit = clamp((f - f2) * 5.0 + 0.42, 0.0, 1.0);
// A cloud base is lit by the whole sky under it, so it is never darker
// than the haze it sits in. Painting it a flat grey made lavender bruises.
vec3 cCol = mix(gAmb * 0.62 + gSunCol * 0.16,
    gSunCol * 1.30 + gAmb * 0.35, lit);

// cirrus, stretched along the wind and much thinner
vec2 xp = cp * vec2(0.10, 0.62) + vec2(t60 * 6.0, 0.0);
float cf = fbm2(xp, vec2(6.0, 512.0), 4.0);
cf = 1.0 - abs(cf * 2.0 - 1.0);
float cd = smoothstep(0.70, 0.97, cf) * uCirrus * 0.62;
vec3 xCol = gSunCol * 1.15 + gAmb * 0.55;

float hz = smoothstep(0.014, 0.135, d.y);     // both decks die into the haze
sky = mix(sky, xCol, clamp(cd, 0.0, 1.0) * hz);
sky = mix(sky, cCol, clamp(dens, 0.0, 1.0) * hz * 0.94);
return sky;
}

// ---- the ranges on the horizon ---------------------------------------------
// Three of them, at 4, 13 and 26 km. The two far ones are drawn by direction
// alone -- no parallax, because 720 metres of travel is nothing at 26 km and
// pinning them to the world would be the one thing in the shader that did not
// close the loop. The near range does parallax, and its lattice wraps at the
// same 720 metres as everything else.
// The BODY of a range, not just its outline. A flat fill inside a jagged edge
// is what "cartoon mountain" means: a real one is mostly internal structure --
// spurs standing out of the shadow of the gullies between them, a sunward
// flank a stop brighter than the other, tops that catch the light and feet
// that drown in haze. dep is how far down the face the point is, 0 at the
// ridge line and 1 at the foot, so it costs nothing to know.
vec3 rangeBody(float az, float dep, float per, float sd2, float side,
   vec3 rock, float hz){
float l1 = lodN(per * 2.0), l2 = lodN(per * 6.0), l3 = lodN(per * 15.0);
float g1 = mix(0.5, vn1(az * (per *  2.0 / TAU) + sd2,       per *  2.0), l1);
float g2 = mix(0.5, vn1(az * (per *  6.0 / TAU) + sd2 * 3.1, per *  6.0), l2);
float g3 = mix(0.5, vn1(az * (per * 15.0 / TAU) + sd2 * 7.7, per * 15.0), l3);
float sh = 0.34 + 0.66 * g1 + 0.36 * g2 + 0.22 * g3;
sh *= 0.70 + 0.60 * side;                      // which flank faces the sun
sh *= 0.80 + 0.42 * (1.0 - dep);               // tops catch it, feet do not
return mix(rock * sh, gHorizon, clamp(hz * (0.12 + 0.88 * dep * dep), 0.0, 1.0));
}

vec3 ranges(vec3 d, vec3 sky, vec3 ro){
if (d.y <= 0.0) return sky;
float az   = atan(d.x, d.z);
float horz = max(length(d.xz), 1e-4);
float tanE = d.y / horz;
float sunAz = atan(gSun.x, gSun.z);
float side = 0.5 + 0.5 * cos(az - sunAz);

// far: almost pure haze, just a shape in it
float e0 = (4400.0 * uMountains * ridge1(az * (12.0 / TAU) + 5.0, 12.0, 0.0) - ro.y) / 26000.0;
vec3  c0 = rangeBody(az, clamp((e0 - tanE) / max(e0, 1e-5), 0.0, 1.0),
     12.0, 0.0, side, vec3(0.300, 0.330, 0.430), 0.92);
sky = mix(sky, c0, smoothstep(0.00030, -0.00030, tanE - e0));

// mid
float e1 = (2450.0 * uMountains * ridge1(az * (20.0 / TAU) + 41.0, 20.0, 0.0) - ro.y) / 13000.0;
vec3  c1 = rangeBody(az, clamp((e1 - tanE) / max(e1, 1e-5), 0.0, 1.0),
     20.0, 11.0, side, vec3(0.272, 0.248, 0.278), 0.80);
sky = mix(sky, c1, smoothstep(0.00030, -0.00030, tanE - e1));

// near, and brown enough to read as rock rather than as distance. Sampled
// by direction only: the parallax version re-drew a different mountain every
// few seconds, which is worse than no parallax by a wide margin.
float e2 = (760.0 * uMountains
   * ridge1(az * (28.0 / TAU) + 17.0, 28.0, 0.0) - ro.y) / 4200.0;
vec3  c2 = rangeBody(az, clamp((e2 - tanE) / max(e2, 1e-5), 0.0, 1.0),
     28.0, 23.0, side, vec3(0.298, 0.232, 0.188), 0.62);
sky = mix(sky, c2, smoothstep(0.00030, -0.00030, tanE - e2));
return sky;
}

// ---- plants ----------------------------------------------------------------
// Every plant is drawn as a signed distance field in the vertical plane that
// faces the camera. That is a billboard, but not a picture of one: it has no
// resolution, it is antialiased against the actual pixel footprint at its own
// distance, and its coverage is exact, so a Joshua tree a quarter of a mile off
// still resolves into separate arms instead of turning into a smear.
//
// Every radius in here is widened to at least a pixel before it is used. A
// twig thinner than the pixel it lands in does not get faint, it gets DASHED --
// it appears and vanishes depending on where the pixel centre falls -- and a
// field of two thousand plants doing that reads as the whole frame boiling.
float gPw;                                       // half a pixel, in world metres

float sdSeg(vec2 p, vec2 a, vec2 b, float ra, float rb){
vec2 pa = p - a, ba = b - a;
float h = clamp(dot(pa, ba) / max(dot(ba, ba), 1e-7), 0.0, 1.0);
return length(pa - ba * h) - max(mix(ra, rb, h), gPw);
}
// A rosette: the polar radius modulated into nl points. Yuccas and Joshua
// trees both end in one of these, a dense ball of stiff dagger leaves, and
// drawing it as a star rather than as a dozen separate leaves costs a tenth as
// much and survives being three pixels tall.
// core says how much of the radius is solid before the spikes start. A Joshua
// tree head is a dense pom-pom and a yucca is a spray of separate daggers; the
// difference between those two numbers is most of the difference between a
// Joshua tree and a palm.
float sdRos(vec2 p, vec2 c, float R, float ph, float nl, float core){
vec2 q = p - c;
float a  = atan(q.y, q.x) * (1.0 / TAU) + ph;
float ai = a * nl;
float sp = abs(fract(ai) * 2.0 - 1.0);
// Every leaf a slightly different length. Without this the head is a cog:
// perfectly equal teeth are the one thing no plant has ever grown.
float j = 0.74 + 0.46 * h11(floor(ai + 0.5) + ph * 37.0);
return length(q) - max(R * (core + (1.0 - core) * sp * j), gPw);
}
float sdBlob(vec2 p, vec2 c, vec2 r){
vec2 q = (p - c) / r;
return (length(q) - 1.0) * min(r.x, r.y);
}

// Size of each habit, in metres, before the per-plant scale.
void plantSize(float kd, float sc, out float H, out float W){
float h = 4.2, w = 1.85;
if      (kd < 0.5) { h = 4.2;  w = 1.85; }   // Joshua tree, young and squat
else if (kd < 1.5) { h = 4.9;  w = 2.20; }   // Joshua tree, the common habit
else if (kd < 2.5) { h = 6.3;  w = 2.55; }   // Joshua tree, an old one
else if (kd < 3.5) { h = 3.0;  w = 1.50; }   // Mojave yucca
else if (kd < 4.5) { h = 1.15; w = 1.15; }   // banana yucca, stemless
else if (kd < 5.5) { h = 1.60; w = 1.35; }   // creosote
else if (kd < 6.5) { h = 2.50; w = 2.00; }   // creosote, an old ring
else if (kd < 7.5) { h = 0.95; w = 1.15; }   // brittlebush
else if (kd < 8.5) { h = 1.70; w = 1.05; }   // cholla
else if (kd < 9.5) { h = 0.95; w = 1.45; }   // granite boulder
else if (kd < 10.5){ h = 1.90; w = 1.30; }   // dead wood
else if (kd < 11.5){ h = 1.60; w = 1.55; }   // prickly pear
else if (kd < 12.5){ h = 0.80; w = 0.78; }   // barrel cactus
else if (kd < 13.5){ h = 0.90; w = 0.85; }   // bunchgrass
else if (kd < 14.5){ h = 0.50; w = 0.62; }   // cow skull
else               { h = 2.90; w = 2.60; }   // large boulder
H = h * sc * uScale;
W = min(w * sc * uScale, 2.70);              // capped so a plant stays in its cell
}

// Coverage of one plant at local point q (metres, x across and y up from the
// base). alb comes back as the albedo, leaf as how much light the material
// passes -- foliage does, bark and granite do not -- and fuzz as the halo a
// cholla wears, which is the only part of any of this that glows.
float plantMask(float kd, vec2 q, float H, float W, float sd,
    out vec3 alb, out float leaf, out float fuzz, out float shade){
float s1 = fract(sd * 7.13), s2 = fract(sd * 31.7);
float s3 = fract(sd * 97.3), s4 = fract(sd * 211.7);
float dw = 1e9, dl = 1e9;                     // wood or rock, and foliage
vec3  cw = vec3(0.20, 0.16, 0.13);
vec3  cl = vec3(0.112, 0.152, 0.092);
fuzz = 0.0;
leaf = 1.0;

if (kd < 2.5){
 // Joshua tree. The thing that makes one recognisable from half a mile
 // away is not the heads, it is the ELBOW: an arm leaves the trunk almost
 // horizontally, runs out, and then turns hard upward, so the silhouette is
 // a set of crooked forearms holding fists up. Straight diagonal limbs --
 // which is what this had -- read as a palm every time. Each arm is
 // therefore two segments, and the further round the fan it sits the more
 // sideways its first one goes.
 float tw = W * 0.085;
 vec2  bs = vec2((s1 - 0.5) * 0.05 * W, -0.06 * H);
 vec2  f0 = vec2((s1 - 0.5) * 0.16 * W, H * (0.28 + 0.08 * s2));
 dw = sdSeg(q, bs, f0, tw * 1.32, tw * 0.94);

 // The dead leaves never fall off a Joshua tree: they stay as a shaggy grey
 // thatch all the way down the trunk and as a collar under every head, and
 // they are a good half of its colour. The trunk one is the capsule again,
 // widened by a scalloped function of height so the edge is ragged.
 float th = sdSeg(q, bs, mix(bs, f0, 0.93), tw * 1.46, tw * 1.20)
   - tw * (0.20 * abs(sin(q.y * 34.0 + s3 * 7.0))
    + 0.09 * abs(sin(q.y * 13.1 - s1 * 5.0)))
   / (1.0 + gPw * gPw * 700.0);

 // A Joshua tree does not carry one head per arm. Every time it flowers the
 // growing tip dies and the branch forks, so an arm that has flowered twice
 // ends in a cluster of heads, and an old tree is a candelabra of them. One
 // head per arm is a young tree at best; the fan of five was really five
 // saplings sharing a trunk. So: three primary arms, each doing the elbow
 // out and up, and each forking at its node into two secondary branches
 // that carry their own head and collar. Six heads in three pairs, which is
 // what reads as a Joshua tree rather than as a candelabrum cactus.
 for (int i = 0; i < 3; i++){
 float fi = float(i);
 float u  = h11(sd * 53.0 + fi * 2.31);
 float v2 = fract(u * 13.7 + s2);
 float an = (fi - 1.0) * (0.72 + 0.16 * s4) + (u - 0.5) * 0.52;
 vec2  d1 = vec2(sin(an), cos(an));
 vec2  el = f0 + vec2(d1.x * W * (0.34 + 0.15 * v2),
      d1.y * H * (0.07 + 0.05 * u));      // out, barely up
 // The three arms fork at clearly different heights. Real ones do, and
 // it is what keeps one arm's pair of heads out of the next arm's.
 vec2  nd = el + vec2(d1.x * W * (0.08 + 0.07 * u),
      H * (0.10 + 0.30 * v2));            // then up, to the fork
 dw = min(dw, sdSeg(q, f0, el, tw * 1.10, tw * 0.86));
 dw = min(dw, sdSeg(q, el, nd, tw * 0.86, tw * 0.64));
 for (int j = 0; j < 2; j++){
  float fj = float(j);
  float w2 = h11(sd * 89.0 + fi * 5.71 + fj * 1.93);
  // A head belongs to ONE branch and should sit clear of its sibling.
  // Two short stubs 0.2 W apart carrying heads 0.22 W across just merge
  // into a single mass, which is what the crown had become. The rule is
  // simple: the tips have to be further apart than the heads are wide,
  // so the branches get longer, the fork gets wider, and the heads get
  // smaller until that is true.
  float a2 = an + (fj - 0.5) * (0.92 + 0.34 * w2);
  vec2  tp = nd + vec2(sin(a2) * W * (0.34 + 0.14 * w2),
       H * (0.13 + 0.11 * w2) * cos(a2 * 0.45));
  dw = min(dw, sdSeg(q, nd, tp, tw * 0.64, tw * 0.44));
  dl = min(dl, sdRos(q, tp, W * (0.158 + 0.042 * w2), w2, 22.0, 0.70));
  // the collar of dead leaves sits UNDER the head and is smaller than
  // it, or it becomes the thing you see instead of the head
  th = min(th, sdRos(q, tp - vec2(0.0, W * 0.120), W * (0.124 + 0.034 * w2),
      w2 + 0.3, 16.0, 0.52));
 }
 }
 if (th < dl){ dl = th; cl = mix(cl, vec3(0.168, 0.140, 0.104), 0.90); }
 cw = vec3(0.232, 0.194, 0.152);
} else if (kd < 3.5){
 // Mojave yucca: a short trunk in a thatch of its own dead leaves, one big
 // spray of daggers on top, and -- on about half of them -- the flower
 // spike, which is a pale panicle on a bare stalk standing clear above the
 // whole plant and is the most visible thing about a yucca in season.
 vec2 tp = vec2((s1 - 0.5) * 0.10 * W, H * 0.44);
 dw = sdSeg(q, vec2(0.0, -0.05 * H), tp, W * 0.34, W * 0.26);
 dl = sdRos(q, tp + vec2(0.0, H * 0.10), W * 0.72, s2, 15.0, 0.34);
 float th = sdSeg(q, vec2(0.0, -0.05 * H), mix(vec2(0.0, 0.0), tp, 0.86),
     W * 0.30, W * 0.24)
   - W * 0.020 * abs(sin(q.y * 38.0 + s1 * 5.0))
      / (1.0 + gPw * gPw * 700.0);
 th = min(th, sdRos(q, tp + vec2(0.0, H * 0.02), W * 0.60, s3 + 0.4, 17.0, 0.40));
 if (th < dl){ dl = th; cl = vec3(0.245, 0.186, 0.120); }
 if (s4 > 0.48){
 vec2 sp = tp + vec2((s2 - 0.5) * 0.12 * W, H * 0.30);
 float fl = sdSeg(q, tp + vec2(0.0, H * 0.06), sp, W * 0.035, W * 0.022);
 fl = min(fl, sdBlob(q, sp + vec2(0.0, H * 0.14), vec2(W * 0.19, H * 0.17)));
 if (fl < dl){ dl = fl; cl = vec3(0.62, 0.60, 0.44); }
 }
 cw = vec3(0.32, 0.25, 0.17);
} else if (kd < 4.5){
 // banana yucca, no stem at all
 dl = sdRos(q, vec2(0.0, H * 0.34), W * 0.86, s1, 17.0, 0.30);
 cl = vec3(0.17, 0.22, 0.115);
} else if (kd < 6.5){
 // Creosote, the plant that actually covers the Mojave: open, wiry, and
 // more gap than plant. The leaflets go on as many SMALL tufts scattered
 // through the crown. Six big ones at a quarter of the bush width each --
 // which is what this was -- do not read as a shrub at all, they read as a
 // handful of flat green stars stuck on a stick, and at nine leaves apiece
 // they were the most conspicuous thing in the frame.
 for (int i = 0; i < 5; i++){
 float fi = float(i);
 float u1 = h11(sd * 131.0 + fi * 3.77);
 float u2 = fract(u1 * 7.13 + s2);
 float an = (fi - 2.0) * 0.44 + (u1 - 0.5) * 0.34;
 vec2  bp = vec2((u1 - 0.5) * 0.10 * W, 0.0);
 vec2  tp = bp + vec2(sin(an) * W * (0.52 + 0.38 * u2),
      H * (0.60 + 0.40 * u1) * cos(an * 0.5));
 dw = min(dw, sdSeg(q, bp, tp, W * 0.042, W * 0.011));
 // Tufts sit AT a fraction along the twig they belong to. Scattered over
 // the crown by their own hash, as they were, most of them ended up
 // hanging in mid air beside the plant with nothing holding them.
 dl = min(dl, sdRos(q, mix(bp, tp, 0.52), W * (0.095 + 0.040 * u1), u1, 7.0, 0.34));
 dl = min(dl, sdRos(q, mix(bp, tp, 0.86), W * (0.105 + 0.050 * u2), u2 + 0.4, 7.0, 0.34));
 }
 cw = vec3(0.24, 0.19, 0.14);
 cl = vec3(0.150, 0.196, 0.088);
} else if (kd < 7.5){
 // Brittlebush: a mound of small silver-grey leaves that throws yellow
 // daisies up on wire stems. Drawn as ONE ellipse at 0.42 albedo it came
 // out as a smooth pale dome the same value as the sand it stands on, with
 // no silhouette and nothing inside it. It is a cluster: overlapping leaf
 // mounds, so the outline is broken and the inside has somewhere for the
 // depth shading to bite.
 for (int i = 0; i < 5; i++){
 float fi = float(i);
 float u  = h11(sd * 67.0 + fi * 2.91);
 float w2 = fract(u * 11.7 + s2);
 dl = min(dl, sdBlob(q, vec2((u - 0.5) * 1.15 * W, H * (0.16 + 0.40 * w2)),
      vec2(W * (0.34 + 0.15 * w2), H * (0.24 + 0.13 * u))));
 }
 cl = vec3(0.243, 0.258, 0.183);
 float fx = 1e9;
 for (int i = 0; i < 4; i++){
 float fi = float(i);
 float a = (h11(sd * 91.0 + fi * 5.31) - 0.5) * 1.7;
 vec2  fp = vec2(a * W, H * (0.78 + 0.42 * h11(sd * 43.0 + fi * 2.13)));
 fx = min(fx, sdSeg(q, vec2(a * W * 0.45, H * 0.42), fp, W * 0.012, W * 0.012));
 fx = min(fx, length(q - fp) - max(W * 0.085, gPw));
 }
 if (fx < dl){ dl = fx; cl = vec3(0.86, 0.66, 0.11); }
} else if (kd < 8.5){
 // Cholla is a CHAIN OF JOINTS, not a smooth Y. Every segment is a stubby
 // cylinder that kinks off the one below it, so the outline is knobbly all
 // the way up -- drawn as three fat capsules in a line it was a cartoon
 // saguaro. A short trunk of two joints, then three limbs of two more,
 // each kinked by its own hash.
 float jr = W * 0.190;
 vec2  p0 = vec2(0.0, -0.02 * H);
 vec2  p1 = vec2((s1 - 0.5) * 0.14 * W, H * 0.21);
 vec2  p2 = p1 + vec2((s2 - 0.5) * 0.18 * W, H * 0.17);
 dl =         sdSeg(q, p0, p1, jr * 1.06, jr * 0.99);
 dl = min(dl, sdSeg(q, p1, p2, jr * 0.99, jr * 0.93));
 for (int i = 0; i < 3; i++){
 float fi = float(i);
 float u  = h11(sd * 41.0 + fi * 3.31);
 float w2 = fract(u * 17.3 + s3);
 float an = (fi - 1.0) * (0.72 + 0.20 * s2) + (u - 0.5) * 0.36;
 vec2  d1 = vec2(sin(an), cos(an));
 vec2  b0 = mix(p1, p2, 0.30 + 0.55 * u);
 vec2  b1 = b0 + vec2(d1.x * W * (0.24 + 0.12 * w2), d1.y * H * (0.15 + 0.08 * u));
 vec2  b2 = b1 + vec2(d1.x * W * (0.05 + 0.07 * w2) + (w2 - 0.5) * 0.10 * W,
      H * (0.14 + 0.09 * w2));
 dl = min(dl, sdSeg(q, b0, b1, jr * 0.93, jr * 0.86));
 dl = min(dl, sdSeg(q, b1, b2, jr * 0.86, jr * 0.70));
 }
 // Spines. Not a fat halo round the outline -- that reads as a bad glow.
 // A cholla is furred all over, so the SILHOUETTE itself is prickly: a fine
 // cross-ripple bitten out of the distance field, damped by the footprint
 // so it never becomes a sub-pixel flicker. The rim on top of it is tight,
 // and only lights up with something behind it.
 float sp = sin(q.x * (34.0 / max(W, 0.30))) * sin(q.y * (31.0 / max(H, 0.30)));
 dl -= W * 0.040 * sp * sp / (1.0 + gPw * gPw * 700.0);
 float rw = max(W * 0.050, gPw * 1.4);
 fuzz = exp(-max(dl, 0.0) / rw) * smoothstep(-rw * 0.5, rw * 0.25, dl);
 cl = vec3(0.205, 0.248, 0.128);
} else if (kd < 9.5){
 // A granite boulder, lumpy in polar radius and cut off at the dirt.
 vec2  c = vec2(0.0, H * 0.40);
 vec2  r = vec2(W, H * 0.60);
 vec2  p = (q - c) / r;
 float a = atan(p.y, p.x);
 float rr = 1.0 + 0.085 * sin(a * 5.0 + s1 * 20.0) + 0.045 * sin(a * 9.0 - s2 * 13.0);
 dw = max((length(p) - rr) * min(r.x, r.y), -q.y - H * 0.05);
 cw = mix(vec3(0.47, 0.425, 0.375), vec3(0.34, 0.295, 0.255), s3);
 leaf = 0.0;
} else if (kd < 10.5){
 // dead wood, bleached the colour of bone
 dw =         sdSeg(q, vec2(0.0, -0.05 * H), vec2(W * (s1 - 0.5) * 0.6, H * 0.55), W * 0.075, W * 0.045);
 dw = min(dw, sdSeg(q, vec2(0.0, H * 0.30), vec2(-W * 0.62, H * 0.86), W * 0.045, W * 0.018));
 dw = min(dw, sdSeg(q, vec2(0.0, H * 0.44), vec2( W * 0.70, H * 0.72), W * 0.042, W * 0.016));
 cw = vec3(0.385, 0.335, 0.268);
 leaf = 0.0;
} else if (kd < 11.5){
 // Prickly pear: overlapping pads climbing out of each other, with tunas
 // sitting on the top rim of the upper ones. The fruit is the reason this
 // plant is worth having -- a magenta dot is the only strong colour for a
 // mile and the eye goes straight to it.
 // Pads have to be a CHAIN, each stepping a good deal more than its own
 // radius from the last, or the union of them is one amorphous green blob
 // with no pads visible in it at all -- which is what this was. A step of
 // about 1.35 radii leaves a waist between each pair, and a waist is the
 // only thing that says "two pads" rather than "one lump".
 float fx = 1e9;
 vec2  pc = vec2(0.0, H * 0.19);
 float pa = (s1 - 0.5) * 0.55;
 for (int i = 0; i < 4; i++){
 float fi = float(i);
 float u  = h11(sd * 37.0 + fi * 3.11);
 vec2  r  = vec2(W * (0.225 + 0.060 * u), H * (0.160 + 0.040 * u));
 dl = min(dl, sdBlob(q, pc, r));
 if (fi > 0.5)
  fx = min(fx, sdBlob(q, pc + vec2((u - 0.5) * r.x * 1.2, r.y * 0.90),
       vec2(W * 0.062, H * 0.048)));
 pa += (u - 0.5) * 1.15;
 pc = vec2(clamp(pc.x + sin(pa) * r.x * 1.35, -W * 0.52, W * 0.52),
    pc.y + r.y * 1.35);
 }
 cl = vec3(0.175, 0.225, 0.125);
 if (fx < dl){ dl = fx; cl = mix(vec3(0.42, 0.055, 0.105),
         vec3(0.72, 0.60, 0.10), step(0.62, s3)); }
} else if (kd < 12.5){
 // Barrel cactus: ribbed, squat, and wearing a ring of flowers on the crown
 vec2  c  = vec2(0.0, H * 0.46);
 vec2  pp = (q - c) / vec2(W * 0.80, H * 0.50);
 float an = atan(pp.y, pp.x);
 dl = (length(pp) - (1.0 + 0.055 * sin(an * 13.0 + s1 * 9.0)))
 * min(W * 0.80, H * 0.50);
 dl = max(dl, -q.y - H * 0.04);
 cl = vec3(0.150, 0.205, 0.110);
 float fx = 1e9;
 for (int i = 0; i < 5; i++){
 float a2 = (float(i) - 2.0) * (0.40 + 0.10 * s2);
 fx = min(fx, length(q - (c + vec2(sin(a2) * W * 0.54,
          cos(a2) * H * 0.50 + H * 0.06)))
    - max(W * 0.10, gPw));
 }
 if (fx < dl){ dl = fx; cl = mix(vec3(0.80, 0.52, 0.06),
         vec3(0.66, 0.14, 0.16), step(0.55, s4)); }
} else if (kd < 13.5){
 // Bunchgrass. Nothing carries the ground like grass does, and a fan of
 // blades is four instructions each.
 for (int i = 0; i < 6; i++){
 float fi = float(i);
 float u  = h11(sd * 29.0 + fi * 4.73);
 float a2 = (fi - 2.5) * 0.30 + (u - 0.5) * 0.26;
 dl = min(dl, sdSeg(q, vec2((u - 0.5) * 0.16 * W, 0.0),
      vec2(sin(a2) * W * (0.50 + 0.50 * u),
       H * (0.55 + 0.45 * u) * cos(a2 * 0.55)),
      W * 0.038, W * 0.007));
 }
 cl = mix(vec3(0.335, 0.290, 0.135), vec3(0.215, 0.235, 0.105), s2);
} else if (kd < 14.5){
 // A steer skull, bleached, lying where it fell. Four things carry it and
 // the first version had two of them wrong. The horns come off the TOP
 // corners of the frontal bone, run out almost level, and then sweep UP at
 // the tips -- out-and-down is a ram, not a steer. And they are dark:
 // horn is grey-brown keratin against white bone, and that two-tone is
 // most of what makes the shape read at forty metres. The rest is a broad
 // flat frontal bone tapering into a long muzzle, with the orbits bitten
 // back out of the sides by a max() against a negated ellipse.
 leaf = 0.0;
 dw =         sdBlob(q, vec2(0.0, H * 0.64), vec2(W * 0.36, H * 0.25));
 dw = min(dw, sdBlob(q, vec2(0.0, H * 0.50), vec2(W * 0.29, H * 0.19)));
 dw = min(dw, sdSeg(q, vec2(0.0, H * 0.46),
     vec2((s1 - 0.5) * 0.10 * W, H * 0.05), W * 0.20, W * 0.115));
 dw = max(dw, -sdBlob(q, vec2(0.0, H * 0.155), vec2(W * 0.070, H * 0.105)));
 dw = max(dw, -sdBlob(q, vec2(-W * 0.285, H * 0.585), vec2(W * 0.115, H * 0.105)));
 dw = max(dw, -sdBlob(q, vec2( W * 0.285, H * 0.585), vec2(W * 0.115, H * 0.105)));
 cw = mix(vec3(0.680, 0.645, 0.575), vec3(0.560, 0.522, 0.455), s3);

 vec2 bl = vec2(-W * 0.30, H * 0.78), br = vec2(W * 0.30, H * 0.78);
 vec2 ml = vec2(-W * 0.72, H * 0.74), mr = vec2(W * 0.72, H * 0.74);
 dl =         sdSeg(q, bl, ml, W * 0.105, W * 0.062);
 dl = min(dl, sdSeg(q, ml, vec2(-W * 0.90, H * 1.00), W * 0.062, W * 0.020));
 dl = min(dl, sdSeg(q, br, mr, W * 0.105, W * 0.062));
 dl = min(dl, sdSeg(q, mr, vec2( W * 0.90, H * 1.00), W * 0.062, W * 0.020));
 cl = mix(vec3(0.250, 0.218, 0.185), vec3(0.150, 0.135, 0.125), s2);
} else {
 // A large boulder. Three harmonics on the polar radius instead of two so
 // the outline has real corners, cut off flat at the dirt, with a second
 // rock fallen against it -- one dome, however lumpy, still reads as a dome.
 vec2  c  = vec2((s1 - 0.5) * 0.10 * W, H * 0.40);
 vec2  r  = vec2(W, H * 0.62);
 vec2  pb = (q - c) / r;
 float ab = atan(pb.y, pb.x);
 float rb = 1.0 + 0.085 * sin(ab * 4.0 + s1 * 19.0)
    + 0.050 * sin(ab * 7.0 - s2 * 11.0)
    + 0.030 * sin(ab * 11.0 + s4 * 7.0);
 dw = max((length(pb) - rb) * min(r.x, r.y), -q.y - H * 0.04);
 dw = min(dw, max(sdBlob(q, c + vec2(W * (s3 - 0.5) * 1.25, -H * 0.24),
       vec2(W * 0.44, H * 0.30)), -q.y - H * 0.04));
 cw = mix(vec3(0.395, 0.358, 0.312), vec3(0.278, 0.244, 0.210), s3);
 leaf = 0.0;
}

// Whichever of the two is nearer decides both the colour and how much light
// the material passes. Bark, bone and granite are opaque; leaves are not.
float d, sh0, sh1;
if (dw < dl){ d = dw; alb = cw; leaf = min(leaf, 0.12); sh0 = 1.06; sh1 = 0.74; }
else        { d = dl; alb = cl;                        sh0 = 1.16; sh1 = 0.44; }
// How deep inside the silhouette the point is. On foliage this is worth
// more than any normal: a rosette is bright along every leaf edge and dark
// down in the middle of the crown, and -d already knows exactly where the
// leaves are, for free. Without it the heads are flat green stickers.
d += (vnz(q * (30.0 / max(W, 0.30)) + sd * 57.0, vec2(4096.0)) - 0.5)
 * W * 0.030 / (1.0 + gPw * gPw * 900.0);
float ins = clamp(-d / max(W * 0.20, gPw * 2.0), 0.0, 1.0);
shade = mix(sh0, sh1, ins * ins);
// Per-plant colour drift. A field of identically green anything is the
// fastest way to look procedural.
alb *= 0.72 + 0.55 * s3;
alb.r *= 0.94 + 0.16 * s4;
alb.g *= 0.95 + 0.12 * s2;
return clamp(smoothstep(gPw, -gPw, d) + fuzz * 0.30, 0.0, 1.0);
}

// ---- the lattice -----------------------------------------------------------
// One cell, one plant, ten metres. Seventy-two cells is 720 metres, which is
// exactly a minute of travel, so wrapping the z index at seventy-two is what
// closes the loop -- change CS and this number changes with it. Everything about the plant falls out of a single exact integer:
// what it is, how big, how far it leans off centre, which way it is turned.
#define CS 10.0

// Every column of cells is slid along z by its own amount. Without this the
// lattice is square-on to the direction of travel and the field arrives in
// ranks -- jitter inside a cell cannot hide that, because a ninth of a cell is
// all the jitter a plant has before it would cross into a neighbour and the
// ray walk would stop looking for it. Sliding whole COLUMNS costs nothing:
// the cells stay square and stay disjoint, they just no longer line up.
float coff(float cx){ return h11(cx * 1.37) * CS; }

bool cellPlant(vec2 cel, float ox, bool needY, out vec3 P, out float kd, out float H, out float W, out float sd){
float key = ihash(mod(cel.x, 512.0) * 313.0 + mod(cel.y, 72.0));
sd = key * (1.0 / IM);
kd = 11.0; H = 0.0; W = 0.0; P = vec3(0.0);
if (fract(sd * 911.0) > uDensity) return false;      // bare gravel
// Sixteen habits on an explicit weight table rather than an even split of
// 24 slots: grass wants to be everywhere and a cow skull wants to be a
// once-in-a-while, and an even split cannot say that. Both selectors are
// multiply-and-fract on an exactly representable seed, so both builds agree
// on every one of them.
if (fract(sd * 71.0) < uJoshua * 0.75){
 kd = floor(mod(key, 3.0));                         // the three Joshua habits
} else {
 float r = fract(sd * 313.0);
 kd = r < 0.085 ?  3.0 : r < 0.150 ?  4.0 : r < 0.265 ?  5.0
 : r < 0.320 ?  6.0 : r < 0.400 ?  7.0 : r < 0.455 ?  8.0
 : r < 0.550 ?  9.0 : r < 0.600 ? 10.0 : r < 0.680 ? 11.0
 : r < 0.735 ? 12.0 : r < 0.915 ? 13.0 : r < 0.923 ? 14.0
 : r < 0.968 ? 15.0 : 16.0;
}
if (kd > 15.5) return false;                         // bare ground
float sc = 0.66 + 0.80 * fract(sd * 173.0);
plantSize(kd, sc, H, W);
// A plant may lean off the middle of its cell, but never far enough to cross
// the boundary: the ray walk visits cells, so anything hanging into the next
// one would simply not be looked for and would blink out at certain angles.
float jit = clamp(CS * 0.5 - W * 1.55 - 0.45, 0.0, 2.10);
vec2  b = vec2((cel.x + 0.5) * CS, ox + (cel.y + 0.5) * CS)
  + (vec2(fract(sd), fract(sd * 53.0)) - 0.5) * (2.0 * jit);
P = vec3(b.x, needY ? terr(b) : 0.0, b.y);
return true;
}

// How much of a ground point is in the shade of something growing up-sun of
// it. Walked along the sun direction rather than over a square of cells,
// because at eleven degrees a Joshua tree throws a shadow five times its own
// height and a 3x3 neighbourhood would miss almost all of them. max() rather
// than a product, so a cell sampled twice does not darken twice.
float plantShadow(vec3 G){
vec2  s2 = normalize(gSun.xz + vec2(1e-5, 0.0));
float slope = max(gSun.y, 0.03) / max(length(gSun.xz), 1e-4);
float sh = 0.0;
for (int i = 0; i < 9; i++){
 vec2 sp = G.xz + s2 * (float(i) * CS * 0.52 + CS * 0.26);
 float cx = floor(sp.x / CS), ox = coff(cx);
 vec3 P; float kd, H, W, sd;
 if (!cellPlant(vec2(cx, floor((sp.y - ox) / CS)), ox, false, P, kd, H, W, sd)) continue;
 vec2  a = P.xz, ba = -s2 * (H / slope);
 vec2  pa = G.xz - a;
 float h = clamp(dot(pa, ba) / max(dot(ba, ba), 1e-6), 0.0, 1.0);
 float d = length(pa - ba * h);
 // Narrow at the near end and wide at the far one, because the near end is
 // the shadow of the trunk and the far end is the shadow of the crown. It
 // was the other way round, which is the one shape a tree shadow is not.
 float rad = W * (0.20 + 0.85 * smoothstep(0.22, 0.92, h));
 sh = max(sh, smoothstep(rad, rad * 0.40, d) * (1.0 - smoothstep(0.72, 1.02, h)));
}
return 1.0 - sh * uShadow * 0.86;
}

// ---- fog -------------------------------------------------------------------
// Clear desert air, so the linear term is small; the quadratic one is the dust
// that sits in the bottom hundred metres and is what actually dissolves the
// far end of the field. Everything fades to the sky colour at its own azimuth,
// which is why the ground can run out at 440 metres without a visible edge.
float fogAmt(float t){
return 1.0 - exp(-(t * 0.0016 + t * t * 1.05e-6) * uFog);
}

// ---- ground ----------------------------------------------------------------
vec3 terrN(vec2 p){
const float e = 1.4;
float h = terr(p);
return normalize(vec3(h - terr(p + vec2(e, 0.0)), e, h - terr(p + vec2(0.0, e))));
}

vec3 groundShade(vec3 G, vec3 rd, float t){
vec2  p  = G.xz;
float px = t * gPxK;                          // metres one pixel covers here
float dg = 1.0 / (1.0 + px * px * 250.0);     // 4cm grain: gone by ~60m

// Two soils: pale sandy alluvium washed down off the fan, and the darker
// varnished pavement it leaves behind where nothing has moved in centuries.
float m1 = stretch(vnz(p * (1.0 / 90.0), vec2(4096.0,   8.0)), 2.6);
float m2 = stretch(vnz(p * (1.0 / 22.5), vec2(4096.0,  32.0)), 2.4);
float gr0 = vnz(p * 22.0,          vec2(4096.0, 15840.0));
float gr = stretch(gr0, 2.6);
float pt = stretch(vnz(p * (1.0 / 6.0),  vec2(4096.0, 120.0)), 2.2);
float mn = stretch(vnz(p * (1.0 / 3.0),  vec2(4096.0, 240.0)), 2.2);
float fn = stretch(vnz(p * 0.9,          vec2(4096.0, 648.0)), 2.2);
float dg2 = 1.0 / (1.0 + px * px * 0.4);      // metre-scale detail lasts longer
// Pink-tan alluvium against grey desert pavement: the two soils differ in
// hue, not just in value, which is what stops a mile of dirt reading as one
// flat wash.
// Three tones, not one wash: pale alluvium that still moves, grey varnished
// pavement that has not moved in centuries, and the red of oxidised iron
// where it has.
vec3  alb = mix(vec3(0.480, 0.348, 0.226), vec3(0.205, 0.178, 0.153),
    smoothstep(0.22, 0.86, m1 * 0.62 + m2 * 0.38));
alb = mix(alb, vec3(0.395, 0.222, 0.134), smoothstep(0.42, 0.88, 1.0 - m1) * 0.50);

// Dry washes: braided, stretched down the fan, and paler than what they cut
// through because the sand in them is the part that still moves.
float wv = vnz(p * vec2(1.0 / 26.0, 1.0 / 144.0) + vec2(m2 * 0.7, 0.0), vec2(4096.0, 5.0));
float wash = 1.0 - smoothstep(0.02, 0.13, abs(wv - 0.5));
alb = mix(alb, vec3(0.575, 0.470, 0.330), wash * uWash * 0.85);
dg *= 1.0 - 0.55 * wash * uWash;               // swept smooth, no gravel in it

alb *= 0.92 + 0.16 * gr * dg;                                  // grain
alb *= 0.84 + 0.32 * pt;                                       // patchy ground
alb *= 0.90 + 0.20 * mn * dg2;                                 // and at a stride
// Where the fines have blown out, what is left is grey and slightly blue.
alb = mix(alb, alb * vec3(0.76, 0.82, 0.91), smoothstep(0.40, 0.88, fn) * 0.42 * dg2);
alb = mix(alb, vec3(0.585, 0.508, 0.252), smoothstep(0.58, 0.90, gr) * 0.45 * dg * uPebbles);

// Windrows in the loose sand. Damped by the pixel footprint rather than by
// distance, so they fade out exactly when they stop being resolvable.
// Same rule as the terrain: the z component of the wave number has to be a
// whole multiple of 2*pi/720 or the windrows do not come back after a
// minute of travel -- 0.570*4.19 was 273.67 turns over the loop, and every
// ripple in the frame was a sixth of a wavelength out of step at t=60.
const float rkz = 274.0 * (TAU / 720.0);
const vec2 rdir = normalize(vec2(3.444, rkz));
float rph  = p.x * 3.444 + p.y * rkz + m2 * 3.0;
float rdmp = 1.0 / (1.0 + px * px * 17.5);

vec3 n = terrN(p);
float ga = vnz(p * 22.0 + vec2(0.16, 0.0), vec2(4096.0, 15840.0)) - gr0;
float gb = vnz(p * 22.0 + vec2(0.0, 0.16), vec2(4096.0, 15840.0)) - gr0;
n = normalize(n + vec3(rdir.x, 0.0, rdir.y) * (cos(rph) * 0.20 * rdmp * uRipple)
    + vec3(ga, 0.0, gb) * (0.30 * dg * uPebbles));

// Loose cobbles, one lattice of their own at 1.2 m. A stone is smaller than
// its cell and cannot leave it, so a single lookup is exact -- no 3x3 -- and
// the normal is the real dome, not a bump: e is the offset in radii, so the
// height is sqrt(1 - e.e) and the normal falls straight out of it.
// Every stone an identical hemisphere on an even lattice catches a low sun
// in exactly the same place, and a field of identical bright crescents reads
// as spilled coins. So: fewer, bigger jitter, each one elongated and turned
// by its own hash, and the dome deliberately shallow -- a stone half buried
// in sand never presents a vertical face anyway, and it was the vertical
// faces going to full sunlight that made them glow.
vec2  rc = floor(p / 1.2);
float rh = h21(vec2(rc.x, mod(rc.y, 600.0)) + 3.7);
if (rh > 0.44){
 float rr = 0.040 + 0.105 * fract(rh * 13.0);
 vec2  rp = (rc + 0.5) * 1.2 + (vec2(fract(rh * 31.0), fract(rh * 97.0)) - 0.5) * 0.80;
 float ra = fract(rh * 53.0) * PI;
 float ca_ = cos(ra), sa_ = sin(ra);
 vec2  el = vec2(rr, rr * (0.60 + 0.45 * fract(rh * 17.0)));
 vec2  ee = p - rp;
 ee = vec2(ee.x * ca_ + ee.y * sa_, ee.y * ca_ - ee.x * sa_) / el;
 float h2 = 1.0 - dot(ee, ee);
 // Contact shadow: the same ellipse tested one stone-width up-sun. Without
 // it a stone is a picture of a stone lying on nothing.
 vec2  so = p + normalize(gSun.xz + vec2(1e-4, 0.0)) * rr * 1.5 - rp;
 so = vec2(so.x * ca_ + so.y * sa_, so.y * ca_ - so.x * sa_) / el;
 alb *= 1.0 - 0.24 * smoothstep(1.25, 0.55, dot(so, so)) * step(h2, 0.0)
     / (1.0 + px * px * 14.0);
 if (h2 > 0.0){
 float rk = smoothstep(0.0, 0.18, h2) / (1.0 + px * px * 14.0);
 // A dome is rotationally symmetric, so turning one changes nothing about
 // how it catches the light -- every stone lit its near flank in exactly
 // the same place, which is what made them read as spilled coins. Tipping
 // the whole stone by its own hash is what breaks that.
 vec2 tl = (vec2(fract(rh * 131.0), fract(rh * 197.0)) - 0.5) * 1.05;
 vec3 sn = normalize(vec3(ee.x * 0.50 + tl.x, sqrt(h2) + 0.42, ee.y * 0.50 + tl.y));
 n = normalize(mix(n, sn, rk * 0.62));
 alb = mix(alb, mix(vec3(0.285, 0.248, 0.210), vec3(0.140, 0.126, 0.114),
      fract(rh * 71.0)), rk * 0.62);
 }
}

float sh  = plantShadow(G);
float ndl = max(dot(n, gSun), 0.0);
vec3  col = alb * (gSunCol * ndl * 6.4 * sh
    + gAmb * (0.34 + 0.66 * max(n.y, 0.0)) * 0.58
    + gSunCol * vec3(1.0, 0.66, 0.36) * 0.13 * sh);
col += gSunCol * pow(1.0 - max(dot(n, -rd), 0.0), 5.0) * 0.09;  // dust sheen

// Beyond the range the billboards are walked out to, the plants have to keep
// going somehow. They stop being individuals and become a tint on the dirt,
// which is what a mile of creosote looks like anyway.
col = mix(col, col * vec3(0.575, 0.625, 0.475),
   smoothstep(120.0, 340.0, t) * uDensity * 0.72);
return col;
}

// ---- drawing one plant -----------------------------------------------------
vec3 drawPlant(vec3 col, vec4 A, vec4 B, vec3 ro, vec3 rd){
vec3  P = A.xyz;
float kd = A.w, H = B.x, W = B.y, sd = B.z, tt = B.w;
gPw = 0.5 * tt * gPxK;

float dP = length(vec2(ro.x - P.x, ro.z - P.z));   // one distance per plant
vec3 nb = normalize(vec3(ro.x - P.x, 0.0, ro.z - P.z));
vec2 rb = vec2(-nb.z, nb.x);
vec3 Q  = ro + rd * tt;
vec2 q  = vec2(dot(Q.xz - P.xz, rb), Q.y - P.y);

// Wind. Bending the sample point instead of the geometry costs one line and
// the tip moves as the square of height, which is roughly what a cantilever
// does. Seven cycles a minute, so it lands back where it started.
float ph = fract(sd * 61.0);
float v  = clamp(q.y / max(H, 0.01), 0.0, 1.0);
// Granite does not move in the wind, and neither does a skull or a dead
// branch. Grass moves more than anything else out here.
float wk = (kd < 8.5) ? 1.0 : ((kd > 12.5 && kd < 13.5) ? 1.7 : 0.0);
q.x -= wk * uWind * 0.30 * (0.5 + 0.9 * ph) * v * v
 * sin(TIME * (TAU * 7.0 / 60.0) + ph * 40.0 + P.z * (TAU * 8.0 / 720.0));

vec3 alb; float leaf, fuzz, shade;
float a = plantMask(kd, q, H, W, sd, alb, leaf, fuzz, shade);
a *= smoothstep(0.35, 1.90, dP) * (1.0 - smoothstep(205.0, 378.0, dP));
if (a <= 0.002) return col;

// The billboard is flat, so the normal has to be invented: a cylinder about
// the plant's own vertical axis, which is close enough to true for anything
// that grows in a rosette and is what stops the whole field reading as
// cardboard.
// The fake normal is a cylinder about the plant's own vertical axis, which
// is close enough to true for anything that grows in a rosette. A boulder is
// not one of those, so it gets the vertical half of the bulge as well and
// comes out domed and lit from above instead of flat.
float sph = step(8.5, kd) * step(kd, 9.5);
vec3 n = normalize(nb * 0.80
    + vec3(rb.x, 0.0, rb.y) * (clamp(q.x / max(W, 0.01), -1.0, 1.0) * 1.15)
    + vec3(0.0, mix(0.30, (v - 0.42) * 2.3, sph), 0.0));
float ndl = max(dot(n, gSun), 0.0);
float ao  = mix(0.52, 1.0, v);                       // dark down in the litter
float bl  = clamp(dot(rd, gSun), 0.0, 1.0);
bl = bl * bl * bl;

vec3 lit = alb * (gSunCol * ndl * 5.6
    + gAmb * (0.34 + 0.42 * n.y) * 0.66
    + gSunCol * vec3(1.0, 0.66, 0.38) * 0.30 * (1.0 - 0.6 * v));
// Light straight through the leaf. At this sun angle it is most of what
// makes a yucca look like a yucca rather than a green blob.
float tex = vnz(q * (13.0 / max(W, 0.35)) + sd * 91.0, vec2(4096.0));
lit *= shade * ao * (0.80 + 0.40 * mix(0.5, tex, 1.0 / (1.0 + gPw * gPw * 900.0)));
// The leaf edges are where the light gets through, so the backlight rides
// the INVERSE of the same depth term.
lit += alb * gSunCol * bl * bl * leaf * 1.9 * uBacklight * (1.55 - 0.55 * shade);
lit += gSunCol * fuzz * (0.10 + 1.5 * bl) * uBacklight * 0.50;
lit += gSunCol * pow(1.0 - abs(dot(n, rd)), 3.0) * 0.22 * (0.3 + 0.7 * bl);

lit = mix(lit, gHorizon, fogAmt(tt));
return mix(col, lit, a);
}

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

// ---- the frame -------------------------------------------------------------
void main(){
vec2 fragCoord = gl_FragCoord.xy;
vec2 R  = RENDERSIZE.xy;
vec2 uv = (fragCoord - 0.5 * R) / R.y;
float T = TIME;

// ---- sun. The colour is Beer-Lambert through however much air the beam has
// to cross to reach the ground, which is the whole reason a low sun is
// orange: the blue end has already been scattered out of it somewhere else.
float se = radians(uSunElev), sa = radians(uSunAzim);
gSun    = normalize(vec3(sin(sa) * cos(se), sin(se), cos(sa) * cos(se)));
gSunCol = exp(-vec3(0.44, 0.94, 1.98) * (1.0 / max(gSun.y, 0.045)) * 0.115) * uSunPower;

// ---- camera. 720 metres a minute, hugging the ground, with a slow sway and
// a little bob -- every rate a whole number of cycles per minute.
vec3 ro = vec3(sin(T * (TAU * 2.0 / 60.0)) * 1.7, 0.0, uSpeed * 12.0 * T);
ro.y = terr(ro.xz) + uCamH * (1.0 + 0.055 * sin(T * (TAU * 9.0 / 60.0)));
float pit = radians(uPitch) + 0.0065 * sin(T * (TAU *  5.0 / 60.0));
float yaw = radians(uYaw)   + 0.0110 * sin(T * (TAU *  3.0 / 60.0));
vec3 fw = normalize(vec3(sin(yaw) * cos(pit), sin(pit), cos(yaw) * cos(pit)));
vec3 rt = normalize(cross(vec3(0.0, 1.0, 0.0), fw));
vec3 up = cross(fw, rt);
// uv runs to +-0.5, so the half angle needs a half in front of it. Get this
// wrong and the shader quietly renders at twice the stated field of view.
float fl = 0.5 / tan(radians(uFov) * 0.5);
vec3 rd = normalize(uv.x * rt + uv.y * up + fl * fw);
gPxK = 2.0 * tan(radians(uFov) * 0.5) / R.y;

// Heat coming off the ground, which only bends what is close to the horizon.
float sm = uShimmer * 0.0026 * exp(-abs(rd.y) * 24.0);
rd.y += sm * (vnz(vec2(uv.x * 85.0, T * 3.0), vec2(4096.0, 180.0)) - 0.5);
rd = normalize(rd);

// ---- background
// Ambient is the sky averaged over the dome, not the zenith. At this sun
// angle the zenith is the DIMMEST part of the sky, and using it left every
// plant in the field a black cut-out.
gAmb = skyAt(vec3(0.0, 1.0, 0.0)) * 0.46
 + skyAt(normalize(vec3( gSun.x, 0.42,  gSun.z))) * 0.22
 + skyAt(normalize(vec3(-gSun.x, 0.42, -gSun.z))) * 0.32;
gHorizon = skyAt(normalize(vec3(rd.x, 0.030, rd.z)));
vec3 sky = skyAt(rd);
sky += sunDisc(rd) * (0.55 + uBloom);
sky = clouds(rd, sky);
sky = ranges(rd, sky, ro);

// ---- ground. A heightfield march with a big geometric step and a linear
// solve on the way out. The relief is gentle enough that a coarse step
// cannot skip a hill, and the alternative -- enough steps to be safe at a
// grazing angle -- is three times the cost for nothing visible.
float tg = -1.0, t = 0.5, lt = 0.0, ld = uCamH;
for (int i = 0; i < 60; i++){
 vec3 pp = ro + rd * t;
 float d = pp.y - terr(pp.xz);
 if (d < 0.0){ tg = lt + (t - lt) * ld / max(ld - d, 1e-5); break; }
 lt = t; ld = d;
 t += 0.085 * t + 0.55;
 if (t > 640.0) break;
}

vec3 col = sky;
if (tg > 0.0){
 gPw = 0.5 * tg * gPxK;
 col = groundShade(ro + rd * tg, rd, tg);
 col = mix(col, gHorizon, fogAmt(tg));
} else if (rd.y < 0.0){
 col = mix(gHorizon, vec3(0.42, 0.36, 0.30) * (gSunCol.r + 0.3), 0.16);
}

// ---- the plant field. A DDA over the lattice in the ground plane, keeping
// the three nearest billboards the ray actually clips. Compositing all of
// them inside the walk would mean inlining the whole plant library forty
// eight times over; three is enough because a plant is confined to its own
// cell, so the only way to stack more than three is to look down a row.
vec2 rd2 = rd.xz;
vec2 sg  = sign(rd2);
vec2 iv  = 1.0 / max(abs(rd2), vec2(1e-5));
vec2 cel = floor(ro.xz / CS);
vec2 tdl = CS * iv;
float ox = coff(cel.x);
cel.y = floor((ro.z - ox) / CS);
vec2 tmx = vec2(((sg.x * 0.5 + 0.5) * CS + sg.x * (cel.x * CS - ro.x)) * iv.x,
    ((sg.y * 0.5 + 0.5) * CS + sg.y * (cel.y * CS + ox - ro.z)) * iv.y);
float far = min(tg > 0.0 ? tg : 1e9, 400.0);

vec4 A0 = vec4(0.0), A1 = vec4(0.0), A2 = vec4(0.0), A3 = vec4(0.0), A4 = vec4(0.0);
vec4 B0 = vec4(0.0, 0.0, 0.0, 1e9), B1 = B0, B2 = B0, B3 = B0, B4 = B0;

for (int i = 0; i < 64; i++){
 vec3 P; float kd, H, W, sd;
 if (cellPlant(cel, ox, true, P, kd, H, W, sd)){
 vec3  nb = vec3(ro.x - P.x, 0.0, ro.z - P.z);
 float dh = length(nb);
 if (dh > 0.35){
  nb /= dh;
  float dn = dot(rd, nb);
  if (dn < -1e-3){
  float tt = dot(P - ro, nb) / dn;
  if (tt > 0.15 && tt < far){
   vec3  Q  = ro + rd * tt;
   vec2  rb = vec2(-nb.z, nb.x);
   float qx = dot(Q.xz - P.xz, rb);
   float qy = Q.y - P.y;
   float prof = profOf(kd, qy / max(H, 0.01));
   if (qy > -0.30 && abs(qx) < W * prof){
   if (tt < B0.w){
    A4 = A3; B4 = B3; A3 = A2; B3 = B2; A2 = A1; B2 = B1; A1 = A0; B1 = B0;
    A0 = vec4(P, kd); B0 = vec4(H, W, sd, tt);
   } else if (tt < B1.w){
    A4 = A3; B4 = B3; A3 = A2; B3 = B2; A2 = A1; B2 = B1;
    A1 = vec4(P, kd); B1 = vec4(H, W, sd, tt);
   } else if (tt < B2.w){
    A4 = A3; B4 = B3; A3 = A2; B3 = B2;
    A2 = vec4(P, kd); B2 = vec4(H, W, sd, tt);
   } else if (tt < B3.w){
    A4 = A3; B4 = B3;
    A3 = vec4(P, kd); B3 = vec4(H, W, sd, tt);
   } else if (tt < B4.w){
    A4 = vec4(P, kd); B4 = vec4(H, W, sd, tt);
   }
   }
  }
  }
 }
 }
 if (tmx.x < tmx.y){
 if (tmx.x > far) break;
 // Stepping into a new column moves the z boundaries with it, so the z
 // cell and its next crossing have to be re-solved from where the ray
 // actually is, not carried over.
 float te = tmx.x;
 cel.x += sg.x; tmx.x += tdl.x;
 ox = coff(cel.x);
 float zz = ro.z + rd.z * te;
 cel.y = floor((zz - ox) / CS);
 tmx.y = te + ((sg.y * 0.5 + 0.5) * CS + sg.y * (cel.y * CS + ox - zz)) * iv.y;
 } else {
 if (tmx.y > far) break; cel.y += sg.y; tmx.y += tdl.y;
 }
}

if (B4.w < 1e8) col = drawPlant(col, A4, B4, ro, rd);
if (B3.w < 1e8) col = drawPlant(col, A3, B3, ro, rd);
if (B2.w < 1e8) col = drawPlant(col, A2, B2, ro, rd);
if (B1.w < 1e8) col = drawPlant(col, A1, B1, ro, rd);
if (B0.w < 1e8) col = drawPlant(col, A0, B0, ro, rd);

// ---- ravens. Six of them, each wheeling on its own circle at a whole
// number of turns a minute so the flock closes with everything else. They
// are billboards like the plants, drawn last but depth-tested against the
// ground so one cannot float over a hill it is behind.
for (int i = 0; i < 6; i++){
 float fi = float(i);
 float u  = h11(fi * 7.13 + 3.0), v = h11(fi * 11.7 + 9.0);
 float ph = T * (TAU * floor(1.0 + 2.0 * u) / 60.0) + v * TAU;
 float rr = 26.0 + 54.0 * v;
 vec3  B  = vec3(cos(ph) * rr + (u - 0.5) * 110.0,
     ro.y + 11.0 + 32.0 * u + sin(ph * 2.0) * 2.5,
     ro.z + 55.0 + 150.0 * v + sin(ph) * rr);
 vec3  tb = B - ro;
 float dB = length(tb);
 vec3  ub = tb / dB;
 float ca = dot(rd, ub);
 if (ca < 0.975 || (tg > 0.0 && dB > tg)) continue;
 vec3  rgt = normalize(cross(vec3(0.0, 1.0, 0.0), ub));
 vec3  upv = cross(ub, rgt);
 vec3  hit = ro + rd * (dB / max(ca, 1e-3));
 vec2  lq  = vec2(dot(hit - B, rgt), dot(hit - B, upv));
 gPw = 0.5 * dB * gPxK;
 float fl = sin(T * TAU + fi * 2.1);                  // 60 wingbeats a minute
 float db = sdSeg(lq, vec2(0.0, 0.0), vec2(-0.86, 0.34 * fl), 0.070, 0.020);
 db = min(db, sdSeg(lq, vec2(0.0, 0.0), vec2(0.86, 0.34 * fl), 0.070, 0.020));
 db = min(db, sdSeg(lq, vec2(0.0, -0.14), vec2(0.0, 0.20), 0.060, 0.038));
 col = mix(col, vec3(0.035, 0.033, 0.040),
   smoothstep(gPw, -gPw, db) * 0.94 * (1.0 - fogAmt(dB) * 0.7));
}

// ---- post
col *= uExposure;
col = aces(col);
col = mix(vec3(dot(col, vec3(0.2126, 0.7152, 0.0722))), col, uSat);
vec2 vuv = (fragCoord / R - 0.5) * vec2(1.86, 2.0);   // aspect independent
col *= 1.0 - uVignette * 0.85 * clamp(pow(length(vuv), 2.6), 0.0, 1.0);
col += (h21(fragCoord + mod(floor(T * 30.0), 1800.0)) - 0.5) * uGrain * 0.055;
gl_FragColor = vec4(clamp(col, 0.0, 1.0), 1.0);
}