Landscape grading: platforms, swales, slope classes and ponding

/** A sloped platform's corner elevations: a plane through the centroid at z0 falling `pct` % towards `dir`. */
export function planeZ(pts, z0, fall) {
  const cx = pts.reduce((s, p) => s + p.x, 0) / pts.length, cy = pts.reduce((s, p) => s + p.y, 0) / pts.length;
  if (!fall) return pts.map((p) => ({ x: p.x, y: p.y, z: z0 }));
  const a = (fall.dir * Math.PI) / 180, g = fall.pct / 100;
  return pts.map((p) => ({ x: p.x, y: p.y, z: z0 - g * ((p.x - cx) * Math.cos(a) + (p.y - cy) * Math.sin(a)) }));
}

/** The mean ground elevation along an outline (vertices and edge midpoints), or null when it is off the surface. */
export function meanGround(surface, pts) {
  if (!surface) return null;
  let s = 0, n = 0;
  for (let i = 0; i < pts.length; i++) {
    const a = pts[i], b = pts[(i + 1) % pts.length];
    for (const q of [a, { x: (a.x + b.x) / 2, y: (a.y + b.y) / 2 }]) { const z = surface.elevationAt(q.x, q.y, 1e-6); if (z !== null) { s += z; n++; } }
  }
  return n ? s / n : null;
}

/** A swale as a gradedSurface patch: a flat-bottomed channel (or a V when bottom = 0) along the line, its sides cast to the ground. */
export function swalePatch(surface, sw, units) {
  const P = sw.pts, problems = [];
  let L = 0;
  const segs = [];
  for (let i = 0; i + 1 < P.length; i++) { const l = Math.hypot(P[i + 1].x - P[i].x, P[i + 1].y - P[i].y); if (l > 1e-9) segs.push({ a: P[i], b: P[i + 1], l, s0: L }); L += l; }
  if (!segs.length) return { problems: [t('ls.g.err.short')], region: null, points: [], breaklines: [], info: null };
  const step = Math.max(units === 'm' ? 0.5 : 1.5, Math.min(L / 200, Math.max(sw.bottom, units === 'm' ? 1 : 3)));
  const stations = [];
  for (let s = 0; s < L; s += step) stations.push(s);
  stations.push(L);
  for (const g of segs) if (g.s0 > 0) stations.push(g.s0);   // every bend gets a section
  stations.sort((a, b) => a - b);
  const at = (s) => { const g = segs.find((q) => s <= q.s0 + q.l + 1e-9) || segs[segs.length - 1]; const f = (s - g.s0) / g.l; return { x: g.a.x + f * (g.b.x - g.a.x), y: g.a.y + f * (g.b.y - g.a.y), ux: (g.b.x - g.a.x) / g.l, uy: (g.b.y - g.a.y) / g.l }; };
  const g0 = surface.elevationAt(P[0].x, P[0].y, 1e-6);
  if (g0 === null) return { problems: [t('ls.g.err.off')], region: null, points: [], breaklines: [], info: null };
  const left = [], right = [], bl = [], br = [], cl = [], rays = [];
  let last = Infinity;
  for (const s of stations) {
    const q = at(s), gz = surface.elevationAt(q.x, q.y, 1e-6);
    if (gz === null) { problems.push(t('ls.g.err.offat', { s: s.toFixed(1) })); continue; }
    let zb = sw.fall ? g0 - sw.depth - (sw.fall / 100) * s : gz - sw.depth;
    if (!sw.fall) zb = Math.min(zb, last);   // following the ground, the bottom never runs uphill
    last = zb;
    const nx = -q.uy, ny = q.ux, h = sw.bottom / 2;
    const pl = { x: q.x + nx * h, y: q.y + ny * h, z: zb }, pr = { x: q.x - nx * h, y: q.y - ny * h, z: zb };
    const rl = castSlope(surface, pl, nx, ny, zb, sw.side, sw.side), rr = castSlope(surface, pr, -nx, -ny, zb, sw.side, sw.side);
    left.push({ x: rl.x, y: rl.y, z: rl.z }); right.push({ x: rr.x, y: rr.y, z: rr.z });
    bl.push(pl); br.push(pr); cl.push({ x: q.x, y: q.y, z: zb });
    rays.push(rl, rr);
  }
  if (left.length < 2) return { problems: problems.length ? problems : [t('ls.g.err.off')], region: null, points: [], breaklines: [], info: null };
  const region = [...left, ...right.slice().reverse()];
  if (Math.abs(polyArea(region)) < 1e-9) return { problems: [t('ls.g.err.flat')], region: null, points: [], breaklines: [], info: null };
  const breaklines = [{ name: 'swale-left', pts: left }, { name: 'swale-right', pts: right }];
  if (sw.bottom > 1e-6) breaklines.push({ name: 'swale-bottom-left', pts: bl }, { name: 'swale-bottom-right', pts: br });
  else breaklines.push({ name: 'swale-bottom', pts: cl });
  breaklines.push({ name: 'swale-end', pts: [left[0], ...(sw.bottom > 1e-6 ? [bl[0], br[0]] : [cl[0]]), right[0]] });
  breaklines.push({ name: 'swale-end', pts: [left.at(-1), ...(sw.bottom > 1e-6 ? [bl.at(-1), br.at(-1)] : [cl.at(-1)]), right.at(-1)] });
  return { problems: [], region, points: [], breaklines, info: { kind: 'swale', bottom: cl, left, right, rays, length: L, fallAt: cl.length > 1 ? (cl[0].z - cl.at(-1).z) / L : 0, offNotes: problems } };
}

/** A grading area: inside its outline the surface is rebuilt from the edge (on the ground below it), the spots and contours. */
function areaPatch(surface, a) {
  const edge = [];
  for (let i = 0; i < a.pts.length; i++) {
    const p = a.pts[i], q = a.pts[(i + 1) % a.pts.length];
    const n = Math.max(1, Math.ceil(Math.hypot(q.x - p.x, q.y - p.y) / Math.max(1e-6, Math.sqrt(Math.abs(polyArea(a.pts))) / 20)));
    for (let k = 0; k < n; k++) { const x = p.x + ((q.x - p.x) * k) / n, y = p.y + ((q.y - p.y) * k) / n; const z = surface.elevationAt(x, y, 1e-6); if (z === null) return { problems: [t('ls.g.err.off')], region: null, points: [], breaklines: [], info: null }; edge.push({ x, y, z }); }
  }
  const inside = (p) => pointInPoly(p.x, p.y, a.pts);
  const spots = a.spots.filter(inside).map((s) => ({ x: s.x, y: s.y, z: s.z }));
  const contours = a.contours.map((c) => ({ name: 'contour ' + c.z, pts: c.pts.filter(inside).map((p) => ({ x: p.x, y: p.y, z: c.z })) })).filter((c) => c.pts.length >= 2);
  return { problems: [], region: a.pts.map((p) => ({ x: p.x, y: p.y })), points: spots, breaklines: [{ name: a.name, pts: [...edge, edge[0]] }, ...contours], info: { kind: 'area', spots: spots.length, contours: contours.length } };
}

/** Triangles by gradient: cls 0 = flatter than limits[0], k = steeper than limits[k-1] (limits as fractions, e.g. 0.05, 0.333). */
export function slopeClasses(surface, limits = [0.05, 1 / 3]) {
  const { X, Y, tris, nT, ox, oy } = surface, counts = new Array(limits.length + 1).fill(0), out = [], area = new Array(limits.length + 1).fill(0);
  for (let k = 0; k < nT; k++) {
    const s = surface.slopeOf(k).slope;
    let c = 0;
    while (c < limits.length && s > limits[c] + 1e-12) c++;
    const a = tris[3 * k], b = tris[3 * k + 1], d = tris[3 * k + 2];
    const ar = Math.abs((X[b] - X[a]) * (Y[d] - Y[a]) - (Y[b] - Y[a]) * (X[d] - X[a])) / 2;
    counts[c]++; area[c] += ar;
    out.push({ cls: c, slope: s, pts: [{ x: X[a] + ox, y: Y[a] + oy }, { x: X[b] + ox, y: Y[b] + oy }, { x: X[d] + ox, y: Y[d] + oy }] });
  }
  return { counts, area, tris: out };
}

/** Ponding: every vertex filled to the lowest level water can drain out at (priority-flood from the surface's outer edge,
 *  Barnes, Lehman & Mulla 2014); a pond is a connected group of vertices filled above their ground by more than `tol`. */
export function ponding(surface, { tol = 1e-3 } = {}) {
  const { X, Y, Z, tris, ox, oy } = surface, n = X.length;
  const used = surface.usedVertices();
  const adj = Array.from({ length: n }, () => []), edgeCount = new Map();
  const key = (a, b) => (a < b ? a * n + b : b * n + a);
  for (let k = 0; k < tris.length; k += 3) for (let e = 0; e < 3; e++) {
    const a = tris[k + e], b = tris[k + ((e + 1) % 3)];
    adj[a].push(b);
    const kk = key(a, b);
    edgeCount.set(kk, (edgeCount.get(kk) || 0) + 1);
  }
  const W = new Float64Array(n).fill(Infinity), done = new Uint8Array(n);
  // a binary heap of [level, vertex]
  const heap = [];
  const push = (lv, v) => { heap.push([lv, v]); let i = heap.length - 1; while (i > 0) { const p = (i - 1) >> 1; if (heap[p][0] <= heap[i][0]) break; [heap[p], heap[i]] = [heap[i], heap[p]]; i = p; } };
  const pop = () => { const top = heap[0], last = heap.pop(); if (heap.length) { heap[0] = last; let i = 0; for (;;) { const l = 2 * i + 1, r = l + 1; let m = i; if (l < heap.length && heap[l][0] < heap[m][0]) m = l; if (r < heap.length && heap[r][0] < heap[m][0]) m = r; if (m === i) break; [heap[m], heap[i]] = [heap[i], heap[m]]; i = m; } } return top; };
  for (const [kk, c] of edgeCount) if (c === 1) { const a = Math.floor(kk / n), b = kk % n; for (const v of [a, b]) if (!done[v]) { done[v] = 1; W[v] = Z[v]; push(Z[v], v); } }   // the outer edge drains
  while (heap.length) {
    const [lv, v] = pop();
    for (const u of adj[v]) if (!done[u] && used[u]) { done[u] = 1; W[u] = Math.max(Z[u], lv); push(W[u], u); }
  }
  const seen = new Uint8Array(n), ponds = [];
  const triArea = (a, b, c) => Math.abs((X[b] - X[a]) * (Y[c] - Y[a]) - (Y[b] - Y[a]) * (X[c] - X[a])) / 2;
  const wet = (v) => W[v] - Z[v] > tol;
  for (let v = 0; v < n; v++) {
    if (seen[v] || !used[v] || !wet(v)) continue;
    const verts = [], stack = [v];
    seen[v] = 1;
    while (stack.length) { const a = stack.pop(); verts.push(a); for (const u of adj[a]) if (!seen[u] && wet(u)) { seen[u] = 1; stack.push(u); } }
    const set = new Set(verts);
    let depth = 0, area = 0, volume = 0, lx = 0, ly = 0;
    for (const a of verts) { depth = Math.max(depth, W[a] - Z[a]); lx += X[a]; ly += Y[a]; }
    for (let k = 0; k < tris.length; k += 3) {   // triangles with all three corners wet: area and (mean depth × area) volume
      const a = tris[k], b = tris[k + 1], c = tris[k + 2];
      if (!set.has(a) && !set.has(b) && !set.has(c)) continue;
      const ar = triArea(a, b, c), d = ((W[a] - Z[a]) + (W[b] - Z[b]) + (W[c] - Z[c])) / 3;
      if (set.has(a) && set.has(b) && set.has(c)) { area += ar; volume += ar * Math.max(0, d); }
    }
    ponds.push({ depth, area, volume, spill: W[verts[0]], x: lx / verts.length + ox, y: ly / verts.length + oy, verts });
  }
  ponds.sort((a, b) => b.volume - a.volume || b.depth - a.depth);
  return { ponds, fill: W };
}