package fluvial

import "math"

// ClampToRepose enforces a maximum slope everywhere: no cell may stand above a neighbour by more than
// talus * distance. It returns the mean thickness removed, in metres.
//
// This replaces iterating thermal.Apply inside the solve, which could not do the job however many passes it
// was given (measured: 3, 10 and 40 passes all left the steepest land slope at 64 degrees against a 22 degree
// repose). The reason is structural rather than a bug. That routine moves half the excess downhill, so on a
// *uniform* over-steep slope every cell sheds exactly as much as it receives, the net change is zero, and the
// slope is a fixed point. It relaxes only where the downhill flux diverges — which is why it cuts a cone,
// whose contours converge, and why it cannot touch a planar hillside.
//
// So the constraint is imposed directly instead. This is the priority-flood mirrored: pop cells in ascending
// elevation, and lower any neighbour standing higher than the repose angle allows. Because a lowered cell is
// set to h[c] + talus*d, which is at or above the elevation just popped, the queue stays monotone and the
// bucket queue works unchanged. One pass, O(n) with the bucket queue, and the constraint holds globally when
// it returns.
//
// It is not mass-conserving: the material is removed rather than piled at the foot of the slope. That is the
// deliberate simplification, because in this landscape the foot of a hillslope is a channel and the channel
// exports the sediment anyway. The mean thickness removed is returned so a run can report it, and a run that
// removes a suspicious amount is saying its uplift and its repose angle disagree.
func (g *Grid) ClampToRepose(h []float32, talus float64) float64 {
	if talus <= 0 {
		return 0
	}
	n := g.W * g.H
	for i := range g.closed {
		g.closed[i] = false
	}
	// Pushed with a jittered bucket, not a plain one. The constraint this pass imposes is isotropic; the
	// order it imposed it in was not. Every cell went in in flat-index order and the queue pops last-in
	// first-out within a bucket, so on ground flat to within a centimetre - which is most of a hillside -
	// cells popped bottom-right to top-left, and whichever popped first decided which of its neighbours got
	// cut. That is where the grid-aligned pyramid faces came from, and it is one hash away from not being
	// there. See bucketpq.go.
	g.pq.reset()
	for i := 0; i < n; i++ {
		x, y := i%g.W, i/g.W
		g.pq.pushJittered(h[i], int32(i), (hashXY(g.seed, g.worldX(x), g.worldY(y), jitterReposeOrder)-0.5)*2*reposeOrderBuckets)
	}

	card := talus * g.CellM
	diag := talus * g.CellM * math.Sqrt2
	var removed float64

	for g.pq.len() > 0 {
		c := g.pq.pop()
		if c < 0 {
			break
		}
		if g.closed[c] {
			continue
		}
		g.closed[c] = true
		cx, cy := int(c)%g.W, int(c)/g.W
		for k := 0; k < 8; k++ {
			nx, ny := cx+dx8[k], cy+dy8[k]
			if nx < 0 || ny < 0 || nx >= g.W || ny >= g.H {
				continue
			}
			ni := int32(ny*g.W + nx)
			if g.closed[ni] || g.fixed[ni] {
				continue
			}
			allow := card
			if dx8[k] != 0 && dy8[k] != 0 {
				allow = diag
			}
			// The same tie-break ComputeReceivers uses and for the same reason: a fixed allowance resolves
			// every near-tie the same way and prints its preferred axis. A tenth of a percent, keyed on the
			// cell being cut, so what a cell is allowed in a direction does not depend on which neighbour
			// reached it first.
			allow *= 1 + 1e-3*(float64(hashXY(g.seed, g.worldX(nx), g.worldY(ny), int32(k)+jitterReposeAllow))-0.5)
			limit := h[c] + float32(allow)
			if h[ni] > limit {
				removed += float64(h[ni] - limit)
				h[ni] = limit
				g.pq.pushJittered(limit, ni, (hashXY(g.seed, g.worldX(nx), g.worldY(ny), jitterReposeOrder)-0.5)*2*reposeOrderBuckets)
			}
		}
	}
	return removed / float64(n)
}