using System;
using System.Collections.Generic;

/// <summary>
/// River bed carving — C0b part 2a (terrain-water task 22). The first river pass
/// that MODIFIES terrain: executes the frozen task-21b plan by carving channel
/// beds for the promoted rivers. NO WATER — part 2b puts water into these beds
/// once the routing-style gate picks SHORT or LOWGROUND.
///
/// Standalone numeric (D-035 family). Runs the task-21 DrainageAnalysis in-
/// pipeline (deterministic: same seed → same eroded surface → same plan), then:
///
///   1. LOWLAND ROUTING (the A/B): each routed giant gets a route from its
///      pooling terminal to the nearest OCEAN cell by deterministic Dijkstra.
///      SHORT     — cost ≈ distance, uphill penalised: heads direct, avoids walls.
///      LOWGROUND — cost ≈ elevation above sea: follows the lowest available
///                  ground and wanders like a real lowland river.
///   2. BED CARVING: every promoted course (trunks, routed giants + their lowland
///      reaches, lake-enders, tributaries) is stamped as a parabolic channel with
///      a smoothstep shoulder — a bed for water to sit in, not a canyon. Width and
///      depth grow downstream with drainage. The bed elevation is made MONOTONE
///      NON-INCREASING toward the outlet (water must flow), and is clamped to
///      sea + margin everywhere — the flood-guard discipline erosion established:
///      below-sea cells are read-only, no carve may create inland below-sea cells,
///      so the rendered coastline cannot move. The bed meets the ocean AT the
///      coast, where the terrain itself descends through sea level.
///   3. Crater: nothing inside the protected core is modified (erosion's rule).
///
/// Output-only: the caller applies this to the RENDER map after classify-side
/// data (biomes, WBID) is already computed — the oracle stays byte-identical.
/// </summary>
public static class RiverCarvePass
{
	public const float M_PER_UNIT = 251f;
	public const byte STYLE_SHORT = 0;
	public const byte STYLE_LOWGROUND = 1;

	private static readonly int[] DX = { -1, -1, -1, 0, 0, 1, 1, 1 };
	private static readonly int[] DY = { -1, 0, 1, -1, 1, -1, 0, 1 };
	private static readonly float[] DIST = {
		1.41421356f, 1f, 1.41421356f, 1f, 1f, 1.41421356f, 1f, 1.41421356f };

	// Routing cost constants. SHORT pays lightly for climbing (8 per metre of rise,
	// so a 10 m wall costs like an 80 px detour — walls are avoided, direction is
	// kept). LOWGROUND pays for BEING high (1 per metre of elevation per px) plus
	// heavily for climbing, so the cheapest corridor is the lowest ground even when
	// that wanders.
	private const float SHORT_UPHILL_PER_M = 8f;
	private const float LOWGROUND_ELEV_PER_M = 1f;
	private const float LOWGROUND_BASE = 0.05f;
	private const float LOWGROUND_UPHILL_PER_M = 50f;

	// Bed geometry: sizes grow downstream from head to mouth, scaled by
	// sqrt(drainage / 1e6) so a 2.3M px giant carves roughly 2.3× deeper/wider at
	// the mouth than a 0.4M px trunk. Kept channel-scale, per the erosion
	// detailing philosophy.
	private const float DEPTH_HEAD_M = 1.0f;
	private const float DEPTH_MOUTH_M = 6.0f;     // × sizeFactor × DepthScale
	private const float HALFWIDTH_HEAD_PX = 2.0f;
	private const float HALFWIDTH_MOUTH_PX = 12.0f; // × sizeFactor × WidthScale
	private const float MIN_BED_SLOPE = 0.002f;   // m per px of enforced descent

	public class Params
	{
		public byte RoutingStyle = STYLE_LOWGROUND;
		public float WidthScale = 1.0f;
		public float DepthScale = 1.0f;
		public float SeaMarginM = 0.2f;           // bed floor above sea, everywhere
		// Task 24: a tributary reach is WET where its along-course flow exceeds this;
		// upstream of that it TAPERS to dry over TribTaperPx rather than ending in a
		// wall of water. 0 = water tributaries end to end.
		public int TribWaterMinFlowPx = 40_000;
		public int TribTaperPx = 120;
		// Lake-enders route to the nearest water body of at least this size — the
		// nearest wet PIXEL was a puddle (task-23 gate finding).
		public int LakeMinTargetPx = 20_000;
		public DrainageAnalysis.Params PlanParams = new();
	}

	public class RiverStat
	{
		public string Name;
		public string Kind;
		public bool SouthernCandidate;
		public long DrainagePx;
		public int CourseLenPx;
		public int RouteLenPx;                    // lowland reach only (routed giants)
		public float RouteStraightPx;
		public float WanderRatio;                 // routeLen / straight-line
		public bool ReachedOcean;
		public float MaxCutM;
		public double VolumeM3;
	}

	/// <summary>The carved geometry the water stage consumes (main rivers only).</summary>
	public class CarvedRiver
	{
		public string Name, Kind;
		public bool Southern;
		public long DrainagePx;
		public List<(float x, float y)> Dense;    // head → mouth, ~1-px samples
		public float[] Bed;                       // raw units, monotone non-increasing
		public float[] HalfW;                     // px
		public bool ReachedWaterTerminal;         // lake-enders: extension reached classify water
		// Task 24: along-course flow (px of drainage) per sample, and the first index
		// that carries water. Between WetFrom-TaperPx and WetFrom the water tapers
		// (narrowing and shallowing to the bed) so a stream head fades out.
		public float[] Flow;
		public int WetFrom;
		public int TaperPx;
	}

	public class Stats
	{
		public List<RiverStat> Rivers = new();
		public List<CarvedRiver> Carved = new();  // for the stepped-water stage (task 23)
		internal float[] PrePass;                 // cumulative-cut baseline
		public long CarvedCells;
		public double CarvedVolumeM3;
		public float MaxCutM;                     // CUMULATIVE vs pre-pass heights (task-22 nit 1 fixed)
		public double AnalysisSeconds, RoutingSeconds, CarveSeconds;
	}

	/// <param name="isSignificantWater">Row-major mask of classify water belonging to
	/// bodies of at least LakeMinTargetPx cells (task 24). Lake-enders route to THIS,
	/// not to any wet pixel: the task-23 build routed one into a puddle a few hundred
	/// px short of the obvious lagoon, because "nearest classify water" is satisfied
	/// by a 3-cell pond.</param>
	public static Stats Apply(float[,] height, int mapSize, bool[] isOcean,
		bool[] isClassifyWater, bool[] isSignificantWater, float southX, float southY,
		float[,] seaMap, float seaFlat,
		float craterCx, float craterCy, float craterCoreRadius,
		Func<double> secondsNow, Params p)
	{
		int n = mapSize;
		var stats = new Stats();
		float SeaAt(int x, int y) => seaMap != null ? seaMap[x, y] : seaFlat;
		float coreSq = craterCoreRadius * craterCoreRadius;

		// --- the frozen plan, recomputed deterministically in-pipeline ---
		double t0 = secondsNow();
		var plan = DrainageAnalysis.Run(height, mapSize, isOcean, isClassifyWater, southX, southY, p.PlanParams);
		stats.AnalysisSeconds = secondsNow() - t0;

		// --- lowland routing for the routed giants (the A/B) ---
		t0 = secondsNow();
		var giantRoutes = new List<List<(float x, float y)>>();
		foreach (var g in plan.Giants)
		{
			if (g.Kind != "routed") { giantRoutes.Add(null); continue; }
			giantRoutes.Add(RouteToOcean(height, n, isOcean,
				(int)g.Terminal.x, (int)g.Terminal.y, p.RoutingStyle, SeaAt));
		}
		stats.RoutingSeconds = secondsNow() - t0;

		// --- carve ---
		t0 = secondsNow();
		// Pre-pass snapshot: max-cut is measured CUMULATIVELY against the heights
		// this pass found, not per-write — overlapping stamps re-cut a cell and the
		// per-write number understated the true deepest cut ~4× (task-22 nit 1).
		float[] pre = new float[n * n];
		for (int x = 0; x < n; x++)
			for (int y = 0; y < n; y++)
				pre[x * n + y] = height[x, y];
		stats.PrePass = pre;
		int riverIdx = 0;
		foreach (var t in plan.Trunks)
		{
			riverIdx++;
			var course = new List<(float x, float y)>(t.Course);
			course.Reverse();                                  // head → mouth
			var rs = CarveRiver($"trunk{riverIdx}", "ocean-trunk", t.DrainageAreaPx,
				course, null, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
			rs.ReachedOcean = true;                            // outlet is on the coast by construction
			foreach (var trib in t.Tributaries)
				CarveTributary(trib, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
		}
		int gi = 0;
		foreach (var g in plan.Giants)
		{
			var route = giantRoutes[gi]; gi++;
			var course = new List<(float x, float y)>(g.Course);
			course.Reverse();                                  // head → terminal

			// Lake-enders (task 23): the stem pools on dry ground short of its lake
			// BECAUSE its pooling point is a local minimum — a blind descent walk
			// dead-ends there immediately (measured: 0 steps). Route to the nearest
			// classify-water cell with the same lowground Dijkstra the routed giants
			// use, so the bed (and then the water) actually joins the lake.
			bool reachedLake = false;
			if (g.Kind == "lake-ender")
			{
				// Target SIGNIFICANT water (task 24). Fall back to any classify water
				// only if no significant body is reachable, so a seed whose lake-ender
				// genuinely has only small ponds still connects rather than dead-ending.
				var ext = RouteToOcean(height, n, isSignificantWater,
					(int)g.Terminal.x, (int)g.Terminal.y, STYLE_LOWGROUND, SeaAt);
				if (ext.Count == 0)
					ext = RouteToOcean(height, n, isClassifyWater,
						(int)g.Terminal.x, (int)g.Terminal.y, STYLE_LOWGROUND, SeaAt);
				if (ext.Count > 0) { route = ext; reachedLake = true; }
			}

			var rs = CarveRiver($"giant{gi}", g.Kind, g.DrainageAreaPx,
				course, route, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
			rs.SouthernCandidate = g.SouthernCandidate;
			rs.ReachedOcean = g.Kind == "routed" ? (route != null && route.Count > 0) : reachedLake;
			if (stats.Carved.Count > 0) stats.Carved[^1].ReachedWaterTerminal = reachedLake;
			foreach (var trib in g.Tributaries)
				CarveTributary(trib, height, n, SeaAt, coreSq, craterCx, craterCy, p, stats);
		}
		stats.CarveSeconds = secondsNow() - t0;
		return stats;
	}

	// ---- Route smoothing (task 23): the road pass's AAA pipeline, numerically ----
	// RDP(4.0) decimation + 4 Chaikin corner-cutting passes, endpoints pinned —
	// the same constants and structure as MapGenerator.SmoothPath, mirrored here
	// because this pass is Godot-free. Kills the 8-connected Dijkstra 45° kinks;
	// the bed then carves along the smoothed centreline.
	private static List<(float x, float y)> SmoothCourse(List<(float x, float y)> raw)
	{
		if (raw.Count < 3) return raw;
		var dec = Rdp(raw, 0, raw.Count - 1, 4.0f);
		if (dec.Count < 3) return raw;
		var sm = dec;
		for (int pass = 0; pass < 4; pass++)
		{
			var nxt = new List<(float x, float y)>(sm.Count * 2) { sm[0] };
			for (int i = 0; i + 1 < sm.Count; i++)
			{
				var a = sm[i]; var b = sm[i + 1];
				nxt.Add((a.x * 0.75f + b.x * 0.25f, a.y * 0.75f + b.y * 0.25f));
				nxt.Add((a.x * 0.25f + b.x * 0.75f, a.y * 0.25f + b.y * 0.75f));
			}
			nxt.Add(sm[^1]);
			sm = nxt;
		}
		return sm;
	}

	private static List<(float x, float y)> Rdp(List<(float x, float y)> pts, int i0, int i1, float tol)
	{
		if (i1 - i0 <= 1) return new List<(float x, float y)> { pts[i0], pts[i1] };
		var a = pts[i0]; var b = pts[i1];
		float abx = b.x - a.x, aby = b.y - a.y;
		float abLen = MathF.Sqrt(abx * abx + aby * aby);
		float maxD = 0f; int maxI = i0;
		for (int i = i0 + 1; i < i1; i++)
		{
			float d = abLen < 1e-6f
				? MathF.Sqrt((pts[i].x - a.x) * (pts[i].x - a.x) + (pts[i].y - a.y) * (pts[i].y - a.y))
				: MathF.Abs(abx * (a.y - pts[i].y) - (a.x - pts[i].x) * aby) / abLen;
			if (d > maxD) { maxD = d; maxI = i; }
		}
		if (maxD <= tol) return new List<(float x, float y)> { pts[i0], pts[i1] };
		var left = Rdp(pts, i0, maxI, tol);
		var right = Rdp(pts, maxI, i1, tol);
		left.RemoveAt(left.Count - 1);
		left.AddRange(right);
		return left;
	}

	/// <summary>
	/// Deterministic Dijkstra from the start cell to the nearest ocean cell under
	/// the selected style's cost model. Returns the path start → ocean (1-px steps),
	/// or an empty list if no path exists (reported upstream, never asserted away).
	/// </summary>
	private static List<(float x, float y)> RouteToOcean(float[,] height, int n,
		bool[] targets, int sx, int sy, byte style, Func<int, int, float> seaAt)
	{
		int total = n * n;
		var gcost = new float[total];
		var parent = new int[total];
		var closed = new bool[total];
		Array.Fill(gcost, float.MaxValue);
		Array.Fill(parent, -1);

		float ElevM(int x, int y) => MathF.Max(0f, (height[x, y] - seaAt(x, y)) * M_PER_UNIT);

		var pq = new PriorityQueue<int, (float c, int i)>();
		int start = sx * n + sy;
		gcost[start] = 0f;
		pq.Enqueue(start, (0f, start));
		int goal = -1;

		while (pq.Count > 0)
		{
			int c = pq.Dequeue();
			if (closed[c]) continue;
			closed[c] = true;
			if (targets[c]) { goal = c; break; }
			int cx = c / n, cy = c % n;
			float hc = height[cx, cy];
			for (int k = 0; k < 8; k++)
			{
				int nx = cx + DX[k], ny = cy + DY[k];
				if (nx < 0 || nx >= n || ny < 0 || ny >= n) continue;
				int ni = nx * n + ny;
				if (closed[ni]) continue;
				float dhM = MathF.Max(0f, (height[nx, ny] - hc) * M_PER_UNIT);
				float step = style == STYLE_SHORT
					? DIST[k] + dhM * SHORT_UPHILL_PER_M
					: DIST[k] * (LOWGROUND_BASE + ElevM(nx, ny) * LOWGROUND_ELEV_PER_M)
						+ dhM * LOWGROUND_UPHILL_PER_M;
				float nc = gcost[c] + step;
				if (nc < gcost[ni])
				{
					gcost[ni] = nc;
					parent[ni] = c;
					pq.Enqueue(ni, (nc, ni));
				}
			}
		}

		var path = new List<(float x, float y)>();
		if (goal >= 0)
		{
			for (int c = goal; c >= 0; c = parent[c])
				path.Add((c / n, c % n));
			path.Reverse();
		}
		return path;
	}

	private static void CarveTributary(DrainageAnalysis.Stream trib, float[,] height, int n,
		Func<int, int, float> seaAt, float coreSq, float craterCx, float craterCy,
		Params p, Stats stats)
	{
		var course = new List<(float x, float y)>(trib.Course);
		course.Reverse();                                      // head → confluence
		// Task 24: tributaries are registered as carved rivers (name "trib") so the
		// water stage sees them — they carved but stayed dry in task 23. They are NOT
		// added to stats.Rivers, so the per-river console table stays the 6 mains.
		CarveRiver("trib", "tributary", trib.DrainageAreaPx, course, null,
			height, n, seaAt, coreSq, craterCx, craterCy, p, stats);
	}

	/// <summary>
	/// Carves one river: densify the course, build a monotone-descending clamped
	/// bed profile, stamp the channel. Returns the per-river stat (also appended
	/// to stats.Rivers unless name is null — tributaries fold into the totals).
	/// </summary>
	private static RiverStat CarveRiver(string name, string kind, long drainagePx,
		List<(float x, float y)> upland, List<(float x, float y)> lowlandRoute,
		float[,] height, int n, Func<int, int, float> seaAt,
		float coreSq, float craterCx, float craterCy, Params p, Stats stats)
	{
		// Full head→mouth polyline: upland stem, then the lowland reach if any.
		// ONLY the lowland reach is smoothed: the Dijkstra 45° kinks live there, on
		// near-flat ground where a rounded corner costs nothing. The upland stems
		// already thread the carved valley FLOORS — smoothing them off-line cut
		// valley walls (measured: max cut 14.6 → 27.3 m before this was split).
		var pts = new List<(float x, float y)>(upland);
		if (lowlandRoute != null && lowlandRoute.Count > 1)
		{
			var smoothedRoute = SmoothCourse(lowlandRoute);
			pts.AddRange(smoothedRoute.GetRange(1, smoothedRoute.Count - 1));
		}

		// Densify to ~1-px samples (plan courses are decimated ×4).
		var dense = new List<(float x, float y)>();
		for (int i = 0; i + 1 < pts.Count; i++)
		{
			var a = pts[i]; var b = pts[i + 1];
			float segLen = MathF.Sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));
			int steps = Math.Max(1, (int)MathF.Ceiling(segLen));
			for (int s2 = 0; s2 < steps; s2++)
				dense.Add((a.x + (b.x - a.x) * s2 / steps, a.y + (b.y - a.y) * s2 / steps));
		}
		if (pts.Count > 0) dense.Add(pts[^1]);
		if (dense.Count < 2) return new RiverStat();

		float sizeFactor = MathF.Sqrt(drainagePx / 1_000_000f);
		var rs = new RiverStat
		{
			Name = name, Kind = kind, DrainagePx = drainagePx,
			CourseLenPx = dense.Count,
			RouteLenPx = lowlandRoute?.Count ?? 0
		};
		if (lowlandRoute != null && lowlandRoute.Count > 1)
		{
			var a = lowlandRoute[0]; var b = lowlandRoute[^1];
			rs.RouteStraightPx = MathF.Sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y));
			// Wander = POLYLINE length over straight-line — cell count undercounts
			// diagonal steps and can read below 1, which is geometrically impossible.
			float polyLen = 0f;
			for (int i = 1; i < lowlandRoute.Count; i++)
			{
				float sdx = lowlandRoute[i].x - lowlandRoute[i - 1].x;
				float sdy = lowlandRoute[i].y - lowlandRoute[i - 1].y;
				polyLen += MathF.Sqrt(sdx * sdx + sdy * sdy);
			}
			rs.RouteLenPx = (int)polyLen;
			rs.WanderRatio = rs.RouteStraightPx > 1f ? polyLen / rs.RouteStraightPx : 1f;
		}

		// Bed profile: raw = terrain − depth(t); then monotone non-increasing
		// downstream; then clamped to sea + margin. The clamp can flatten the tail
		// near the mouth — allowed: non-increasing is what water needs, and the
		// flood guard is absolute.
		int m = dense.Count;
		var bed = new float[m];
		var depth = new float[m];
		var halfW = new float[m];
		for (int i = 0; i < m; i++)
		{
			float t = m > 1 ? (float)i / (m - 1) : 1f;
			depth[i] = (DEPTH_HEAD_M + (DEPTH_MOUTH_M * sizeFactor - DEPTH_HEAD_M) * t) * p.DepthScale;
			if (depth[i] < 0.5f) depth[i] = 0.5f;
			halfW[i] = (HALFWIDTH_HEAD_PX + (HALFWIDTH_MOUTH_PX * sizeFactor - HALFWIDTH_HEAD_PX) * t) * p.WidthScale;
			if (halfW[i] < 1.5f) halfW[i] = 1.5f;
			int cx = (int)dense[i].x, cy = (int)dense[i].y;
			bed[i] = height[cx, cy] - depth[i] / M_PER_UNIT;
		}
		for (int i = 1; i < m; i++)
		{
			float maxAllowed = bed[i - 1] - MIN_BED_SLOPE / M_PER_UNIT;
			if (bed[i] > maxAllowed) bed[i] = maxAllowed;
		}
		for (int i = 0; i < m; i++)
		{
			int cx = (int)dense[i].x, cy = (int)dense[i].y;
			float floor = seaAt(cx, cy) + p.SeaMarginM / M_PER_UNIT;
			if (bed[i] < floor) bed[i] = floor;
		}

		// Stamp: parabolic channel to the rim, smoothstep shoulder back to terrain.
		for (int i = 0; i < m; i++)
		{
			float hw = halfW[i];
			float outer = hw * 2f;
			int cx0 = (int)MathF.Floor(dense[i].x - outer), cx1 = (int)MathF.Ceiling(dense[i].x + outer);
			int cy0 = (int)MathF.Floor(dense[i].y - outer), cy1 = (int)MathF.Ceiling(dense[i].y + outer);
			float rimH = bed[i] + depth[i] / M_PER_UNIT;
			for (int x = cx0; x <= cx1; x++)
			{
				if (x < 0 || x >= n) continue;
				for (int y = cy0; y <= cy1; y++)
				{
					if (y < 0 || y >= n) continue;
					float rx = x - dense[i].x, ry = y - dense[i].y;
					float r = MathF.Sqrt(rx * rx + ry * ry);
					if (r > outer) continue;
					float ddx = x - craterCx, ddy = y - craterCy;
					if (ddx * ddx + ddy * ddy < coreSq) continue;   // crater core protected
					float sea = seaAt(x, y);
					float old = height[x, y];
					if (old < sea) continue;                        // below-sea cells read-only
					float target;
					if (r <= hw)
					{
						float f = r / hw;
						target = bed[i] + (depth[i] / M_PER_UNIT) * f * f;
					}
					else
					{
						float f = (r - hw) / hw;                    // 0..1 across the shoulder
						f = f * f * (3f - 2f * f);                  // smoothstep
						target = rimH + (old - rimH) * f;
					}
					float floor = sea + p.SeaMarginM / M_PER_UNIT;
					if (target < floor) target = floor;
					if (target < old)
					{
						float cutM = (old - target) * M_PER_UNIT;
						// Cumulative depth vs the PRE-PASS surface (nit 1): the
						// honest "how deep did we cut here in total" number.
						float cumM = (stats.PrePass[x * n + y] - target) * M_PER_UNIT;
						height[x, y] = target;
						stats.CarvedCells++;
						stats.CarvedVolumeM3 += cutM;
						if (cumM > stats.MaxCutM) stats.MaxCutM = cumM;
						if (cumM > rs.MaxCutM) rs.MaxCutM = cumM;
						rs.VolumeM3 += cutM;
					}
				}
			}
		}

		if (name != null)
		{
			if (kind != "tributary") stats.Rivers.Add(rs);
			// Along-course flow: drainage grows from a headwater trickle to the full
			// figure at the mouth. Quadratic in t so the substantial lower half
			// dominates — a first-order stand-in for real accumulation, which the
			// plan only records per-river.
			var flow = new float[m];
			for (int i = 0; i < m; i++)
			{
				float t = m > 1 ? (float)i / (m - 1) : 1f;
				flow[i] = drainagePx * (0.05f + 0.95f * t * t);
			}
			int wetFrom = 0;
			if (kind == "tributary" && p.TribWaterMinFlowPx > 0)
			{
				wetFrom = m;                       // dry unless the threshold is met
				for (int i = 0; i < m; i++)
					if (flow[i] >= p.TribWaterMinFlowPx) { wetFrom = i; break; }
			}
			stats.Carved.Add(new CarvedRiver
			{
				Name = name, Kind = kind, DrainagePx = drainagePx,
				Dense = dense, Bed = bed, HalfW = halfW,
				Flow = flow, WetFrom = wetFrom, TaperPx = p.TribTaperPx
			});
		}
		return rs;
	}

	/// <summary>
	/// The stepped-water builder (task 23, part 2b): segments each carved main
	/// river into REACHES — flat water bodies stepping down the bed toward the
	/// outlet — and stamps their ids into the WBID grid. Reuses the existing
	/// levels-not-cells water model exactly: one body per reach, one flat level
	/// each; the writer derives WSRF from body levels as it always has. The step
	/// drops are the smoothing dial (smaller drop = more, finer steps); tilted
	/// water is the deferred model B and is NOT built here.
	///
	/// Emission is plain data (no engine types): the caller turns reaches into
	/// WBTB entries. Wet cells: inside the channel half-width, currently dry in
	/// WBID, at/above sea (below-sea cells belong to the ocean/crater-seam rule),
	/// bed below the reach level. Existing water bodies are never overwritten —
	/// a river MEETS a lake or the sea, it does not repaint them.
	/// </summary>
	public class Reach
	{
		public ushort Id;
		public string River;
		public float Level;                        // raw units
		public int PixelCount;
		public double Cx, Cy;                      // centroid accumulators → mean
	}

	public static List<Reach> AddSteppedWater(float[,] height, int mapSize,
		ushort[,] wbid, ushort firstId, List<CarvedRiver> rivers,
		float[,] seaMap, float seaFlat, float craterCx, float craterCy,
		float craterCoreRadius, float stepDropM, float waterDepthM)
	{
		int n = mapSize;
		float SeaAt(int x, int y) => seaMap != null ? seaMap[x, y] : seaFlat;
		float coreSq = craterCoreRadius * craterCoreRadius;
		var reaches = new List<Reach>();
		ushort nextId = firstId;

		foreach (var r in rivers)
		{
			int m = r.Dense.Count;
			if (m < 2) continue;
			// Task 24: start at the wet-from index (tributary threshold; 0 for mains),
			// but back up by the taper length so the transition is a FADE, not a wall.
			int i = Math.Max(0, r.WetFrom - r.TaperPx);
			if (i >= m) continue;                       // entirely below threshold: dry
			int taperStart = i, taperEnd = Math.Min(m - 1, r.WetFrom);
			float lastLevel = float.MaxValue;
			while (i < m)
			{
				// Reach spans from i while the bed stays within stepDropM of the
				// reach's starting bed; its flat level sits waterDepthM above that
				// start (deepening toward the next step — the pool behind a riffle).
				float startBed = r.Bed[i];
				float level = startBed + waterDepthM / M_PER_UNIT;
				if (level >= lastLevel)                    // enforce strict descent
					level = lastLevel - 0.01f / M_PER_UNIT;
				int j = i;
				while (j < m && r.Bed[j] > startBed - stepDropM / M_PER_UNIT) j++;

				var reach = new Reach { Id = nextId, River = r.Name, Level = level };
				for (int k2 = i; k2 < j; k2++)
				{
					// Taper: 0 at the dry end of the fade zone → 1 at full water. The
					// water narrows AND its surface drops to the bed, so a stream head
					// thins out and disappears instead of ending in a flat wall.
					float taper = 1f;
					if (taperEnd > taperStart && k2 < taperEnd)
						taper = (float)(k2 - taperStart) / (taperEnd - taperStart);
					if (taper <= 0.02f) continue;
					float hw = r.HalfW[k2] * (0.25f + 0.75f * taper);
					int x0 = (int)MathF.Floor(r.Dense[k2].x - hw), x1 = (int)MathF.Ceiling(r.Dense[k2].x + hw);
					int y0 = (int)MathF.Floor(r.Dense[k2].y - hw), y1 = (int)MathF.Ceiling(r.Dense[k2].y + hw);
					for (int x = x0; x <= x1; x++)
					{
						if (x < 0 || x >= n) continue;
						for (int y = y0; y <= y1; y++)
						{
							if (y < 0 || y >= n) continue;
							if (wbid[x, y] != 0) continue;          // never repaint existing water
							float rx = x - r.Dense[k2].x, ry = y - r.Dense[k2].y;
							if (rx * rx + ry * ry > hw * hw) continue;
							float ddx = x - craterCx, ddy = y - craterCy;
							if (ddx * ddx + ddy * ddy < coreSq) continue;
							float h = height[x, y];
							float sea = SeaAt(x, y);
							if (h < sea) continue;                  // ocean/seam territory
							// The tapered surface sits between the bed and the reach
							// level, so the wet strip shallows as it narrows.
							float localLevel = taper >= 1f
								? level
								: r.Bed[k2] + (level - r.Bed[k2]) * taper;
							if (h >= localLevel) continue;          // bank above the water line
							wbid[x, y] = nextId;
							reach.PixelCount++;
							reach.Cx += x; reach.Cy += y;
						}
					}
				}
				if (reach.PixelCount > 0)
				{
					reach.Cx /= reach.PixelCount; reach.Cy /= reach.PixelCount;
					reaches.Add(reach);
					nextId++;
					lastLevel = level;
				}
				i = j;
			}
		}
		return reaches;
	}
}