package uplift

import (
	"math"

	"salty/terrain/internal/noise"
	"salty/terrain/internal/world"
)

// Faults on a painted planet: the structure inside a range that the painting cannot draw.
//
// A painted class is one rate over every cell of a colour, and a massif block breaks that into plain and
// upland. Neither can produce the thing a real range has running through it - an escarpment, a block tilted
// against its neighbour, a valley that is straight for fifteen kilometres because it is following a break in
// the rock. That is a fault, and a fault is not a shape to paint: it is a *difference in uplift rate across a
// line*, steep on one side and gentle on the other, which erosion then carves into a scarp. Paint it as
// terrain and the solve erodes it away; apply it as a rate and it maintains itself, which is the same
// argument as "paint the uplift, never the height" one scale down.
//
// The procedural path has had this since the beginning and none of it could be carried over as it stood,
// for two reasons that are by now familiar.
//
// **Placement is in world metres, not map fractions.** `uplift.Build` draws a trace centre as two calls to
// `s.Float()`, which are fractions of the grid it happens to be filling. On a decomposed planet that is a
// different place in every region. Here the whole set is drawn once for the planet, in metres east of the
// seam and metres south of the top painted row, and a region filters it to the traces that reach into its
// own frame - so a fault that crosses a region boundary is one fault, and two different decompositions of the
// same planet produce the same escarpment.
//
// **And a trace is placed in the ground the author painted for it.** A class carries `faults`, so an author
// says *this range is faulted and that plain is not*, which is both the control they want and what is
// actually true of the world: faults belong to orogens. The influence is not restricted to the class,
// because a range-front fault runs along the edge of a range by definition and its scarp faces the lowland.
//
// The four defects recorded against the procedural version in Terrain-Next 4.A2 are fixed here rather than
// carried, because writing a new implementation with a known fault list is cheaper than porting one and
// fixing it afterwards:
//
//   - The trace is a walk with a *perturbed heading* rather than one 8-point parabola, so the distance field
//     around it has no polygonal contours.
//   - The throw is **tapered to zero over the last sixth at each tip** rather than stopping dead where the
//     last segment ends, which is what made a fault cut abruptly across a summit.
//   - A long fault is broken into overlapping en-echelon segments, which is how long faults actually step.
//   - And nothing is clamped to a fraction of a global rate. The procedural version flattens its strongest
//     throws with `if r > convergent*1.6`, which turns exactly the faults that matter most into plateaus.
//     The ceiling here is the one with a physical meaning - the rate at which a divide reaches the angle of
//     repose - and the count of cells that reach it is reported rather than hidden.

const (
	srcPaintFaults = 27
	srcPaintGrain  = 28
)

// The escarpment's shape, in metres, and the reason it is that shape (D-62).
//
// The first version of this made the rate difference a *step*: the whole throw on one side of the trace and
// the whole throw negated on the other, one cell apart, the positive side decaying to nothing over six
// hundred metres and the negative side over six kilometres. Measured on Bake_013 that is a wall of two
// throws across a single 8 m cell, standing in a welt six hundred metres wide. Two things follow, and both
// are visible in a hillshade before any number is taken.
//
// **A step in the rate is a painted cliff.** "Paint the uplift, never the height" is a claim about what a
// solve can undo, and a discontinuity in the rate field is precisely what it cannot: the surface has
// nowhere to put the difference but into a scarp at the angle of repose, so the trace comes out as a
// facetted line at any throw, and turning throw_m down only lowers the same artefact.
//
// **And six hundred metres is narrower than one hillslope.** Bake_013's drainage density is 0.45 channels
// per kilometre, so a divide sits about 1.1 km from its channel. Nothing can dissect a block six hundred
// metres wide - there is no drainage area at that width for stream power to work with, and hillslope
// diffusion only smooths what is already there - so the uplift profile is *printed* onto the surface
// rather than eroded into a landform. That is why every fault in that bake reads as a smooth ruled ridge
// running through terrain dissected everywhere else: it is the one part of the map erosion never touched.
//
// So the profile is antisymmetric, continuous through the trace, and kilometres wide on both flanks.
//
//   - faultRampM is how far the rate takes to cross from the hanging wall to the footwall - about one
//     hillslope length, which makes the mountain front the sharpest thing this landscape can express
//     without it being a cliff nobody solved for.
//   - faultFootwallM and faultHangingM are how far each flank reaches. Several hillslope lengths, so a
//     drainage network fits on the block and cuts it into spurs and valleys, which is what a range front
//     is and what an extruded cross-section will never be.
//
// Neither width scales with the trace's length, and that is deliberate twice over. Physically, the width
// of flexural footwall uplift is set by how the crust bends rather than by the fault in it, so a short
// fault on the same lithosphere makes a *lower* range, not a narrower one - which is what a throw does
// here already. And practically, the floor is the one that matters: a flank has to be wide enough for a
// drainage network whatever the trace is, and scaling it down for the planet's shortest traces - 1.7 km
// on the shipped template against a 6.4 km median - would put the printing artefact straight back on
// exactly those. `length_km` says how far a fault runs along strike; it does not say how wide a belt it
// deforms.
//
// The anomaly is therefore zero *on the trace itself*, which is also the honest reading: a rate difference
// across a line says one side rises relative to the other, and at the line the two average to the regional
// rate. The old profile asserted +throw and -throw at the same point.
const (
	faultRampM     = 900.0
	faultFootwallM = 6000.0
	faultHangingM  = 4000.0
)

// faultReachM is where the influence is cut off.
//
// The flank envelope reaches zero *with zero gradient* at its own width, so unlike the old
// exponential-minus-a-floor there is nothing to subtract and no step at the box edge to hide: the cut-off
// is the support of the function rather than a truncation of it.
const faultReachM = faultFootwallM

// faultShape is the unnormalised profile at a signed distance from the trace, positive on the upthrown
// side: an odd saturating ramp across the line, times a flank envelope.
//
// The ramp is d/sqrt(R*R+d*d) rather than tanh and the envelope is (1-u*u)^2 rather than an exponential,
// because this runs at every cell of every fault's box - a few hundred million times on a planet - and
// neither transcendental buys anything over the algebraic pair.
func faultShape(d float64) float64 {
	w := faultFootwallM
	if d < 0 {
		w = faultHangingM
	}
	u := math.Abs(d) / w
	if u >= 1 {
		return 0
	}
	e := 1 - u*u
	return d / math.Sqrt(faultRampM*faultRampM+d*d) * e * e
}

// faultNorm scales the profile so the whole step across a fault - the footwall crest less the hanging wall
// trough - is exactly the throw the author asked for, which is what the word means: the vertical
// displacement across the fault. The old profile put a full throw on each side and so built two.
//
// Measured over the profile rather than written down, so that changing a width above cannot silently
// change what throw_m means.
var faultNorm = func() float64 {
	up, down := 0.0, 0.0
	for d := 1.0; d < faultReachM; d++ {
		if v := faultShape(d); v > up {
			up = v
		}
		if v := -faultShape(-d); v > down {
			down = v
		}
	}
	return up + down
}()

// faultWeight is the escarpment profile, scaled so the crest-to-trough step across it is one throw.
func faultWeight(d float64) float64 { return faultShape(d) / faultNorm }

// faultStackBonus is how much more than its strongest single fault a whole stack of them may build,
// and it is the answer to the defect D-62 caused (D-63).
//
// Faults are rasterised with `+=`, which was harmless while a fault reached six hundred metres: two of
// them almost never met. At six kilometres they meet constantly, and a fault set is *sub-parallel by
// construction* - traces within one cell of the orientation grain share a strike, and a belt fault takes
// its strike from the plate margin - so where they meet they are all pushing the same way. Measured on
// Bake_018's region 11, twenty-two kilometres across with thirteen traces at strikes spanning fourteen
// degrees: **75 % of the faulted ground had two or more faults on it**, the sum was a median 1.77 times
// the largest single contribution there and up to 4.46 times, and the landmass came out at 221 m against
// 75 m for the same ground before D-62. Thirteen per cent of it asked for more uplift than the repose
// ceiling allows on its own, so the hard clamp downstream fired on 160 289 cells - 4.1 % of the region,
// against 0.15 % over the whole planet before - and a hard clamp makes plateaus.
//
// The knee is the **largest single contribution at that cell**, not the largest throw in the set. A
// planet-wide throw would not bite: the biggest throw on this template is 744 m while the biggest single
// contribution anywhere in region 11 is 209 m, because a fault's own taper and falloff have already
// reduced it by the time it reaches anywhere. Keyed per cell, one fault passes through untouched and only
// the stacking is bent.
const faultStackBonus = 0.6

// softStack combines the summed anomaly at a cell with the largest single contribution there.
//
// Below the knee it is the identity, so a cell reached by one fault gets exactly what that fault asked
// for and D-62's calibration - the step across a fault is its throw - is unchanged. Above it the excess
// is bent through a tanh onto an asymptote of (1+faultStackBonus) times the knee, so a belt still stands
// higher than an unfaulted one, which is the point of a belt, but five parallel faults cannot deliver
// five throws. Odd in `sum`, so a stack of hanging walls is bounded on the same terms.
//
// It is continuous: `peak` is a max of continuous functions and the join at the knee has gradient 1 on
// both sides. And it is frame-independent, which is what keeps the decomposition honest - `sum` and
// `peak` at a cell depend only on the faults within reach of it, and a fault too far away to be in a
// frame's box contributes nothing to either.
func softStack(sum, peak float64) float64 {
	if peak <= 0 {
		return 0
	}
	a := math.Abs(sum)
	if a <= peak {
		return sum
	}
	head := faultStackBonus * peak
	v := peak + head*math.Tanh((a-peak)/head)
	if sum < 0 {
		return -v
	}
	return v
}

// enEchelonM is the length past which a fault is drawn as overlapping segments instead of one line, and
// enEchelonSpan is how long each segment is as a fraction of the parent.
const (
	enEchelonM    = 12000.0
	enEchelonSpan = 0.55
)

// FaultSpec is what one painted class asks for. A class with no spec has no faults.
type FaultSpec struct {
	// Per1000Km2 is how many traces to place in every thousand square kilometres of this class. It is a
	// density rather than a count because a class covers whatever an author painted, and a count would mean
	// something different on every template.
	Per1000Km2 float64

	// ThrowM is the whole step across the fault over the run, low to high: the footwall crest less the
	// hanging wall trough, which is what the word throw means. It becomes a rate - the solve integrates it
	// for `steps * dt_yr` years - so what an author picks is how much higher the upthrown side would stand
	// than the downthrown one if erosion never touched either. Before D-62 the profile put a whole throw on
	// each flank and so built two of them, and a legend written against that asks for half what it did.
	ThrowM [2]float64

	// LengthKm is how long a trace is, low to high.
	LengthKm [2]float64
}

// Wanted reports whether this spec asks for anything.
func (s FaultSpec) Wanted() bool {
	return s.Per1000Km2 > 0 && s.LengthKm[1] > 0 && s.ThrowM[1] > 0
}

// FaultTrace is one fault, in world metres.
//
// X is **unwrapped**: a trace that crosses the seam has X running past the circumference or below zero rather
// than jumping, so that every segment of it is a straight line between neighbouring points and no consumer
// has to special-case the meridian. Whoever draws or tests it wraps by the circumference.
type FaultTrace struct {
	PointsM [][2]float64 `json:"points_m"`
	ThrowM  float64      `json:"throw_m"`
	LengthM float64      `json:"length_m"`
	Class   int          `json:"class"`

	// Reverse flips which side goes up. Half of them do, drawn from the same stream, because a fault set in
	// which every block tilts the same way reads as corduroy.
	Reverse bool `json:"reverse"`
}

// BuildFaults draws the planet's whole fault set, once, deterministically from the seed.
//
// candidates[c] is a strided sample of the planet cells belonging to class c: the list a trace centre is
// drawn from, so a fault lands in the ground its class was painted on. Sampled rather than enumerated because
// the full list for a class covering a seventh of a 76-million-cell planet is ten million entries, and the
// only thing asked of it is a uniform draw. areaCells is the class's *exact* cell count, which the projection
// already counted, so the density is not estimated from the sample.
//
// grainKm is the wavelength of the orientation field. Faults within one of its cells come out sub-parallel
// and the set swings gradually across the world, which is what a fault set looks like and what a single
// global strike angle - the procedural path's `grainAngle` - does not: Terrain-Next 4.A3 records that one
// running as straight corduroy across a whole map.
func BuildFaults(p world.Planet, seed int64, grainKm float64, specs []FaultSpec,
	candidates [][]int32, areaCells []int) []FaultTrace {

	if grainKm <= 0 {
		return nil
	}
	any := false
	for _, s := range specs {
		any = any || s.Wanted()
	}
	if !any {
		return nil
	}

	grainCells := int(p.NoisePeriodM/(grainKm*1000) + 0.5)
	if grainCells < 1 {
		grainCells = 1
	}
	gs := noise.NewSource(seed, srcPaintGrain)
	// Two lattices read as a vector rather than one read as an angle. A value lattice runs 0..1 and an angle
	// taken straight from it jumps by a whole turn wherever it crosses its own wrap, which would put a hard
	// seam through the fault set along a contour nobody can see. atan2 of two fields is continuous.
	gx := noise.NewLattice(grainCells, gs)
	gy := noise.NewLattice(grainCells, gs)
	// The bend lattice is finer, so a trace curves within the province its strike came from.
	bend := noise.NewLattice(grainCells*4, gs)

	s := noise.NewSource(seed, srcPaintFaults)
	cellArea := p.CellM * p.CellM

	var out []FaultTrace
	for c := range specs {
		spec := specs[c]
		if !spec.Wanted() || c >= len(candidates) || len(candidates[c]) == 0 || c >= len(areaCells) {
			continue
		}
		areaKm2 := float64(areaCells[c]) * cellArea / 1e6
		want := spec.Per1000Km2 * areaKm2 / 1000
		// Stochastic rounding, so a class too small for one whole fault still gets one sometimes and the
		// density means what it says when averaged over a world rather than being floored to zero.
		n := int(want)
		if s.Float() < want-float64(n) {
			n++
		}
		for i := 0; i < n; i++ {
			cell := candidates[c][s.IntN(len(candidates[c]))]
			cx := float64(int(cell)%p.W) * p.CellM
			cy := p.YM(int(cell) / p.W)
			lengthM := spec.LengthKm[0] + (spec.LengthKm[1]-spec.LengthKm[0])*s.Float()
			lengthM *= 1000
			throw := spec.ThrowM[0] + (spec.ThrowM[1]-spec.ThrowM[0])*s.Float()
			reverse := s.Float() < 0.5
			a := strikeAt(p, gx, gy, grainCells, cx, cy)
			out = append(out, traceSet(p, bend, grainCells, s, cx, cy, a, lengthM, throw, reverse, c)...)
		}
	}
	return out
}

// traceSet turns one drawn fault into the one or more traces it is actually made of.
//
// The strike is an argument rather than something this function looks up. Where a fault points is the whole
// difference between the two placements that call it: a class fault takes its angle from a noise grain, and a
// belt fault takes it from the plate boundary it belongs to. Everything below - the en-echelon step, the
// walk, the taper - is the same fault either way.
func traceSet(p world.Planet, bend *noise.Lattice, grainCells int, s *noise.Source,
	cx, cy, a, lengthM, throw float64, reverse bool, class int) []FaultTrace {

	if lengthM <= enEchelonM {
		return []FaultTrace{walkTrace(p, bend, grainCells, cx, cy, a, lengthM, throw, reverse, class)}
	}
	// A long fault steps. Two or three overlapping segments, each a little over half the parent's length,
	// staggered along strike and offset across it - which is what a long fault does in the ground and is also
	// the difference between a fifteen-kilometre ruled line and something that reads as structure.
	n := 2
	if s.Float() < 0.5 {
		n = 3
	}
	segLen := lengthM * enEchelonSpan
	out := make([]FaultTrace, 0, n)
	for i := 0; i < n; i++ {
		// Centres spread along the parent, from -0.5 to +0.5 of its length.
		t := (float64(i)/float64(n-1) - 0.5) * (lengthM - segLen)
		lateral := (s.Float() - 0.5) * 0.12 * lengthM
		sx := cx + math.Cos(a)*t - math.Sin(a)*lateral
		sy := cy + math.Sin(a)*t + math.Cos(a)*lateral
		out = append(out, walkTrace(p, bend, grainCells, sx, sy, a, segLen, throw, reverse, class))
	}
	return out
}

// strikeAt is the fault grain's direction at a world position, as an angle.
func strikeAt(p world.Planet, gx, gy *noise.Lattice, cells int, xM, yM float64) float64 {
	u := xM / p.NoisePeriodM * float64(cells)
	v := yM / p.NoisePeriodM * float64(cells)
	return math.Atan2(float64(gy.Sample(u, v))-0.5, float64(gx.Sample(u, v))-0.5)
}

// walkTrace steps outward from the centre in both directions, turning a little each step.
//
// A walk rather than a formula. The procedural path draws one parabola through eight points, and a distance
// field built from eight long straight segments has visibly polygonal contours - which is the first of the
// four things Terrain-Next 4.A2 lists against it. Short steps with a heading that wanders make the same
// gentle curve with none of that.
func walkTrace(p world.Planet, bend *noise.Lattice, grainCells int,
	cx, cy, a, lengthM, throw float64, reverse bool, class int) FaultTrace {

	segs := int(lengthM / 500)
	if segs < 8 {
		segs = 8
	}
	if segs > 64 {
		segs = 64
	}
	if segs%2 == 1 {
		segs++ // even, so the centre is a point rather than the middle of a segment
	}
	step := lengthM / float64(segs)
	// The heading turns by at most this much per step. Correlated through the lattice, so consecutive steps
	// see nearly the same value and the walk integrates into a smooth arc rather than a jitter.
	const maxTurn = 0.12

	pts := make([][2]float64, segs+1)
	mid := segs / 2
	pts[mid] = [2]float64{cx, cy}
	for dir := -1; dir <= 1; dir += 2 {
		x, y, ang := cx, cy, a
		if dir < 0 {
			ang += math.Pi
		}
		for k := 1; k <= mid; k++ {
			u := x / p.NoisePeriodM * float64(grainCells*4)
			v := y / p.NoisePeriodM * float64(grainCells*4)
			ang += (float64(bend.Sample(u, v))*2 - 1) * maxTurn * float64(dir)
			x += math.Cos(ang) * step
			y += math.Sin(ang) * step
			pts[mid+dir*k] = [2]float64{x, y}
		}
	}
	return FaultTrace{PointsM: pts, ThrowM: throw, LengthM: lengthM, Class: class, Reverse: reverse}
}

// FaultDelta is the uplift-rate change the fault set contributes over one frame, in metres a year, or nil
// when none of them reaches it.
//
// Rasterised per fault into its own box rather than per cell over every fault: the set is planet-wide, so a
// per-cell loop over all of it would be the whole planet's faults tested at every cell of every region. A
// fault's box is its trace plus faultReachM on all sides, clipped to the frame, and inside it each cell tests
// only the segments whose own boxes contain it.
func FaultDelta(f world.Frame, faults []FaultTrace, runYears float64) []float32 {
	if len(faults) == 0 || runYears <= 0 {
		return nil
	}
	cellM := f.P.CellM
	circ := f.P.CircumferenceM()
	x0M, y0M := f.OriginXM(), f.OriginYM()
	x1M, y1M := x0M+float64(f.W)*cellM, y0M+float64(f.H)*cellM

	// peak is the largest single contribution at each cell, which is the knee softStack bends the sum
	// over. Carried alongside rather than derived afterwards because a second pass over every fault would
	// cost exactly what the first one did.
	var out, peak []float32
	for fi := range faults {
		ft := &faults[fi]
		if len(ft.PointsM) < 2 {
			continue
		}
		// Shift the trace by whole turns of the planet so it sits nearest this frame. X is unwrapped in the
		// stored trace, so this is the one place the seam is dealt with, once per fault instead of per cell.
		pts := shiftToFrame(ft.PointsM, (x0M+x1M)/2, circ)

		lo, hi := traceBounds(pts)
		if lo[0]-faultReachM > x1M || hi[0]+faultReachM < x0M ||
			lo[1]-faultReachM > y1M || hi[1]+faultReachM < y0M {
			continue
		}
		if out == nil {
			out = make([]float32, f.W*f.H)
			peak = make([]float32, f.W*f.H)
		}
		cxa := clampInt(int((lo[0]-faultReachM-x0M)/cellM), 0, f.W-1)
		cxb := clampInt(int((hi[0]+faultReachM-x0M)/cellM)+1, 0, f.W-1)
		cya := clampInt(int((lo[1]-faultReachM-y0M)/cellM), 0, f.H-1)
		cyb := clampInt(int((hi[1]+faultReachM-y0M)/cellM)+1, 0, f.H-1)

		sign := 1.0
		if ft.Reverse {
			sign = -1
		}
		rate := ft.ThrowM / runYears
		// One box per segment, computed here rather than inside the cell loop. The rejection below runs for
		// every segment at every cell in the fault's box - a few hundred million times on a large region -
		// and recomputing four min/max per test was most of what the pass cost.
		boxes := segmentBoxes(pts)
		for y := cya; y <= cyb; y++ {
			py := y0M + float64(y)*cellM
			row := y * f.W
			for x := cxa; x <= cxb; x++ {
				px := x0M + float64(x)*cellM
				d, along, ok := nearestOnTrace(px, py, pts, boxes)
				if !ok {
					continue
				}
				v := float32(rate * faultWeight(d*sign) * tipTaper(along))
				out[row+x] += v
				if v < 0 {
					v = -v
				}
				if v > peak[row+x] {
					peak[row+x] = v
				}
			}
		}
	}
	// Bend the stacking. Unconditional whenever anything was rasterised, including when a single trace
	// reached this frame: skipping it there would make a cell's value depend on which frame it was asked
	// about, which is the one thing FaultDelta is not allowed to do.
	for i := range out {
		if peak[i] > 0 {
			out[i] = float32(softStack(float64(out[i]), float64(peak[i])))
		}
	}
	return out
}

// tipTaper is how much of its throw a fault carries at a fraction along its length.
//
// It was a ramp over the last sixth at each end with a flat top over the middle two thirds. That dies out
// at the tips, which is what it was written for, and leaves the cross-section above *extruded* unchanged
// along two thirds of every trace - which is the other half of why a fault reads as a ruled line. An
// extrusion has no along-strike structure, so erosion has no reason to head a valley in one place rather
// than another and the ridge stays as smooth as the function that drew it.
//
// A real fault carries most of its displacement near the middle and none at either tip, in a profile
// somewhere between elliptical and a linear taper. This is that: an ellipse in q, lifted by q*(2-q) so the
// middle four fifths keeps a body rather than coming to a point. Polynomial on purpose - it is evaluated
// at every cell of every fault's box.
func tipTaper(along float64) float64 {
	if along <= 0 || along >= 1 {
		return 0
	}
	q := 4 * along * (1 - along)
	return q * (2 - q)
}

// nearestOnTrace is the signed perpendicular distance to a polyline and how far along it the nearest point
// sits, 0 at one tip and 1 at the other. ok is false beyond the ends, where a fault has no effect.
func nearestOnTrace(px, py float64, pts [][2]float64, boxes [][4]float64) (dist, along float64, ok bool) {
	best := math.Inf(1)
	sign := 1.0
	at := 0.0
	n := len(pts) - 1
	for j := 0; j < n; j++ {
		// Cheap rejection first: this loop runs for every cell in the fault's box, and on most of them every
		// segment misses.
		b := &boxes[j]
		if px < b[0] || px > b[2] || py < b[1] || py > b[3] {
			continue
		}
		ax, ay := pts[j][0], pts[j][1]
		bx, by := pts[j+1][0], pts[j+1][1]
		dx, dy := bx-ax, by-ay
		l2 := dx*dx + dy*dy
		if l2 < 1e-9 {
			continue
		}
		t := ((px-ax)*dx + (py-ay)*dy) / l2
		if t < 0 || t > 1 {
			continue // beyond this segment; a neighbouring one may still claim the point
		}
		projx, projy := ax+t*dx, ay+t*dy
		d := math.Hypot(px-projx, py-projy)
		if d < best {
			best = d
			at = (float64(j) + t) / float64(n)
			if (px-ax)*dy-(py-ay)*dx < 0 {
				sign = -1
			} else {
				sign = 1
			}
		}
	}
	if math.IsInf(best, 1) {
		return 0, 0, false
	}
	return best * sign, at, true
}

// segmentBoxes is each segment's own box, grown by the reach: the rejection test in nearestOnTrace.
func segmentBoxes(pts [][2]float64) [][4]float64 {
	out := make([][4]float64, len(pts)-1)
	for j := range out {
		ax, ay := pts[j][0], pts[j][1]
		bx, by := pts[j+1][0], pts[j+1][1]
		out[j] = [4]float64{
			math.Min(ax, bx) - faultReachM, math.Min(ay, by) - faultReachM,
			math.Max(ax, bx) + faultReachM, math.Max(ay, by) + faultReachM,
		}
	}
	return out
}

// shiftToFrame moves a trace by whole circumferences so its middle is nearest a given longitude.
func shiftToFrame(pts [][2]float64, centreXM, circ float64) [][2]float64 {
	mid := pts[len(pts)/2][0]
	k := math.Round((centreXM - mid) / circ)
	if k == 0 {
		return pts
	}
	out := make([][2]float64, len(pts))
	for i, p := range pts {
		out[i] = [2]float64{p[0] + k*circ, p[1]}
	}
	return out
}

func traceBounds(pts [][2]float64) (lo, hi [2]float64) {
	lo = [2]float64{math.Inf(1), math.Inf(1)}
	hi = [2]float64{math.Inf(-1), math.Inf(-1)}
	for _, p := range pts {
		lo[0] = math.Min(lo[0], p[0])
		lo[1] = math.Min(lo[1], p[1])
		hi[0] = math.Max(hi[0], p[0])
		hi[1] = math.Max(hi[1], p[1])
	}
	return lo, hi
}

func clampInt(v, lo, hi int) int {
	if v < lo {
		return lo
	}
	if v > hi {
		return hi
	}
	return v
}