143 lines
4.7 KiB
Go
143 lines
4.7 KiB
Go
// 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) }
|