// Package thermal is mass-conserving thermal weathering at an angle of repose: where a cell stands above a
// neighbour by more than the angle allows, material slides down.
//
// It is not decoration, and this is worth stating because it looks like decoration. Stream power puts *no
// upper bound* on hillslope angle: the steady state S = (U/K)^(1/n) * A^(-m/n) is evaluated per cell, and a
// cell next to a divide drains one cell, so A is tiny and S is enormous. Raising the intraplate uplift rate
// to give the plains some relief therefore made every hillslope on the map proportionally steeper, and the
// median slope went from 10 to 39 degrees. What limits a real hillslope is not fluvial incision, it is
// landsliding, and this is landsliding.
//
// Ported from heightmap_erosion.thermal, including the lesson attached to it: shed half the *largest* excess
// rather than half the mean, because the mean converges far slower. Mass is conserved, so a cliff keeps its
// face and scree gathers at its foot instead of evaporating.
package thermal

import (
	"math"

	"salty/terrain/internal/field"
)

var (
	dx8 = [8]int{-1, 0, 1, -1, 1, -1, 0, 1}
	dy8 = [8]int{-1, -1, -1, 0, 0, 1, 1, 1}
)

// Apply runs `passes` iterations over h in place. talus is the maximum slope as rise over run (tan of the
// angle of repose). fixed marks cells that may not move, which is base level.
//
// Each pass reads h and writes a delta, then applies it: no cell is read after any cell has been written, so
// the result does not depend on the order the rows were processed and the parallel version is identical to
// the serial one.
func Apply(h []float32, w, hgt int, cellM, talus float64, passes int, fixed []bool, scratch []float32) {
	if passes <= 0 || talus <= 0 {
		return
	}
	// The scratch is the caller's so a pass inside a loop does not allocate a grid every time; it is optional
	// because every other caller of this package passes one and the one that did not spent its first run in a
	// panic (slice bounds out of range) rather than in a weather simulation.
	if cap(scratch) < len(h) {
		scratch = make([]float32, len(h))
	}
	delta := scratch[:len(h)]
	card := cellM
	diag := cellM * math.Sqrt2

	for p := 0; p < passes; p++ {
		for i := range delta {
			delta[i] = 0
		}
		// Deltas are accumulated per source cell into a private view; each row range owns its own writes for
		// the cell it is standing on, and the neighbour credits are gathered rather than scattered, which is
		// what keeps this free of races without a lock.
		field.Rows(hgt, func(y0, y1 int) {
			for y := y0; y < y1; y++ {
				for x := 0; x < w; x++ {
					i := y*w + x
					// Gather: how much this cell RECEIVES from higher neighbours, and how much it sheds.
					var give, take float64
					if fixed == nil || !fixed[i] {
						give = shed(h, w, hgt, x, y, card, diag, talus)
					}
					take = 0
					for k := 0; k < 8; k++ {
						nx, ny := x+dx8[k], y+dy8[k]
						if nx < 0 || ny < 0 || nx >= w || ny >= hgt {
							continue
						}
						ni := ny*w + nx
						if fixed != nil && fixed[ni] {
							continue
						}
						// The share this cell gets of that neighbour's shed material.
						if share := shedShare(h, w, hgt, nx, ny, i, card, diag, talus); share > 0 {
							take += share
						}
					}
					delta[i] = float32(take - give)
				}
			}
		})
		for i := range h {
			h[i] += delta[i]
		}
	}
}

// shed is the total a cell gives away this pass: half the largest excess over the angle of repose.
func shed(h []float32, w, hgt, x, y int, card, diag, talus float64) float64 {
	i := y*w + x
	maxExcess := 0.0
	for k := 0; k < 8; k++ {
		nx, ny := x+dx8[k], y+dy8[k]
		if nx < 0 || ny < 0 || nx >= w || ny >= hgt {
			continue
		}
		d := card
		if dx8[k] != 0 && dy8[k] != 0 {
			d = diag
		}
		if e := float64(h[i]-h[ny*w+nx]) - talus*d; e > maxExcess {
			maxExcess = e
		}
	}
	return maxExcess / 2
}

// shedShare is how much of cell (x, y)'s shed material lands on cell target, proportional to target's share
// of the total excess below it.
func shedShare(h []float32, w, hgt, x, y, target int, card, diag, talus float64) float64 {
	i := y*w + x
	maxExcess, total, mine := 0.0, 0.0, 0.0
	for k := 0; k < 8; k++ {
		nx, ny := x+dx8[k], y+dy8[k]
		if nx < 0 || ny < 0 || nx >= w || ny >= hgt {
			continue
		}
		d := card
		if dx8[k] != 0 && dy8[k] != 0 {
			d = diag
		}
		ni := ny*w + nx
		e := float64(h[i]-h[ni]) - talus*d
		if e <= 0 {
			continue
		}
		if e > maxExcess {
			maxExcess = e
		}
		total += e
		if ni == target {
			mine = e
		}
	}
	if total <= 0 || mine <= 0 {
		return 0
	}
	return (maxExcess / 2) * (mine / total)
}

// TalusFromDegrees converts an angle of repose to the rise-over-run the solver wants.
func TalusFromDegrees(deg float64) float64 { return math.Tan(deg * math.Pi / 180) }