hex_terrain source

← hex_terrain · API reference · version 0.1.8

3 file(s), 1613 lines. This is the source the published version ships, and it is loft — the library is written in the same language you call it from.

Public items

A worked example is a call site the author chose as the thing to read first. Call sites are every line naming the item, found mechanically over the source with comments excluded — so an item with no worked example is not an item nobody uses.

src/hex_terrain.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// hex_terrain — the OVERLAND terrain layer of the hex_* family.
//
// A coarse overland hex lattice (one cell ~ hundreds of meters, the moros
// convention via hex_grid) is the terrain AUTHORITY; the fine, walked world is
// derived from it by pure functions.  The model:
//
//  * CELLS carry height / material / moisture / water-flow direction +
//    accumulation / relief — the same record shape at every scale.
//  * SEAMLESS MERGE: continuous fields at a fine point are a smooth kernel
//    blend of the overland cells around it — two adjacent tiles merge by
//    construction, no stitching anywhere.
//  * TERRAIN TYPES are content: each has TWO independent numbers — RISE (the
//    general lift it takes from tile relief) and STEEP (its internal detail
//    jaggedness, which also decides who WINS a boundary: steeper wins).
//    Between any two types runs a crisp fractal LINE, not a wide crossfade.
//  * RELIEF: a cell's relief = its max height delta to neighbors, COMPOUNDED
//    by how many neighbors are steep; water-flow cells contribute ZERO relief,
//    so the dry land rises around the rivers and they cut deep for free.
//  * RIVERS: coarse hydrology (priority-flood fill, steepest descent, flow
//    accumulation) gives the truth; the geometry is three levels of diversity:
//    a canonical left/right EDGE BIT per shared edge (never exits mid-side),
//    a hash-jittered interior CONTROL point per tile (confluences for free),
//    then fractal midpoint displacement whose amplitude GROWS as the land
//    flattens toward the sea.  Courses are TRIMMED to attach exactly to the
//    sea/lake shorelines.  Every course carves its own valley (additive;
//    slope-shaped: gorges in steep land, wide floodplains on the flats) and a
//    near-sea-level depth bonus floods the carved mouth: the sea comes inland.
//
// THE INVARIANT, AND THE HALF IT DOES NOT COVER.  SAMPLING is pure: a sampled
// value is a function of (terrain data, params, types, x, y) and nothing else —
// no hidden state, no sequential RNG, no raster propagation — so the fine world
// derives on demand, in any order, from any viewport.  That is what "window
// independence" names, and the tests hold it by re-deriving samples from scratch.
//
// ⚠ THE RECORD ITSELF IS NOT INDEPENDENT OF THE MAP'S EXTENT, and the phrase
// invites the opposite reading.  Hydrology is a GLOBAL pass and accumulation is
// a catchment property, so widening the map raises accumulation everywhere,
// which flips cells to water, which zeroes their relief, which moves the fine
// surface through the blend kernel.  Measured on one authored ramp: a 5x5 and an
// 8x8 window differ at 24 of 27 interior sample points, by up to 15.1 m, and the
// smaller map has NO rivers where the bigger one has five.  So generate the
// overland map ONCE, at its full extent — a window is a unit of SAMPLING, never
// a unit of generation.
//
// Style: globally-unique pub names (terrain_*); pre-allocated vectors with
// index writes for hot data.

use hex_grid::*;

// ── Tuning records ─────────────────────────────────────────────────────────

// One terrain type = content (a bundle row in the consuming game).
// tt_rise and tt_steep are INDEPENDENT; tt_wet marks water bodies (sea/lake):
// they take no relief and no boundary contest.
pub struct TerrainType {
  pub tt_key:   text,
  pub tt_rise:  float,   // m of ridge lift at the 100 m reference relief
  pub tt_steep: float,   // m of detail amplitude; sqrt(steep) wins boundaries
  pub tt_wet:   boolean,
}

// Algorithm parameters (engine side).  terrain_params() gives the blueprint's
// user-tuned defaults; every field is meters unless noted.
pub struct TerrainParams {
  tp_seed:        integer,
  tp_tile:        float,    // overland hex pitch
  pub tp_sea:         float,    // sea level
  tp_lake_fill:   float,    // flood-fill raise that makes a cell a lake
  pub tp_acc_min:     float,    // accumulation (tiles) before a stream is a river
  tp_meander:     float,    // fractal displacement strength K
  tp_seg_min:     float,    // stop subdividing below this segment length
  pub tp_nb_base:     float,    // relief multiplier with zero steep neighbors
  pub tp_nb_boost:    float,    // ADDED per steep neighbor (compounding)
  pub tp_steep_delta: float,    // neighbor height delta that counts as steep
  tp_carve_max:   float,    // max single-course valley depth
  tp_valley_w:    float,    // smallest valley falloff width
  tp_valley_max:  float,    // biggest river's falloff width
  tp_detail_wl:   float,    // detail-noise base wavelength
  tp_ridge_wl:    float,    // ridge-field wavelength (peaks COMBINE)
  tp_warp:        float,    // domain-warp displacement of the detail
}

// The knobs a generated terrain is built from: a seed and a tile size, plus the thresholds
// that follow.  Hold one and the terrain is reproducible — same params, same world.
pub fn terrain_params(seed: integer, tile: float) -> TerrainParams {
  TerrainParams { tp_seed: seed, tp_tile: tile, tp_sea: 0.0, tp_lake_fill: 4.0,
    tp_acc_min: 6.0, tp_meander: 0.34, tp_seg_min: 28.0,
    tp_nb_base: 0.6, tp_nb_boost: 0.2, tp_steep_delta: 60.0,
    tp_carve_max: 22.0, tp_valley_w: 46.0, tp_valley_max: 130.0,
    tp_detail_wl: 420.0, tp_ridge_wl: 950.0, tp_warp: 240.0 }
}

// ── The overland lattice ───────────────────────────────────────────────────

// Parallel per-cell vectors, index = r * tr_nx + c (odd-r offset coords).
pub struct Terrain {
  pub tr_nx:     integer,
  pub tr_ny:     integer,
  pub tr_h:      vector<float>,   // height (m); after hydrology = the FILLED surface
  pub tr_mat:    vector<integer>, // index into the caller's TerrainType vector
  pub tr_moist:  vector<float>,
  pub tr_flow:   vector<integer>, // steepest-descent neighbor dir 0..5, -1 = none
  pub tr_acc:    vector<float>,   // flow accumulation (tiles)
  pub tr_relief: vector<float>,   // compounded neighbor relief (0 on wet cells)
  pub tr_wet:    vector<boolean>, // water-flow cell: wet material OR river acc
}

// An empty nx-by-ny terrain: height, material and moisture per cell, all parallel, so a
// cell index means the same thing in each.
pub fn terrain_new(nx: integer, ny: integer) -> Terrain {
  h: vector<float> = [];
  mat: vector<integer> = [];
  moist: vector<float> = [];
  flow: vector<integer> = [];
  acc: vector<float> = [];
  relief: vector<float> = [];
  wet: vector<boolean> = [];
  for _ in 0..nx * ny {
    h += [0.0];
    mat += [0];
    moist += [0.5];
    flow += [-1];
    acc += [1.0];
    relief += [0.0];
    wet += [false];
  }
  Terrain { tr_nx: nx, tr_ny: ny, tr_h: h, tr_mat: mat, tr_moist: moist,
            tr_flow: flow, tr_acc: acc, tr_relief: relief, tr_wet: wet }
}

// The flat index of column c, row r.  The one place the row-major layout is written down.
pub fn terrain_idx(t: Terrain, c: integer, r: integer) -> integer {
  r * t.tr_nx + c
}

// Cell center in world meters — hex_grid's odd-r lattice scaled so the
// horizontal pitch equals tp_tile (hex_grid's unit pitch is SQRT3).
//
// Example: @HXT-001
pub fn terrain_center(p: TerrainParams, c: integer, r: integer) -> (float, float) {
  (hx, hy) = hex_to_px(c, r);
  s = p.tp_tile / SQRT3;
  (hx * s, hy * s)
}

// The downstream cell of (c, r) — itself when it has no flow direction.
pub fn terrain_flow_to(t: Terrain, c: integer, r: integer) -> (integer, integer) {
  d = t.tr_flow[r * t.tr_nx + c] ?? -1;
  if d < 0 {
    return (c, r);
  }
  hex_neighbor(c, r, d)
}

// ── Hashed noise (deterministic, overflow-trap-safe) ───────────────────────
// Every intermediate stays under 2^31 via masks, so products stay < 2^63.

fn th_hash(seed: integer, ix: integer, iy: integer, ch: integer) -> integer {
  a = (ix * 73856093) & 0x7FFFFFFF;
  b = (iy * 19349663) & 0x7FFFFFFF;
  c = (ch * 83492791) & 0x7FFFFFFF;
  s = (seed * 1597334677) & 0x7FFFFFFF;
  n = a ^ b ^ c ^ s;
  n = n ^ (n >> 13);
  n = (n * 1274126177) & 0x7FFFFFFF;
  n ^ (n >> 16)
}

// A reproducible value in 0..1 for this seed and cell.  A hash, not a random: the same
// inputs give the same answer on every machine and every run, which is what makes a
// generated world shareable by seed alone.
pub fn terrain_hash01(seed: integer, ix: integer, iy: integer, ch: integer) -> float {
  (th_hash(seed, ix, iy, ch) & 0x7FFFFFFF) as float / 2147483648.0
}

// Smooth value noise at a world point, at the given wavelength.  Built on `terrain_hash01`,
// so it inherits the same reproducibility.
pub fn terrain_vnoise(seed: integer, x: float, y: float, wl: float, ch: integer) -> float {
  u = (x / wl) ?? 0.0;
  v = (y / wl) ?? 0.0;
  fu = floor(u);
  fv = floor(v);
  iu = fu as integer;
  iv = fv as integer;
  du = u - fu;
  dv = v - fv;
  su = du * du * (3.0 - 2.0 * du);
  sv = dv * dv * (3.0 - 2.0 * dv);
  a = terrain_hash01(seed, iu, iv, ch);
  b = terrain_hash01(seed, iu + 1, iv, ch);
  c = terrain_hash01(seed, iu, iv + 1, ch);
  d = terrain_hash01(seed, iu + 1, iv + 1, ch);
  (a * (1.0 - su) + b * su) * (1.0 - sv) + (c * (1.0 - su) + d * su) * sv
}

// Fractal sum, roughly in [-1, 1].
pub fn terrain_fbm(seed: integer, x: float, y: float, wl0: float, octaves: integer, ch: integer) -> float {
  out = 0.0;
  amp = 1.0;
  wl = wl0;
  tot = 0.0;
  for o in 0..octaves {
    out = out + amp * (terrain_vnoise(seed, x, y, wl, ch * 16 + o) * 2.0 - 1.0);
    tot = tot + amp;
    amp = amp * 0.5;
    wl = wl * 0.5;
  }
  (out / tot) ?? 0.0
}

// Detail noise with DOMAIN WARP (the sampling space is swirled by a second
// noise so the value-noise lattice never reads as a repeating pattern).
pub fn terrain_detail_at(p: TerrainParams, x: float, y: float) -> float {
  wx = terrain_fbm(p.tp_seed, x, y, 700.0, 2, 80);
  wy = terrain_fbm(p.tp_seed, x, y, 700.0, 2, 81);
  terrain_fbm(p.tp_seed, x + p.tp_warp * wx, y + p.tp_warp * wy, p.tp_detail_wl, 8, 3)
}

// LOW-frequency ridged field, crest-sharpened: roughly ONE combined crest per
// steep zone — peaks consolidate into ridgelines instead of speckling.
pub fn terrain_ridge_at(p: TerrainParams, x: float, y: float) -> float {
  wx = terrain_fbm(p.tp_seed, x, y, 1500.0, 2, 84);
  wy = terrain_fbm(p.tp_seed, x, y, 1500.0, 2, 85);
  r = 1.0 - abs(terrain_fbm(p.tp_seed, x + 300.0 * wx, y + 300.0 * wy, p.tp_ridge_wl, 2, 86));
  r * r * r
}

// ── Hydrology on the lattice ───────────────────────────────────────────────

// PRIORITY-FLOOD pit fill + steepest-descent flow + accumulation.  After this
// every land cell drains to the sea; cells the flood raised by more than
// tp_lake_fill become lakes (material set to lake_mat).  Selection-based
// flood (O(n^2)) — overland grids are small and this runs once at gen time.
//
// Example: @HXT-003
pub fn terrain_hydrology(t: Terrain, p: TerrainParams, lake_mat: integer) {
  nx = t.tr_nx;
  ny = t.tr_ny;
  n = nx * ny;
  th = &t.tr_h;       // live refs: index-writes go through to the record (C86: a
  tm = &t.tr_mat;     // plain bind now COPIES, so `&` is required to write through)
  tf = &t.tr_flow;
  ta = &t.tr_acc;
  orig: vector<float> = [];
  level: vector<float> = [];
  state: vector<integer> = [];   // 0 = untouched, 1 = open, 2 = closed
  for i in 0..n {
    orig += [th[i] ?? 0.0];
    level += [1.0e18];
    state += [0];
  }
  for r in 0..ny {
    for c in 0..nx {
      i = r * nx + c;
      edge = c == 0 || r == 0 || c == nx - 1 || r == ny - 1;
      seeded = edge;
      if (th[i] ?? 0.0) <= p.tp_sea {
        seeded = true;
      }
      if seeded {
        if i < len(level) {
          level[i] = th[i] ?? 0.0;
        }
        if i < len(state) {
          state[i] = 1;
        }
      }
    }
  }
  done = false;
  while !done {
    best = -1;
    bl = 1.0e18;
    for i in 0..n {
      if (state[i] ?? 2) == 1 && (level[i] ?? 1.0e18) < bl {
        bl = level[i] ?? 1.0e18;
        best = i;
      }
    }
    if best < 0 {
      done = true;
    } else {
      if best < len(state) { state[best] = 2; }
      bc = (best % nx) ?? 0;
      br = ((best - bc) / nx) ?? 0;
      for d in 0..6 {
        (nc, nr) = hex_neighbor(bc, br, d);
        if nc >= 0 && nc < nx && nr >= 0 && nr < ny {
          j = nr * nx + nc;
          if (state[j] ?? 2) != 2 {
            nl = max(orig[j] ?? 0.0, bl + 0.05);
            if nl < (level[j] ?? 1.0e18) {
              if j < len(level) {
                level[j] = nl;
              }
            }
            if j < len(state) {
              state[j] = 1;
            }
          }
        }
      }
    }
  }
  for i in 0..n {
    if i < len(th) {
      th[i] = level[i] ?? 0.0;
    }
    raised = (level[i] ?? 0.0) > (orig[i] ?? 0.0) + p.tp_lake_fill;
    if raised && (level[i] ?? 0.0) > p.tp_sea {
      if i < len(tm) {
        tm[i] = lake_mat;
      }
    }
  }
  // steepest-descent flow on the filled surface — no pits remain
  for r in 0..ny {
    for c in 0..nx {
      i = r * nx + c;
      if i < len(tf) { tf[i] = -1; }
      if i < len(ta) { ta[i] = 1.0; }
      if (th[i] ?? 0.0) > p.tp_sea {
        bh = th[i] ?? 0.0;
        bi = -1;
        for d in 0..6 {
          (nc, nr) = hex_neighbor(c, r, d);
          if nc >= 0 && nc < nx && nr >= 0 && nr < ny && (th[nr * nx + nc] ?? 1.0e18) < bh {
            bh = th[nr * nx + nc] ?? 0.0;
            bi = d;
          }
        }
        if i < len(tf) { tf[i] = bi; }
      }
    }
  }
  // accumulation: high cells drain first (selection over a worklist)
  taken: vector<boolean> = [];
  for _ in 0..n {
    taken += [false];
  }
  for _ in 0..n {
    best = -1;
    bh = 0.0 - 1.0e18;
    for i in 0..n {
      if !(taken[i] ?? true) && (th[i] ?? 0.0) > bh {
        bh = th[i] ?? 0.0;
        best = i;
      }
    }
    if best >= 0 {
      if best < len(taken) { taken[best] = true; }
      bc = (best % nx) ?? 0;
      br = ((best - bc) / nx) ?? 0;
      bi = tf[best] ?? -1;
      if bi >= 0 && (th[best] ?? 0.0) > p.tp_sea {
        (nc, nr) = hex_neighbor(bc, br, bi);
        if nc >= 0 && nc < nx && nr >= 0 && nr < ny {
          j = nr * nx + nc;
          if j < len(ta) {
            ta[j] = (ta[j] ?? 1.0) + (ta[best] ?? 1.0);
          }
        }
      }
    }
  }
}

// Wet flags + the compounded neighbor relief.  Call AFTER the caller has
// classified land materials (relief reads the final wet set).
//
// Example: @HXT-004
pub fn terrain_relief_pass(t: Terrain, types: vector<TerrainType>, p: TerrainParams) {
  nx = t.tr_nx;
  ny = t.tr_ny;
  tw = &t.tr_wet;
  tl = &t.tr_relief;
  for r in 0..ny {
    for c in 0..nx {
      i = r * nx + c;
      wetmat = false;
      mi = t.tr_mat[i] ?? 0;
      if mi < len(types) {
        ty = types[mi];
        if ty != null {
          wetmat = ty.tt_wet;
        }
      }
      if i < len(tw) {
        tw[i] = wetmat || (t.tr_h[i] ?? 0.0) <= p.tp_sea
             || (t.tr_acc[i] ?? 1.0) >= p.tp_acc_min;
      }
    }
  }
  for r in 0..ny {
    for c in 0..nx {
      i = r * nx + c;
      if i < len(tl) { tl[i] = 0.0; }
      if !(t.tr_wet[i] ?? false) {
        m = 0.0;
        nsteep = 0;
        for d in 0..6 {
          (nc, nr) = hex_neighbor(c, r, d);
          if nc >= 0 && nc < nx && nr >= 0 && nr < ny {
            dh = abs((t.tr_h[i] ?? 0.0) - (t.tr_h[nr * nx + nc] ?? 0.0));
            m = max(m, dh);
            if dh > p.tp_steep_delta {
              nsteep = nsteep + 1;
            }
          }
        }
        if i < len(tl) {
          tl[i] = m * (p.tp_nb_base + p.tp_nb_boost * (nsteep as float));
        }
      }
    }
  }
}

// ── The seamless kernel blend ──────────────────────────────────────────────

// One blended surface record at a fine point.
pub struct TerrainSurf {
  pub sf_h:      float,    // fine height BEFORE river carving
  pub sf_blend:  float,    // the coarse blended height alone
  pub sf_mat:    integer,  // winning land material (the boundary contest)
  pub sf_mix:    float,    // 0.5..1: how decisively it wins (edge softness)
  pub sf_relief: float,
  pub sf_moist:  float,
  pub sf_watery: float,    // blended weight of wet-material cells
  pub sf_lakey:  float,    // blended weight of wet cells ABOVE sea level
}

// Smooth kernel blend of the overland cells around (x, y) + the boundary-line
// contest + the fine height.  Pure function — THE seam invariant lives here.
//
// Example: @HXT-002
// Example: @HXT-008
pub fn terrain_surface_at(t: Terrain, types: vector<TerrainType>, p: TerrainParams, x: float, y: float) -> TerrainSurf {
  nt = len(types);
  wm: vector<float> = [];
  for _ in 0..nt {
    wm += [0.0];
  }
  s = p.tp_tile / SQRT3;
  vs = 1.5 * s;
  rk = p.tp_tile * 1.02;
  rk2 = rk * rk;
  ry0 = round((y / vs) ?? 0.0) as integer;
  wsum = 0.0;
  hsum = 0.0;
  relsum = 0.0;
  moisum = 0.0;
  watsum = 0.0;
  laksum = 0.0;
  for dr in 0..3 {
    r = ry0 + dr - 1;
    off = 0.5 * ((r & 1) as float);
    c0 = floor(((x / p.tp_tile) ?? 0.0) - off) as integer;
    for dc in 0..4 {
      c = c0 + dc - 1;
      cx = ((c as float) + off) * p.tp_tile;
      cy = (r as float) * vs;
      d2 = (x - cx) * (x - cx) + (y - cy) * (y - cy);
      if d2 < rk2 {
        w = 1.0 - ((d2 / rk2) ?? 0.0);
        w = w * w;
        ci = min(max(c, 0), t.tr_nx - 1);
        ri = min(max(r, 0), t.tr_ny - 1);
        i = ri * t.tr_nx + ci;
        wsum = wsum + w;
        hsum = hsum + w * (t.tr_h[i] ?? 0.0);
        relsum = relsum + w * (t.tr_relief[i] ?? 0.0);
        moisum = moisum + w * (t.tr_moist[i] ?? 0.5);
        mi = t.tr_mat[i] ?? 0;
        if mi < len(wm) {
          wm[mi] = (wm[mi] ?? 0.0) + w;
        }
        wetm = false;
        if mi < len(types) {
          tym = types[mi];
          if tym != null {
            wetm = tym.tt_wet;
          }
        }
        if wetm || (t.tr_h[i] ?? 0.0) <= p.tp_sea {
          watsum = watsum + w;
          if (t.tr_h[i] ?? 0.0) > p.tp_sea {
            laksum = laksum + w;
          }
        }
      }
    }
  }
  iw = (1.0 / max(wsum, 1.0e-12)) ?? 1.0;
  hb = hsum * iw;
  rel = relsum * iw;
  // the boundary-line contest: per-type weights, fractally jittered scores,
  // STEEPER WINS the contested band (bias by sqrt of the type's STEEP)
  m1 = 0;
  s1 = -1.0;
  s2 = -1.0;
  steep1 = 0.0;
  steep2 = 0.0;
  rise1 = 0.0;
  rise2 = 0.0;
  for cm in 0..nt {
    tyc = types[cm];
    if tyc != null && !tyc.tt_wet && (wm[cm] ?? 0.0) > 0.0 {
      jit = 0.72 + 0.56 * (terrain_fbm(p.tp_seed, x, y, 340.0, 3, 40 + cm) * 0.5 + 0.5);
      steep_root = sqrt(max(tyc.tt_steep, 0.01));
      sc = (wm[cm] ?? 0.0) * iw * (steep_root ?? 0.0) * jit;
      if sc > s1 {
        s2 = s1;
        steep2 = steep1;
        rise2 = rise1;
        m1 = cm;
        s1 = sc;
        steep1 = tyc.tt_steep;
        rise1 = tyc.tt_rise;
      } else {
        if sc > s2 {
          s2 = sc;
          steep2 = tyc.tt_steep;
          rise2 = tyc.tt_rise;
        }
      }
    }
  }
  if s2 < 0.0 {
    s2 = s1;
    steep2 = steep1;
    rise2 = rise1;
  }
  tt = 1.0;
  if s1 > 0.0 {
    tt = min(max(((s1 - s2) / (s1 * 0.18)) ?? 1.0, 0.0), 1.0);
  }
  mt = 0.5 + 0.5 * tt;
  steep_eff = steep1 * mt + steep2 * (1.0 - mt);
  rise_eff = rise1 * mt + rise2 * (1.0 - mt);
  detail = terrain_detail_at(p, x, y);
  ridge = terrain_ridge_at(p, x, y);
  hf = hb + steep_eff * detail + rise_eff * (rel / 100.0) * ridge;
  TerrainSurf { sf_h: hf, sf_blend: hb, sf_mat: m1, sf_mix: mt, sf_relief: rel,
                sf_moist: moisum * iw, sf_watery: watsum * iw, sf_lakey: laksum * iw }
}

// Lightweight coarse-blend height only (for slopes along river courses).
//
// Example: @HXT-005
// Example: @HXT-006
pub fn terrain_blend_h(t: Terrain, p: TerrainParams, x: float, y: float) -> float {
  s = p.tp_tile / SQRT3;
  vs = 1.5 * s;
  rk = p.tp_tile * 1.02;
  rk2 = rk * rk;
  ry0 = round((y / vs) ?? 0.0) as integer;
  wsum = 0.0;
  hsum = 0.0;
  for dr in 0..3 {
    r = ry0 + dr - 1;
    off = 0.5 * ((r & 1) as float);
    c0 = floor(((x / p.tp_tile) ?? 0.0) - off) as integer;
    for dc in 0..4 {
      c = c0 + dc - 1;
      cx = ((c as float) + off) * p.tp_tile;
      cy = (r as float) * vs;
      d2 = (x - cx) * (x - cx) + (y - cy) * (y - cy);
      if d2 < rk2 {
        w = 1.0 - ((d2 / rk2) ?? 0.0);
        w = w * w;
        ci = min(max(c, 0), t.tr_nx - 1);
        ri = min(max(r, 0), t.tr_ny - 1);
        wsum = wsum + w;
        hsum = hsum + w * (t.tr_h[ri * t.tr_nx + ci] ?? 0.0);
      }
    }
  }
  (hsum / max(wsum, 1.0e-12)) ?? 0.0
}

// The fractally jittered lake contour field (>= 0.5 = lake water).
pub fn terrain_lake_field(p: TerrainParams, x: float, y: float, lakey: float) -> float {
  lakey + 0.11 * terrain_fbm(p.tp_seed, x, y, 330.0, 4, 70)
}

// Is (x, y) open water (sea by the fine height contour, or lake)?  PRE-carve —
// this is the shoreline rivers attach to.
//
// Example: @HXT-007
pub fn terrain_water_at(t: Terrain, types: vector<TerrainType>, p: TerrainParams, x: float, y: float) -> boolean {
  sf = terrain_surface_at(t, types, p, x, y);
  sf.sf_h <= p.tp_sea || terrain_lake_field(p, x, y, sf.sf_lakey) >= 0.5
}

// ── Rivers ─────────────────────────────────────────────────────────────────

// One river course PIECE: fine polyline (ends ON the shorelines) with
// per-vertex water width, carve depth and valley falloff width.
pub struct TerrainRiver {
  pub rv_x:  vector<float>,
  pub rv_y:  vector<float>,
  rv_w:  vector<float>,    // water half-width test: dist < rv_w/2
  rv_d:  vector<float>,    // carve depth
  rv_vw: vector<float>,    // carve falloff width
  rv_x0: float,   // bbox (+ tp_valley_max*4 = the influence apron)
  rv_y0: float,
  rv_x1: float,
  rv_y1: float,
}

fn th_cellhash(seed: integer, a1: integer, a2: integer, b1: integer, b2: integer, ch: integer) -> float {
  n = (a1 * 73856093) & 0x7FFFFFFF;
  n = n ^ ((a2 * 19349663) & 0x7FFFFFFF);
  n = n ^ ((b1 * 83492791) & 0x7FFFFFFF);
  n = n ^ ((b2 * 1297357571) & 0x7FFFFFFF);
  n = n ^ ((ch * 668265263) & 0x7FFFFFFF);
  n = n ^ ((seed * 1597334677) & 0x7FFFFFFF);
  n = n ^ (n >> 13);
  n = (n * 1274126177) & 0x7FFFFFFF;
  ((n ^ (n >> 16)) & 0x7FFFFFFF) as float / 2147483648.0
}

// THE EDGE BIT (user rule): a river never exits mid-side; each shared edge
// carries one canonical left/right bit.  Both tiles hash the SORTED cell pair,
// so they derive the same crossing point — seamless.
pub fn terrain_edge_cross(p: TerrainParams, ca: integer, ra: integer, cb: integer, rb: integer) -> (float, float) {
  c1 = ca;
  r1 = ra;
  c2 = cb;
  r2 = rb;
  if ca > cb || (ca == cb && ra > rb) {
    c1 = cb;
    r1 = rb;
    c2 = ca;
    r2 = ra;
  }
  (x1, y1) = terrain_center(p, c1, r1);
  (x2, y2) = terrain_center(p, c2, r2);
  mx = (x1 + x2) * 0.5;
  my = (y1 + y2) * 0.5;
  d = sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
  px = 0.0 - (((y2 - y1) / d) ?? 0.0);
  py = ((x2 - x1) / d) ?? 0.0;
  side = p.tp_tile / SQRT3;          // hex side length = circumradius
  off = 0.30 * side;
  if th_cellhash(p.tp_seed, c1, r1, c2, r2, 5) < 0.5 {
    off = 0.0 - off;
  }
  (mx + px * off, my + py * off)
}

// The interior routing point of a tile (hash-jittered off-center; the natural
// confluence spot — tributaries entering the tile share it).
pub fn terrain_control_pt(p: TerrainParams, c: integer, r: integer) -> (float, float) {
  (cx, cy) = terrain_center(p, c, r);
  side = p.tp_tile / SQRT3;
  jx = (th_cellhash(p.tp_seed, c, r, c, r, 6) - 0.5) * 0.8 * side;
  jy = (th_cellhash(p.tp_seed, c, r, c, r, 7) - 0.5) * 0.8 * side;
  (cx + jx, cy + jy)
}

// Meander factor: flat land winds hard, steep land runs straight; extra wind
// near sea level (the user rule: more curves in the flat lands by the sea).
fn th_meander(t: Terrain, p: TerrainParams, x: float, y: float) -> float {
  e = 60.0;
  hx = terrain_blend_h(t, p, x + e, y) - terrain_blend_h(t, p, x - e, y);
  hy = terrain_blend_h(t, p, x, y + e) - terrain_blend_h(t, p, x, y - e);
  sl = (sqrt(hx * hx + hy * hy) / (2.0 * e)) ?? 0.0;
  f = 0.12 + 0.88 * exp(0.0 - sl / 0.016);
  h = terrain_blend_h(t, p, x, y);
  if h > 0.0 && h < 60.0 {
    f = min(1.25, f * (1.0 + 0.5 * (1.0 - h / 60.0)));
  }
  f
}

// Fractal midpoint displacement, recursive over one anchor segment.  Appends
// the start point + the displaced interior to (ox, oy, ov) — never the end
// (the caller closes the chain).  Hashes the QUANTIZED endpoints, no
// traversal index: window-independent by construction.
fn th_subdiv(t: Terrain, p: TerrainParams, ax: float, ay: float, av: float, bx: float, by: float, bv: float, lvl: integer, ox: vector<float>, oy: vector<float>, ov: vector<float>) {
  seg = sqrt((bx - ax) * (bx - ax) + (by - ay) * (by - ay));
  if lvl > 7 || seg < p.tp_seg_min * 2.0 {
    ox += [ax];
    oy += [ay];
    ov += [av];
  } else {
    px = 0.0 - (((by - ay) / seg) ?? 0.0);
    py = ((bx - ax) / seg) ?? 0.0;
    mx = (ax + bx) * 0.5;
    my = (ay + by) * 0.5;
    hq = th_cellhash(p.tp_seed, round(ax) as integer, round(ay) as integer,
                     round(bx) as integer, round(by) as integer, 100 + lvl);
    d = (hq - 0.5) * 2.0 * seg * p.tp_meander * th_meander(t, p, mx, my);
    mv = (av + bv) * 0.5;
    th_subdiv(t, p, ax, ay, av, mx + px * d, my + py * d, mv, lvl + 1, ox, oy, ov);
    th_subdiv(t, p, mx + px * d, my + py * d, mv, bx, by, bv, lvl + 1, ox, oy, ov);
  }
}

// Bisect the exact shoreline crossing between a dry and a wet point.
fn th_shore_cross(t: Terrain, types: vector<TerrainType>, p: TerrainParams, dx0: float, dy0: float, wx0: float, wy0: float) -> (float, float) {
  ddx = dx0;
  ddy = dy0;
  wwx = wx0;
  wwy = wy0;
  for _ in 0..9 {
    mx = (ddx + wwx) * 0.5;
    my = (ddy + wwy) * 0.5;
    if terrain_water_at(t, types, p, mx, my) {
      wwx = mx;
      wwy = my;
    } else {
      ddx = mx;
      ddy = my;
    }
  }
  ((ddx + wwx) * 0.5, (ddy + wwy) * 0.5)
}

// Per-vertex shaping: water width from accumulation; carve depth slope-shaped
// (gorges steep, floodplains flat) + the ESTUARY bonus near sea level (the
// carved mouth dips below the sea -> the sea floods inland, size-proportional).
fn th_shape_river(t: Terrain, p: TerrainParams, vx: vector<float>, vy: vector<float>, vv: vector<float>) -> TerrainRiver {
  ww: vector<float> = [];
  dd: vector<float> = [];
  vw: vector<float> = [];
  x0 = vx[0] ?? 0.0;
  y0 = vy[0] ?? 0.0;
  x1 = x0;
  y1 = y0;
  e = 60.0;
  for k in 0..len(vx) {
    kx = vx[k] ?? 0.0;
    ky = vy[k] ?? 0.0;
    kv = vv[k] ?? 1.0;
    // `sqrt` answers `float?` (a negative argument has no real root) and so does a
    // division, so every value below that flows from one is nullable until it is
    // discharged.  Left undischarged they reach the `vector<float>` stores at the
    // bottom of the loop, where a null would be written as the float null sentinel
    // and read back as a depth or a width -- so each is discharged at the point it
    // is created, the same way `kx` / `ky` / `kv` are one line up.
    //
    // An accumulation is a count of drained tiles, so `kv` is never negative and
    // this root always exists; 0.0 is the value a zero accumulation would give.
    skv = sqrt(kv) ?? 0.0;
    hx = terrain_blend_h(t, p, kx + e, ky) - terrain_blend_h(t, p, kx - e, ky);
    hy = terrain_blend_h(t, p, kx, ky + e) - terrain_blend_h(t, p, kx, ky - e);
    sl = (sqrt(hx * hx + hy * hy) / (2.0 * e)) ?? 0.0;
    // Both divisors below are non-zero constants, so neither quotient can be null;
    // the fallbacks say what a null WOULD mean -- a flat slope, and no drop below
    // the reference height -- rather than being a value picked to silence a check.
    fs = min(1.0, (sl / 0.05) ?? 0.0);
    h0 = terrain_blend_h(t, p, kx, ky);
    fall = (max(h0, 0.0) / 14.0) ?? 0.0;
    base = 3.0 + 1.7 * skv;
    dep = base * (0.45 + 0.55 * fs) + base * 0.90 * exp(0.0 - fall);
    dd += [min(p.tp_carve_max, dep)];
    vw += [min(p.tp_valley_max, (24.0 + 11.0 * skv) * (1.6 - 0.6 * fs))];
    ww += [min(26.0, 2.2 + 1.35 * skv)];
    x0 = min(x0, kx);
    y0 = min(y0, ky);
    x1 = max(x1, kx);
    y1 = max(y1, ky);
  }
  ap = p.tp_valley_max * 4.0;
  TerrainRiver { rv_x: vx, rv_y: vy, rv_w: ww, rv_d: dd, rv_vw: vw,
                 rv_x0: x0 - ap, rv_y0: y0 - ap, rv_x1: x1 + ap, rv_y1: y1 + ap }
}

// All river course pieces, fully shaped: anchors (edge bits + control points)
// -> fractal displacement -> corner-cut smoothing -> shoreline trim -> width/
// depth/valley-width per vertex (slope-shaped + the estuary depth bonus).
//
// Example: @HXT-002
pub fn terrain_rivers(t: Terrain, types: vector<TerrainType>, p: TerrainParams) -> vector<TerrainRiver> {
  nx = t.tr_nx;
  ny = t.tr_ny;
  n = nx * ny;
  out: vector<TerrainRiver> = [];
  isr: vector<boolean> = [];
  feeds: vector<boolean> = [];
  for i in 0..n {
    isr += [(t.tr_acc[i] ?? 1.0) >= p.tp_acc_min && (t.tr_h[i] ?? 0.0) > p.tp_sea];
    feeds += [false];
  }
  for r in 0..ny {
    for c in 0..nx {
      i = r * nx + c;
      if (isr[i] ?? false) && (t.tr_flow[i] ?? -1) >= 0 {
        (nc, nr) = hex_neighbor(c, r, t.tr_flow[i] ?? 0);
        if nc >= 0 && nc < nx && nr >= 0 && nr < ny && (isr[nr * nx + nc] ?? false) {
          j = nr * nx + nc;
          if j < len(feeds) {
            feeds[j] = true;
          }
        }
      }
    }
  }
  for hr in 0..ny {
    for hc in 0..nx {
      hi = hr * nx + hc;
      if (isr[hi] ?? false) && !(feeds[hi] ?? true) {
        // 1) the cell chain head -> mouth
        cells_c: vector<integer> = [];
        cells_r: vector<integer> = [];
        wc = hc;
        wr = hr;
        walking = true;
        while walking {
          cells_c += [wc];
          cells_r += [wr];
          bi = t.tr_flow[wr * nx + wc] ?? -1;
          if bi < 0 {
            walking = false;        // lake/pit mouth
          } else {
            (nc, nr) = hex_neighbor(wc, wr, bi);
            if nc < 0 || nc >= nx || nr < 0 || nr >= ny {
              walking = false;
            } else {
              if (t.tr_h[nr * nx + nc] ?? 0.0) <= p.tp_sea {
                cells_c += [nc];     // one step into the sea
                cells_r += [nr];
                walking = false;
              } else {
                wc = nc;
                wr = nr;
              }
            }
          }
        }
        if len(cells_c) >= 3 {
          // 2) anchor polyline: control pt -> edge crossing -> control pt ...
          axs: vector<float> = [];
          ays: vector<float> = [];
          avs: vector<float> = [];
          for k in 0..len(cells_c) {
            kc = cells_c[k] ?? 0;
            kr = cells_r[k] ?? 0;
            if k > 0 {
              pc = cells_c[k - 1] ?? 0;
              pr = cells_r[k - 1] ?? 0;
              (ex, ey) = terrain_edge_cross(p, pc, pr, kc, kr);
              axs += [ex];
              ays += [ey];
              avs += [t.tr_acc[pr * nx + pc] ?? 1.0];
            }
            (px2, py2) = terrain_control_pt(p, kc, kr);
            axs += [px2];
            ays += [py2];
            avs += [t.tr_acc[kr * nx + kc] ?? 1.0];
          }
          // 3) fractal displacement over the anchors
          fdx: vector<float> = [];
          fdy: vector<float> = [];
          fdv: vector<float> = [];
          for k in 0..len(axs) - 1 {
            th_subdiv(t, p, axs[k] ?? 0.0, ays[k] ?? 0.0, avs[k] ?? 1.0,
                      axs[k + 1] ?? 0.0, ays[k + 1] ?? 0.0, avs[k + 1] ?? 1.0, 0, fdx, fdy, fdv);
          }
          fdx += [axs[len(axs) - 1] ?? 0.0];
          fdy += [ays[len(ays) - 1] ?? 0.0];
          fdv += [avs[len(avs) - 1] ?? 1.0];
          // 4) one corner-cut smoothing pass (kills residual sharp corners)
          smx: vector<float> = [];
          smy: vector<float> = [];
          smv: vector<float> = [];
          smx += [fdx[0] ?? 0.0];
          smy += [fdy[0] ?? 0.0];
          smv += [fdv[0] ?? 1.0];
          for k in 0..len(fdx) - 1 {
            ka = fdx[k] ?? 0.0;
            kb = fdx[k + 1] ?? 0.0;
            la = fdy[k] ?? 0.0;
            lb = fdy[k + 1] ?? 0.0;
            va = fdv[k] ?? 1.0;
            vb = fdv[k + 1] ?? 1.0;
            smx += [ka * 0.75 + kb * 0.25];
            smy += [la * 0.75 + lb * 0.25];
            smv += [va * 0.75 + vb * 0.25];
            smx += [ka * 0.25 + kb * 0.75];
            smy += [la * 0.25 + lb * 0.75];
            smv += [va * 0.25 + vb * 0.75];
          }
          smx += [fdx[len(fdx) - 1] ?? 0.0];
          smy += [fdy[len(fdy) - 1] ?? 0.0];
          smv += [fdv[len(fdv) - 1] ?? 1.0];
          // 5) trim to the shorelines: split into land pieces whose ends lie
          //    exactly ON the sea/lake contour (bisection per crossing)
          px3: vector<float> = [];
          py3: vector<float> = [];
          pv3: vector<float> = [];
          for k in 0..len(smx) {
            kx = smx[k] ?? 0.0;
            ky = smy[k] ?? 0.0;
            kv = smv[k] ?? 1.0;
            wetk = terrain_water_at(t, types, p, kx, ky);
            if !wetk {
              if len(px3) == 0 && k > 0 {
                (cx2, cy2) = th_shore_cross(t, types, p, kx, ky,
                                            smx[k - 1] ?? 0.0, smy[k - 1] ?? 0.0);
                px3 += [cx2];
                py3 += [cy2];
                pv3 += [kv];
              }
              px3 += [kx];
              py3 += [ky];
              pv3 += [kv];
            } else {
              if len(px3) > 0 {
                (qx2, qy2) = th_shore_cross(t, types, p, smx[k - 1] ?? 0.0, smy[k - 1] ?? 0.0, kx, ky);
                px3 += [qx2];
                py3 += [qy2];
                pv3 += [kv];
                if len(px3) >= 3 {
                  out += [th_shape_river(t, p, px3, py3, pv3)];
                }
                px3 = [];
                py3 = [];
                pv3 = [];
              }
            }
          }
          if len(px3) >= 3 {
            out += [th_shape_river(t, p, px3, py3, pv3)];
          }
        }
      }
    }
  }
  out
}

// ── Fine sampling (the game-facing surface) ────────────────────────────────

// The terrain at one point, as `terrain_sample` answers it: the carved height, the land
// material, the water there (0 dry, 1 river, 2 sea, 3 lake) and the distance to the nearest
// river.
pub struct TerrainSample {
  pub ts_h:     float,    // carved fine height
  pub ts_mat:   integer,  // winning land material
  pub ts_water: integer,  // 0 dry, 1 river, 2 sea, 3 lake
  ts_rdist: float,    // distance to the nearest river course (1e9 = none)
}

// The full sample: surface + ADDITIVE valley carving from every course near
// the point + the water classification (sea floods the carved mouth).  Pure
// function of (t, types, p, rivers, x, y) — THE window-independence invariant.
//
// Example: @HXT-007
pub fn terrain_sample(t: Terrain, types: vector<TerrainType>, p: TerrainParams, rivers: vector<TerrainRiver>, x: float, y: float) -> TerrainSample {
  sf = terrain_surface_at(t, types, p, x, y);
  carve = 0.0;
  gdist = 1.0e9;
  gw = 0.0;
  for ri in 0..len(rivers) {
    rv = rivers[ri];
    if rv != null && x >= rv.rv_x0 && x <= rv.rv_x1 && y >= rv.rv_y0 && y <= rv.rv_y1 {
      cd = 1.0e9;
      cdep = 0.0;
      cvw = p.tp_valley_w;
      for k in 0..len(rv.rv_x) - 1 {
        ax = rv.rv_x[k] ?? 0.0;
        ay = rv.rv_y[k] ?? 0.0;
        bx = rv.rv_x[k + 1] ?? 0.0;
        by = rv.rv_y[k + 1] ?? 0.0;
        ex = bx - ax;
        ey = by - ay;
        ee = ex * ex + ey * ey;
        if ee > 0.000000001 {
          tt = min(max((((x - ax) * ex + (y - ay) * ey) / ee) ?? 0.0, 0.0), 1.0);
          ddx = x - (ax + tt * ex);
          ddy = y - (ay + tt * ey);
          d = sqrt(ddx * ddx + ddy * ddy) ?? 0.0;
          if d < cd {
            cd = d;
            cdep = (rv.rv_d[k] ?? 0.0) * (1.0 - tt) + (rv.rv_d[k + 1] ?? 0.0) * tt;
            cvw = (rv.rv_vw[k] ?? p.tp_valley_w) * (1.0 - tt) + (rv.rv_vw[k + 1] ?? p.tp_valley_w) * tt;
            if d < gdist {
              gdist = d;
              gw = (rv.rv_w[k] ?? 0.0) * (1.0 - tt) + (rv.rv_w[k + 1] ?? 0.0) * tt;
            }
          }
        }
      }
      if cd < cvw * 4.0 {
        carve = carve + cdep * exp(0.0 - ((cd / cvw) ?? 0.0)) * (1.0 - ((cd / (4.0 * cvw)) ?? 0.0));
      }
    }
  }
  h = sf.sf_h - min(carve, p.tp_carve_max * 1.25);
  water = 0;
  if h <= p.tp_sea {
    water = 2;     // sea — including the flooded estuary
  } else {
    if terrain_lake_field(p, x, y, sf.sf_lakey) >= 0.5 {
      water = 3;   // lake (a FLAT filled surface — render without hillshade)
    } else {
      if gdist < gw * 0.5 && sf.sf_h > p.tp_sea && sf.sf_watery < 0.5 {
        water = 1; // river
      }
    }
  }
  TerrainSample { ts_h: h, ts_mat: sf.sf_mat, ts_water: water, ts_rdist: gdist }
}

tests/01-hex-terrain.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
// hex_terrain: hydrology, relief compounding, river attachment + meander,
// carving, and THE window-independence invariant (samples re-derived from a
// rebuilt world are identical).
use hex_terrain::*;
use hex_grid;

// Indices into the test's type table.
SEA = 0;
LAKE = 1;
PLAIN = 2;
HILL = 3;
ROCK = 4;

fn test_types() -> vector<hex_terrain::TerrainType> {
  ts: vector<hex_terrain::TerrainType> = [];
  ts += [hex_terrain::TerrainType { tt_key: "sea", tt_rise: 0.0, tt_steep: 1.5, tt_wet: true }];
  ts += [hex_terrain::TerrainType { tt_key: "lake", tt_rise: 0.0, tt_steep: 2.0, tt_wet: true }];
  ts += [hex_terrain::TerrainType { tt_key: "plain", tt_rise: 10.0, tt_steep: 2.5, tt_wet: false }];
  ts += [hex_terrain::TerrainType { tt_key: "hill", tt_rise: 80.0, tt_steep: 50.0, tt_wet: false }];
  ts += [hex_terrain::TerrainType { tt_key: "rock", tt_rise: 140.0, tt_steep: 98.0, tt_wet: false }];
  ts
}

// A small deterministic island: radial falloff + fbm + one massif, sea rim.
// The continent recipe is CALLER side (content); the lib provides mechanism.
fn build_world(seed: integer) -> (hex_terrain::Terrain, hex_terrain::TerrainParams) {
  p = hex_terrain::terrain_params(seed, 400.0);
  p.tp_acc_min = 3.0;
  nx = 16;
  ny = 13;
  t = hex_terrain::terrain_new(nx, ny);
  th = &t.tr_h;
  tmo = &t.tr_moist;
  ww = (nx as float) * 400.0;
  wh = ((ny as float) * 1.5 * 400.0 / hex_grid::SQRT3) ?? 0.0;
  for r in 0..ny {
    for c in 0..nx {
      (x, y) = hex_terrain::terrain_center(p, c, r);
      n = hex_terrain::terrain_fbm(seed, x, y, 3200.0, 4, 7);
      dx = ((x - ww * 0.45) / (ww * 0.42)) ?? 0.0;
      dy = ((y - wh * 0.5) / (wh * 0.46)) ?? 0.0;
      rad = sqrt(dx * dx + dy * dy);
      gx = ((x - ww * 0.35) / (ww * 0.2)) ?? 0.0;
      gy = ((y - wh * 0.42) / (wh * 0.3)) ?? 0.0;
      g = exp(0.0 - gx * gx - gy * gy);
      i = hex_terrain::terrain_idx(t, c, r);
      if i < len(th) {
        th[i] = 170.0 * n + 90.0 - 480.0 * (pow(rad, 2.4) ?? 0.0) + 520.0 * g;
      }
      if i < len(tmo) {
        tmo[i] = 0.5;
      }
    }
  }
  hex_terrain::terrain_hydrology(t, p, LAKE);
  // classify land (caller content): sea by height, lakes already set
  tmat = &t.tr_mat;
  for ci in 0..nx * ny {
    if ci < len(tmat) {
      if (th[ci] ?? 0.0) <= p.tp_sea {
        tmat[ci] = SEA;
      } else {
        if (tmat[ci] ?? 0) != LAKE {
          if (th[ci] ?? 0.0) > 320.0 {
            tmat[ci] = ROCK;
          } else {
            if (th[ci] ?? 0.0) > 160.0 {
              tmat[ci] = HILL;
            } else {
              tmat[ci] = PLAIN;
            }
          }
        }
      }
    }
  }
  hex_terrain::terrain_relief_pass(t, test_types(), p);
  (t, p)
}

fn main() {
  println("=== hex_terrain package test ===");
  (t, p) = build_world(1337);
  types = test_types();

  // 1) HYDROLOGY: after the flood every land cell drains (flow >= 0) and
  //    following the flow strictly descends to the sea — no pits, no cycles.
  land = 0;
  drains = true;
  descends = true;
  for r in 0..t.tr_ny {
    for c in 0..t.tr_nx {
      i = hex_terrain::terrain_idx(t, c, r);
      edge = c == 0 || r == 0 || c == t.tr_nx - 1 || r == t.tr_ny - 1;
      if (t.tr_h[i] ?? 0.0) > p.tp_sea && !edge {
        land = land + 1;
        if (t.tr_flow[i] ?? -1) < 0 { drains = false; }
        if (t.tr_flow[i] ?? -1) >= 0 {
          (nc, nr) = hex_terrain::terrain_flow_to(t, c, r);
          if (t.tr_h[nr * t.tr_nx + nc] ?? 0.0) >= (t.tr_h[i] ?? 0.0) { descends = false; }
        }
      }
    }
  }
  assert(land > 30, "the island has a real interior ({land} land cells)");
  assert(drains, "every interior land cell has a flow direction (pits filled)");
  assert(descends, "flow strictly descends on the filled surface");

  // 2) RELIEF: wet cells contribute zero; some land cell carries real relief.
  relmax = 0.0;
  wetzero = true;
  for wi in 0..t.tr_nx * t.tr_ny {
    if (t.tr_wet[wi] ?? false) && (t.tr_relief[wi] ?? 0.0) != 0.0 { wetzero = false; }
    relmax = max(relmax, t.tr_relief[wi] ?? 0.0);
  }
  assert(wetzero, "water-flow cells blend ZERO relief (they cut deep for free)");
  assert(relmax > 40.0, "the massif carries compounded relief ({relmax})");

  // 3) RIVERS exist, never run straight, and ATTACH to the shoreline.
  rivers = hex_terrain::terrain_rivers(t, types, p);
  assert(len(rivers) >= 1, "at least one river course ({len(rivers)})");
  crooked = true;
  attached = true;
  for ri in 0..len(rivers) {
    rv = rivers[ri];
    n = len(rv.rv_x);
    pathlen = 0.0;
    for k in 0..n - 1 {
      ddx = (rv.rv_x[k + 1] ?? 0.0) - (rv.rv_x[k] ?? 0.0);
      ddy = (rv.rv_y[k + 1] ?? 0.0) - (rv.rv_y[k] ?? 0.0);
      pathlen = pathlen + (sqrt(ddx * ddx + ddy * ddy) ?? 0.0);
    }
    sx = (rv.rv_x[n - 1] ?? 0.0) - (rv.rv_x[0] ?? 0.0);
    sy = (rv.rv_y[n - 1] ?? 0.0) - (rv.rv_y[0] ?? 0.0);
    straight = sqrt(sx * sx + sy * sy);
    if pathlen < straight * 1.02 { crooked = false; }
    // the LAST vertex sits ON the water line: just beyond it (extrapolate the
    // final segment a few meters) is water, while the vertex itself is land.
    ex = (rv.rv_x[n - 1] ?? 0.0) - (rv.rv_x[n - 2] ?? 0.0);
    ey = (rv.rv_y[n - 1] ?? 0.0) - (rv.rv_y[n - 2] ?? 0.0);
    el = sqrt(ex * ex + ey * ey);
    if el > 0.001 {
      bx = (rv.rv_x[n - 1] ?? 0.0) + ((ex / el) ?? 0.0) * 4.0;
      by = (rv.rv_y[n - 1] ?? 0.0) + ((ey / el) ?? 0.0) * 4.0;
      ax = (rv.rv_x[n - 1] ?? 0.0) - ((ex / el) ?? 0.0) * 4.0;
      ay = (rv.rv_y[n - 1] ?? 0.0) - ((ey / el) ?? 0.0) * 4.0;
      if !hex_terrain::terrain_water_at(t, types, p, bx, by) { attached = false; }
      if hex_terrain::terrain_water_at(t, types, p, ax, ay) { attached = false; }
    }
  }
  assert(crooked, "rivers never go straight (fractal course)");
  assert(attached, "every course end lies ON the shoreline (trim + bisection)");

  // 4) CARVING: at a mid-course vertex the carved height is below the
  //    pre-carve surface, and the vertex reads as water (river or flooded).
  rv0 = rivers[0];
  mid = (len(rv0.rv_x) / 2) ?? 0;
  mx = rv0.rv_x[mid] ?? 0.0;
  my = rv0.rv_y[mid] ?? 0.0;
  sf = hex_terrain::terrain_surface_at(t, types, p, mx, my);
  smp = hex_terrain::terrain_sample(t, types, p, rivers, mx, my);
  assert(smp.ts_h < sf.sf_h, "the course carves its valley ({smp.ts_h} < {sf.sf_h})");
  assert(smp.ts_water != 0, "the course centerline is water (kind {smp.ts_water})");

  // 5) RISE and STEEP are independent: with both zero the fine height equals
  //    the coarse blend exactly; steep alone adds detail, rise alone needs relief.
  flat: vector<hex_terrain::TerrainType> = [];
  flat += [hex_terrain::TerrainType { tt_key: "sea", tt_rise: 0.0, tt_steep: 1.5, tt_wet: true }];
  flat += [hex_terrain::TerrainType { tt_key: "lake", tt_rise: 0.0, tt_steep: 2.0, tt_wet: true }];
  flat += [hex_terrain::TerrainType { tt_key: "plain", tt_rise: 0.0, tt_steep: 0.0, tt_wet: false }];
  flat += [hex_terrain::TerrainType { tt_key: "hill", tt_rise: 0.0, tt_steep: 0.0, tt_wet: false }];
  flat += [hex_terrain::TerrainType { tt_key: "rock", tt_rise: 0.0, tt_steep: 0.0, tt_wet: false }];
  (qx, qy) = hex_terrain::terrain_center(p, 7, 6);
  sf0 = hex_terrain::terrain_surface_at(t, flat, p, qx + 13.0, qy + 7.0);
  assert(abs(sf0.sf_h - sf0.sf_blend) < 0.0000001,
         "rise 0 + steep 0 -> fine height == the coarse blend");
  sfr = hex_terrain::terrain_surface_at(t, types, p, qx + 13.0, qy + 7.0);
  assert(sfr.sf_h != sf0.sf_blend || sfr.sf_relief == 0.0,
         "the tuning numbers shape the surface");

  // 6) THE INVARIANT: rebuild the whole world from the same seed — samples at
  //    arbitrary points are IDENTICAL (pure functions all the way down).
  (t2, p2) = build_world(1337);
  rivers2 = hex_terrain::terrain_rivers(t2, types, p2);
  same = true;
  for k in 0..5 {
    px = 900.0 + (k as float) * 731.0;
    py = 700.0 + (k as float) * 577.0;
    a = hex_terrain::terrain_sample(t, types, p, rivers, px, py);
    b = hex_terrain::terrain_sample(t2, types, p2, rivers2, px, py);
    if a.ts_h != b.ts_h || a.ts_mat != b.ts_mat || a.ts_water != b.ts_water {
      same = false;
    }
  }
  assert(same, "window independence: re-derived samples are byte-identical");

  // 7) DIFFERENT SEED -> a different world (the seed actually steers).
  (t3, p3) = build_world(99);
  r3 = hex_terrain::terrain_rivers(t3, types, p3);
  diff = false;
  for k in 0..5 {
    px = 900.0 + (k as float) * 731.0;
    py = 700.0 + (k as float) * 577.0;
    a = hex_terrain::terrain_sample(t, types, p, rivers, px, py);
    b = hex_terrain::terrain_sample(t3, types, p3, r3, px, py);
    if a.ts_h != b.ts_h { diff = true; }
  }
  assert(diff, "a different seed builds a different world");

  println("=== HEX_TERRAIN OK ===");
}

tests/02-worked-examples.loft

// Copyright (c) 2026 Jurjen Stellingwerff
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// Worked examples for `hex_terrain` (@PLN141).  Each `test_*` below is a REAL,
// run-every-CI demonstration of a contract the signature cannot state; the
// function documents itself by CITING the tag with `// Example: @HXT-00N`, and
// `check_doc_drift.sh examples` fails if a citation goes dangling.  These are
// teaching examples, distinct from the conformance suite in `01-hex-terrain.loft`.
//
// WHAT THIS PACKAGE ASKS OF A CALLER, and why a signature cannot carry it.
// `hex_terrain` is the OVERLAND layer: one cell is hundreds of metres and is the
// terrain AUTHORITY, while the world you walk through is fine.  Everything hard
// about using it is that boundary:
//
//   the lattice   cells are `(c, r)` on a coarse hex grid — the THIRD lattice in
//                 this repo, and its integers mean nothing to the other two
//   the sampling  every `*_at` / `terrain_sample` call takes METRES, and there
//                 is no coordinate it will refuse
//   the passes    hydrology and relief are GLOBAL, run once over the whole map,
//                 and they rewrite the record the pure sampling reads
//
// The fixture below is a hand-authored ramp: `h(c, r) = 100c + 10r` with a sea
// column at `c = 0`, 500 m tiles, so every number in these tests is arithmetic
// you can redo on paper.
use hex_terrain::*;
use hex_grid::*;

TILE = 500.0;
SEED = 7;

// Material indices into the type table below.
SEA   = 0;
LAKE  = 1;
PLAIN = 2;
ROCK  = 3;

fn types4() -> vector<TerrainType> {
  ts: vector<TerrainType> = [];
  ts += [TerrainType { tt_key: "sea",   tt_rise: 0.0,   tt_steep: 0.0,  tt_wet: true }];
  ts += [TerrainType { tt_key: "lake",  tt_rise: 0.0,   tt_steep: 0.0,  tt_wet: true }];
  ts += [TerrainType { tt_key: "plain", tt_rise: 10.0,  tt_steep: 2.0,  tt_wet: false }];
  ts += [TerrainType { tt_key: "rock",  tt_rise: 140.0, tt_steep: 98.0, tt_wet: false }];
  ts
}

// The fixture: a ramp climbing east, a sea column at c = 0, optionally one PIT.
// The CONTENT is the caller's — this library supplies the mechanism only — so a
// terrain is authored by writing the vectors and then running the passes.
fn ramp(nx: integer, ny: integer, pit: boolean) -> Terrain {
  t = terrain_new(nx, ny);
  th = &t.tr_h;
  tm = &t.tr_mat;
  for r in 0..ny {
    for c in 0..nx {
      i = r * nx + c;
      if c == 0 {
        th[i] = 0.0 - 50.0;
        tm[i] = SEA;
      } else {
        th[i] = 100.0 * (c as float) + 10.0 * (r as float);
        tm[i] = PLAIN;
      }
    }
  }
  if pit { th[2 * nx + 2] = 5.0; }
  t
}

// The same fixture with both global passes run — what a caller actually samples.
fn cooked(nx: integer, ny: integer, pit: boolean) -> Terrain {
  t = ramp(nx, ny, pit);
  terrain_hydrology(t, terrain_params(SEED, TILE), LAKE);
  terrain_relief_pass(t, types4(), terrain_params(SEED, TILE));
  t
}

// The smallest call that does something real: author a ramp, run the passes, and
// sample a point in metres.
fn main() {
  p = terrain_params(SEED, TILE);
  t = cooked(5, 5, false);
  (x, y) = terrain_center(p, 2, 2);           // cell (2,2) in METRES
  s = terrain_surface_at(t, types4(), p, x, y);
  assert(s.sf_mat == PLAIN, "the middle of the ramp is plain land");
  assert(s.sf_h > 200.0 && s.sf_h < 240.0, "at about the cell's own 220 m, got {s.sf_h}");
}

// @HXT-001 — THREE LATTICES IN ONE REPO, AND ONLY ONE OF THEM IS METRES.
// `hex_grid` walks a fine lattice; `hex_world`'s axial coordinates are a second;
// this package's `(c, r)` is a third — a coarse OVERLAND grid whose pitch is
// `tp_tile` metres.  The integers look identical and mean nothing to each other,
// and `terrain_center` is the ONE bridge across.
//
// Two things bite.  First, "500 m tiles" does not make a 500 m square grid: the
// column pitch is exactly the tile, and the ROW pitch is `tile * sqrt(3)/2` with
// odd rows offset by half a tile — this is a hex lattice, not a raster.  Second,
// every sampling entry point takes METRES and will refuse nothing, so a
// coordinate that crossed the boundary — cell indices, or a walked position in
// world units — returns a real height from a real place that is not the place
// you meant.
fn test_the_overland_lattice_is_a_third_one_and_it_is_in_metres() {
  p = terrain_params(SEED, TILE);

  (x0, y0) = terrain_center(p, 0, 0);
  (x1, y1) = terrain_center(p, 1, 0);
  (x2, y2) = terrain_center(p, 0, 1);
  assert(x0 == 0.0 && y0 == 0.0, "cell (0,0) is the origin");
  assert(abs(x1 - TILE) < 1e-9 && y1 == 0.0, "the COLUMN pitch is exactly the tile, got {x1}");
  assert(abs(y2 - TILE * 0.8660254037844386) < 1e-9,
         "the ROW pitch is tile * sqrt(3)/2 = 433.01 m, not 500, got {y2}");
  assert(abs(x2 - TILE * 0.5) < 1e-9, "and an odd row is offset half a tile, got {x2}");

  // The same integers in the WALKED lattice are a different place entirely.
  (cx, cy) = terrain_center(p, 3, 2);
  (wx, wy) = hex_to_px(3, 2);
  assert(abs(cx - 1500.0) < 1e-9 && abs(cy - 866.0254037844386) < 1e-9,
         "overland (3,2) is 1500 m east and 866 m north, got ({cx},{cy})");
  assert(wx < 6.0 && wy < 4.0, "hex_grid's own (3,2) is ({wx},{wy}) — its own units, its own lattice");

  // What a crossed coordinate costs: a height, from somewhere else, with no complaint.
  t = cooked(5, 5, false);
  here = terrain_blend_h(t, p, cx, cy);
  as_index = terrain_blend_h(t, p, 3.0, 2.0);
  as_walked = terrain_blend_h(t, p, wx, wy);
  assert(here > 300.0, "the cell's own centre reads about 320 m, got {here}");
  assert(as_index < 0.0, "the INDEX read as metres lands in the sea column, got {as_index}");
  assert(as_walked < 0.0, "and so does a walked position, got {as_walked}");
  assert(as_index != as_walked, "two different wrong places, both perfectly plausible");
}

// @HXT-002 — "WINDOW INDEPENDENCE" IS DETERMINISM, NOT EXTENT-INDEPENDENCE, AND
// THE DIFFERENCE DECIDES HOW YOU GENERATE A WORLD.  What holds is that sampling
// is a PURE function of the terrain record: the same point answers the same
// value whatever else you sampled before it, in any order, from any viewport.
// That is what makes the fine world derivable on demand and it is worked first
// below.
//
// What does NOT hold is that the record is independent of the map's EXTENT.
// Hydrology is a global pass and flow ACCUMULATION is a catchment property — it
// counts every cell upstream — so widening the map raises accumulation
// everywhere, which flips cells to water, which zeroes their relief, which moves
// the fine surface through the blend kernel.  The same authored content in a
// 5x5 and an 8x8 window is not the same terrain.
//
// So: generate the overland map ONCE, at its full extent.  A window is a unit of
// SAMPLING, never a unit of generation.
fn test_sampling_is_pure_but_the_record_depends_on_the_whole_map() {
  p = terrain_params(SEED, TILE);
  ty = types4();
  big = cooked(8, 8, false);
  (qx, qy) = terrain_center(p, 2, 2);

  // The half that holds: order and viewport independence.
  first = terrain_surface_at(big, ty, p, qx + 41.0, qy + 17.0);
  elsewhere = terrain_surface_at(big, ty, p, 12345.0, 6789.0);
  again = terrain_surface_at(big, ty, p, qx + 41.0, qy + 17.0);
  assert(first.sf_h == again.sf_h, "the same point answers the same value, exactly");
  assert(elsewhere.sf_h != first.sf_h, "and an unrelated sample in between changes nothing");

  // The half that does not: the same content, two extents.
  small = cooked(5, 5, false);
  ia = terrain_idx(small, 2, 1);
  ib = terrain_idx(big, 2, 1);
  assert((small.tr_acc[ia] ?? 0.0) == 5.0 && (big.tr_acc[ib] ?? 0.0) == 11.0,
         "cell (2,1) drains 5 tiles in the small map and 11 in the big one");
  assert(!(small.tr_wet[ia] ?? true) && (big.tr_wet[ib] ?? false),
         "which crosses tp_acc_min = 6, so the big map calls that cell water");
  assert((small.tr_relief[ia] ?? 0.0) == 154.0 && (big.tr_relief[ib] ?? 0.0) == 0.0,
         "and a water cell contributes ZERO relief");

  // …and that reaches the fine surface, well inside both windows.
  worst = 0.0;
  ndiff = 0;
  nsame = 0;
  for r in 1..4 {
    for c in 1..4 {
      (bx, by) = terrain_center(p, c, r);
      for k in 0..3 {
        px = bx + (k as float) * 90.0;
        py = by + (k as float) * 60.0;
        sa = terrain_surface_at(small, ty, p, px, py);
        sb = terrain_surface_at(big, ty, p, px, py);
        d = abs(sa.sf_h - sb.sf_h);
        if d > worst { worst = d; }
        if sa.sf_h == sb.sf_h { nsame = nsame + 1; } else { ndiff = ndiff + 1; }
      }
    }
  }
  assert(nsame + ndiff == 27, "27 interior samples");
  assert(ndiff == 24, "24 of them differ between the two extents, got {ndiff}");
  assert(abs(worst - 15.103122697077197) < 1e-9, "by up to 15.1 m, got {worst}");

  // The number that makes it concrete: the small map has no rivers at all, and
  // the big one has five — same content, same seed, same parameters.
  assert(len(terrain_rivers(small, ty, p)) == 0, "no course reaches tp_acc_min in the 5x5");
  assert(len(terrain_rivers(big, ty, p)) == 5, "five do in the 8x8");
}

// @HXT-003 — THE PASSES HAVE AN ORDER, AND HYDROLOGY REWRITES THE HEIGHTS YOU
// AUTHORED.  `terrain_hydrology` is a priority-flood: it raises every pit until
// the whole map drains, and it writes that FILLED surface back into `tr_h`.  The
// height you wrote is gone, and nothing in the record says a pass has run — no
// flag, no version, no second vector.
//
// The order is hydrology, then the caller's land classification, then
// `terrain_relief_pass`, then `terrain_rivers`.  Relief reads the wetness
// hydrology computed; rivers read the flow and accumulation it computed.  Run
// them out of order and every one of them answers plausibly from stale data.
fn test_hydrology_rewrites_the_heights_and_the_passes_have_an_order() {
  p = terrain_params(SEED, TILE);
  t = ramp(5, 5, true);
  i = terrain_idx(t, 2, 2);
  assert((t.tr_h[i] ?? 0.0) == 5.0, "the pit is authored at 5 m");
  assert((t.tr_flow[i] ?? 0) == -1, "and before hydrology nothing has a flow direction");
  assert((t.tr_mat[i] ?? 0) == PLAIN, "nor is it water");

  terrain_hydrology(t, p, LAKE);

  // The pit's lowest rim is cell (1,1) at 110 m, and the flood adds its epsilon.
  assert(abs((t.tr_h[i] ?? 0.0) - 110.05) < 1e-9,
         "the pit comes back at its lowest rim + 0.05, got {t.tr_h[i] ?? 0.0}");
  assert((t.tr_mat[i] ?? 0) == LAKE, "and a raise that large makes it a lake");
  assert((t.tr_flow[i] ?? -1) >= 0, "everything drains now");

  // Every land cell drains, and following the flow strictly descends.
  descends = true;
  for r in 1..4 {
    for c in 1..4 {
      j = terrain_idx(t, c, r);
      (nc, nr) = terrain_flow_to(t, c, r);
      if (t.tr_h[terrain_idx(t, nc, nr)] ?? 0.0) >= (t.tr_h[j] ?? 0.0) { descends = false; }
    }
  }
  assert(descends, "flow strictly descends on the FILLED surface — no pits, no cycles");
}

// @HXT-004 — RELIEF IS COMPOUNDED BY THE COUNT OF STEEP NEIGHBOURS, AND WATER
// CONTRIBUTES NOTHING.  A cell's relief is its largest height difference to a
// neighbour, multiplied by `tp_nb_base + tp_nb_boost * (steep neighbours)` — so
// a cell in a busy neighbourhood gets more relief than the same drop in a quiet
// one, which is what makes ridges gather instead of speckling.
//
// A water-flow cell takes ZERO relief, and that is the rivers-cut-deep-for-free
// rule: the dry land rises around a course while the course itself stays put.
fn test_relief_compounds_by_the_number_of_steep_neighbours() {
  p = terrain_params(SEED, TILE);
  t = cooked(5, 5, true);
  assert(p.tp_nb_base == 0.6 && p.tp_nb_boost == 0.2 && p.tp_steep_delta == 60.0,
         "the defaults this arithmetic uses");

  // Cell (1,0) = 100 m.  Its neighbours are the sea at -50 (twice, a 150 m drop,
  // both steep), 200 to the east (100 m, steep) and 110 to the south-east (10 m,
  // not steep): max 150, three steep -> 150 * (0.6 + 3 * 0.2) = 180.
  a = terrain_idx(t, 1, 0);
  assert(abs((t.tr_relief[a] ?? 0.0) - 180.0) < 1e-9,
         "150 * (0.6 + 3*0.2) = 180, got {t.tr_relief[a] ?? 0.0}");

  // Cell (3,2) = 320 m, beside the filled lake at 110.05: max 209.95 with four
  // steep neighbours -> 209.95 * (0.6 + 4 * 0.2) = 293.93.
  b = terrain_idx(t, 3, 2);
  assert(abs((t.tr_relief[b] ?? 0.0) - 293.93) < 1e-9,
         "209.95 * 1.4 = 293.93, got {t.tr_relief[b] ?? 0.0}");

  // The lake itself, and the sea column, take none.
  lake = terrain_idx(t, 2, 2);
  sea = terrain_idx(t, 0, 2);
  assert((t.tr_wet[lake] ?? false) && (t.tr_relief[lake] ?? 1.0) == 0.0, "the lake takes no relief");
  assert((t.tr_wet[sea] ?? false) && (t.tr_relief[sea] ?? 1.0) == 0.0, "nor does the sea");
}

// @HXT-005 — SAMPLING A CELL'S OWN CENTRE DOES NOT RETURN THAT CELL'S HEIGHT.
// The blend kernel has radius `1.02 * tile`, so it always overlaps the six
// neighbours — that overlap is exactly what makes two adjacent tiles merge with
// no seam and no stitching anywhere.  The price is that a cell's authored value
// is diluted toward its neighbours, at its own centre, by about 0.9%.
//
// It is not an error term to be tuned away: it is the same weight that carries
// the merge.  Author a value expecting to read it back and you will be wrong by
// a fixed, computable fraction, and by more where the neighbours differ more.
fn test_a_cell_centre_does_not_read_the_cells_own_height() {
  p = terrain_params(SEED, TILE);
  t = terrain_new(5, 5);
  th = &t.tr_h;
  for r in 0..5 { for c in 0..5 { th[r * 5 + c] = 0.0; } }
  th[2 * 5 + 2] = 100.0;

  // The kernel weight of a neighbour at exactly one tile: (1 - 1/1.02^2)^2.
  w = 1.0 - (1.0 / (1.02 * 1.02) ?? 0.0);
  w = w * w;
  assert(abs(w - 0.001507863550543432) < 1e-15, "the neighbour weight, got {w}");

  (cx, cy) = terrain_center(p, 2, 2);
  own = terrain_blend_h(t, p, cx, cy);
  assert(abs(own - (100.0 / (1.0 + 6.0 * w) ?? 0.0)) < 1e-12,
         "its own centre reads 100/(1+6w) = 99.103, not 100 — got {own}");
  assert(own < 100.0 && own > 99.0, "a 0.9% dilution");

  // …and the neighbour picks up exactly what this cell lost.
  (nx1, ny1) = terrain_center(p, 3, 2);
  nb = terrain_blend_h(t, p, nx1, ny1);
  assert(abs(nb - (100.0 * w / (1.0 + 6.0 * w) ?? 0.0)) < 1e-12,
         "the neighbour's centre reads 100w/(1+6w) = 0.149, got {nb}");
  assert(abs((own + 6.0 * nb) - 100.0) < 1e-9, "the seven readings account for the whole 100 m");
}

// @HXT-006 — OUTSIDE THE MAP THERE IS NO EDGE, ONLY A PLATEAU.  The kernel
// clamps its cell indices, so a coordinate beyond the authored extent reads the
// nearest edge cell — forever, in every direction, with no complaint.  Ten
// kilometres past a five-cell map you still get a height, and it is a perfectly
// ordinary number.
//
// So bound your own queries.  The library has no concept of "outside", which is
// the same reason a crossed coordinate (@HXT-001) is silent: there is no
// coordinate this API can refuse.
fn test_outside_the_map_is_a_plateau_and_never_a_refusal() {
  p = terrain_params(SEED, TILE);
  t = cooked(5, 5, false);
  edge = terrain_idx(t, 4, 2);
  (ex, ey) = terrain_center(p, 4, 2);
  assert((t.tr_h[edge] ?? 0.0) == 420.0, "the eastern edge cell is 420 m");

  inside = terrain_blend_h(t, p, ex, ey);
  assert(inside < 420.0, "at its own centre the blend is still diluted, got {inside}");
  for k in [1.0, 2.0, 4.0, 20.0] {
    out = terrain_blend_h(t, p, ex + k * TILE, ey);
    assert(abs(out - 420.0) < 1e-9, "{k} tiles past the edge still reads 420 m, got {out}");
  }
  // …flatter outside than inside: every cell the kernel finds is the same clamped one.
  far = terrain_blend_h(t, p, ex + 20.0 * TILE, ey);
  assert(far > inside, "the plateau is the edge cell UNDILUTED, {far} against {inside}");
}

// @HXT-007 — THERE ARE TWO WATER ANSWERS AND THEY DISAGREE ON PURPOSE.
// `terrain_water_at` is the PRE-CARVE shoreline: the coast and lakes as the
// blended surface has them, which is the line river courses are trimmed to
// attach to.  `terrain_sample(...).ts_water` is the finished world: it includes
// the river itself and the mouth its carve floods.
//
// At a point on a river both are correct and they answer differently.  Ask the
// first when you need the shoreline, the second when you need to know whether
// something can stand there.
fn test_the_shoreline_and_the_finished_water_are_two_questions() {
  p = terrain_params(SEED, TILE);
  ty = types4();
  t = cooked(8, 8, false);
  rivers = terrain_rivers(t, ty, p);
  assert(len(rivers) > 0, "the 8x8 map carries courses");

  r0 = rivers[0];
  n = len(r0.rv_x);
  mid = (n / 2) ?? 0;
  mx = r0.rv_x[mid] ?? 0.0;
  my = r0.rv_y[mid] ?? 0.0;

  sf = terrain_surface_at(t, ty, p, mx, my);
  sm = terrain_sample(t, ty, p, rivers, mx, my);
  assert(sm.ts_h < sf.sf_h, "the course carves its valley: {sm.ts_h} below {sf.sf_h}");
  assert(!terrain_water_at(t, ty, p, mx, my), "the PRE-carve shoreline says this is land");
  assert(sm.ts_water != 0, "and the finished sample says it is water, kind {sm.ts_water}");
}

// @HXT-008 — A TYPE BOUNDARY IS A LINE, NOT A CROSSFADE, AND THE STEEPER TYPE
// WINS IT.  Two terrain types meeting do not blend over a wide band: the winner
// is decided per sample by `sqrt(steep)`-weighted proximity, so the transition
// is a crisp fractal line with a narrow softening zone around it.  `sf_mix`
// reports how decisively the winner won — 1.0 is well inside a type, 0.5 is the
// contested line itself.
fn test_a_type_boundary_is_a_line_and_the_steeper_type_wins() {
  p = terrain_params(SEED, TILE);
  ty = types4();
  t = cooked(8, 8, false);
  tm = &t.tr_mat;
  for r in 0..8 { for c in 3..8 { if c > 0 { tm[r * 8 + c] = ROCK; } } }

  (ax, ay) = terrain_center(p, 2, 2);
  (bx, _by) = terrain_center(p, 3, 2);
  flips = 0;
  soft = 0;
  prev = 0 - 1;
  for k in 0..41 {
    x = ax + ((bx - ax) * (k as float) / 40.0 ?? 0.0);
    s = terrain_surface_at(t, ty, p, x, ay);
    if prev >= 0 && s.sf_mat != prev { flips = flips + 1; }
    prev = s.sf_mat;
    if s.sf_mix < 0.999 { soft = soft + 1; }
    assert(s.sf_mix >= 0.5 && s.sf_mix <= 1.0, "mix stays in 0.5..1, got {s.sf_mix}");
  }
  assert(flips == 1, "41 samples across one tile cross the boundary exactly once, got {flips}");
  assert(soft <= 4, "and only {soft} of them are inside the softening zone — a line, not a fade");
}