package check

import (
	"math"
	"testing"

	"salty/terrain/internal/fluvial"
	"salty/terrain/internal/uplift"
	"salty/terrain/internal/world"
)

// The claim the whole fault feature rests on, end to end: a difference in uplift rate across a line survives
// the solve as an escarpment, on the side the fault raises.
//
// It is here rather than in internal/uplift because everything up there tests the *rate* field - that it is
// asymmetric, that two frames agree about it, that it tapers at the tips - and none of that says the solve
// leaves anything behind. A fault is applied as a rate precisely so that erosion cannot remove it, and
// "erosion cannot remove it" is a statement about a thousand steps of stream power, not about a weight
// function. Measured on the real planet it comes out at 2.7 to 50 m of scarp for throws of 139 to 399 m, all
// five facing the right way; this is that in miniature and fast enough to run every time.
func TestAFaultLeavesAScarpAfterTheSolve(t *testing.T) {
	const w, h = 400, 400
	const cellM = 8.0
	const steps = 400
	const dtYr = 1500.0
	const runYears = steps * dtYr

	p := world.Planet{CellM: cellM, W: w, H: h, PadY: 0, NoisePeriodM: float64(w) * cellM}
	if err := p.Validate(); err != nil {
		t.Fatal(err)
	}
	f := world.Whole(p)

	// One straight east-west trace across the middle of the grid. Straight on purpose: the question is what
	// the solve does to the step, and a curve would only make the measurement harder to read.
	midM := float64(h) * cellM / 2
	pts := make([][2]float64, 17)
	for i := range pts {
		pts[i] = [2]float64{float64(i) * float64(w) * cellM / 16, midM}
	}
	trace := uplift.FaultTrace{PointsM: pts, ThrowM: 300, LengthM: float64(w) * cellM}
	delta := uplift.FaultDelta(f, []uplift.FaultTrace{trace}, runYears)
	if delta == nil {
		t.Fatal("the trace reached nothing")
	}

	// A quiet landscape to put it in: the sea along the left edge as base level, and a low uniform rate
	// everywhere else so that anything standing up is the fault's doing and not the background's.
	base := make([]bool, w*h)
	rate := make([]float32, w*h)
	height := make([]float32, w*h)
	const backgroundMYr = 4.5e-5 // 0.045 mm/yr, the shipped highland foreland
	for y := 0; y < h; y++ {
		for x := 0; x < w; x++ {
			i := y*w + x
			if x < 12 {
				base[i] = true
				continue
			}
			r := backgroundMYr + float64(delta[i])
			if r < 0 {
				r = 0
			}
			rate[i] = float32(r)
			height[i] = float32(20 + 4*math.Sin(float64(x)/23)*math.Cos(float64(y)/31))
		}
	}

	g := fluvial.NewGrid(w, h, cellM, base)
	g.SetElevationRange(-2000, 4000)
	g.Run(height, rate, nil, fluvial.Params{
		K: 5e-5, M: 0.5, N: 1, DtYr: dtYr, Steps: steps, Diffusion: 0.02, FillEvery: 1,
		TalusSlope: math.Tan(35 * math.Pi / 180), ThermalEvery: 4, ThermalPasses: 24,
		CriticalSlope: math.Tan(35 * math.Pi / 180), SlopeCap: 0.9, MaxHillslopeSub: 24,
	}, nil)

	// The trace runs east-west, so the two sides are north and south of it. nearestOnTrace signs a point by
	// the cross product, which for a west-to-east trace puts the *north* side at d > 0 - the steep, upthrown
	// side of a fault that is not reversed.
	const offCells = 75 // 600 m either side, the same offset the planet-scale measurement used
	midCell := h / 2
	mean := func(row int) float64 {
		sum, n := 0.0, 0
		for x := 40; x < w-40; x++ {
			sum += float64(height[row*w+x])
			n++
		}
		return sum / float64(n)
	}
	up := mean(midCell - offCells)
	down := mean(midCell + offCells)

	if up <= down {
		t.Fatalf("no scarp: the upthrown side averages %.1f m and the downthrown side %.1f m", up, down)
	}
	// Big enough to be terrain rather than noise, and well under the throw, because erosion takes most of a
	// fault's displacement away - which is the whole reason a fault has to be applied as a rate and not as a
	// shape. The planet-scale measurement puts the survivor at a few per cent to a fifth of the throw.
	if step := up - down; step < 5 {
		t.Errorf("the scarp is only %.1f m across a 300 m throw; that is not an escarpment", step)
	} else if step > trace.ThrowM {
		t.Errorf("the scarp is %.1f m against a %.0f m throw; nothing should exceed its own displacement",
			step, trace.ThrowM)
	}

	// And it is *at the fault*, not a general tilt of the map: the step across the trace has to be far
	// sharper than the same distance measured entirely on one side of it.
	across := up - down
	within := math.Abs(mean(midCell-offCells) - mean(midCell-2*offCells))
	if across <= within {
		t.Errorf("the step across the trace is %.1f m and a step of the same span on one side of it is "+
			"%.1f m; that is a tilted map, not a fault", across, within)
	}
}