package uplift

import (
	"math"

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

// Faults belonging to a plate boundary rather than to a painted class.
//
// The specification for this file is a map of the Alpide belt - Spain through the Maghreb, Italy, Greece,
// Turkey, Iran, Afghanistan, the Pamir, the Himalaya and into Burma - with every mapped fault trace on it.
// Six things are true of that picture and not one of them is true of a fault set scattered inside a painted
// colour:
//
//  1. **The traces are in swarms along a belt, and everywhere else is blank.** The Sahara has none. Arabia's
//     interior has none. Peninsular India, Kazakhstan, Ukraine: none. A craton is not lightly faulted, it is
//     unfaulted, and the belt next to it is saturated. Density is therefore a function of distance to a
//     boundary and of nothing else - not of which colour the ground was painted.
//  2. **The belt is wide, and how wide varies enormously.** Through Italy and Greece it is a hundred
//     kilometres; across Iran and Tibet it is well over a thousand, a fan of parallel traces from the Zagros
//     to the Alborz. So the zone is not a fixed halo: it scales with what the margin is doing.
//  3. **Faults are near the line, not on it.** Almost none of those traces *is* the plate boundary. They sit
//     tens to hundreds of kilometres either side of it, thickest near it and thinning outwards - deformation
//     is distributed across a zone, and the boundary is only where it is centred.
//  4. **Within a swarm they are sub-parallel**, to each other and to the belt, and they follow it round its
//     bends: the Turkish arc, the Zagros arc, the Himalayan arc, the fan at the Burma syntaxis. The strike
//     comes from the local tangent of the boundary, which is why the whole set curves where the margin does.
//  5. **There is a second, conjugate direction** in the wide interiors - Tibet and Mongolia show two sets
//     crossing at a high angle. One direction alone reads as corduroy, which is the defect Terrain-Next 4.A3
//     records against the procedural path at a different scale.
//  6. **They splay and anastomose** rather than running as isolated segments, which the en-echelon stepping
//     in painted_faults.go already produces and which is kept here unchanged.
//
// Everything about a *trace* - the walked heading, the taper over the last sixth, the en-echelon step past
// twelve kilometres, the escarpment profile, the repose ceiling - is shared with the class-based set through
// traceSet. What is new here is only where a fault is put and which way it points, which is exactly the part
// that was wrong.

const srcBeltFaults = 29

// beltWidth is how wide each kind of margin's deformation zone is, as a multiple of the configured width.
//
// These are ratios between kinds of boundary rather than tuning, which is why they are constants and not
// manifest keys. A continental collision has nowhere to put the convergence except into the crust on both
// sides, so it deforms a belt a thousand kilometres across; a subduction margin puts most of it down the slab
// and deforms an arc and a forearc; a transform is a narrow braid however long it runs, because the motion is
// taken up by sliding rather than by shortening; a rift deforms its two shoulders; and a mid-ocean ridge is
// the narrowest of all, an axis a few tens of kilometres wide.
var beltWidth = map[plates.Kind]float64{
	plates.Collision:  1.00,
	plates.Subduction: 0.55,
	plates.Rift:       0.35,
	plates.Transform:  0.30,
	plates.Ridge:      0.15,
}

// beltFalloff shapes how the traces thin out away from the line.
//
// An offset drawn as zone*u would spread them evenly across the whole zone, which is not what the map shows:
// the swarm is dense at the margin and trails off. Raising a uniform draw to this power biases it towards
// zero, so the density falls smoothly outwards and the zone edge is a fading-out rather than a line where
// faults stop.
const beltFalloff = 1.8

// beltThrowFloor and beltThrowCeil bound how far the closing rate is allowed to scale a throw. A margin that
// has almost stopped still has inherited structure in it, and one going twice as fast as the reference does
// not build scarps four times the size, because the repose ceiling is waiting either way.
const (
	beltThrowFloor = 0.35
	beltThrowCeil  = 2.0
)

// beltLandProbes is how many positions are sampled across the zone to find out how much of it is land.
//
// It is measured rather than assumed because the density has to keep meaning what it says. A margin running
// down the middle of an ocean and one running along a continent have the same length and the same zone area,
// and if the count came from the zone area alone the first would ask for as many traces as the second and
// then fail to place them - so the density would quietly mean something different on every boundary.
const beltLandProbes = 512

// beltPlaceTries is how many times a trace is redrawn when it lands in the sea before giving up on it.
const beltPlaceTries = 12

// beltLandShare is how much of a trace has to be on land for it to be kept. Half rather than all, because a
// fault that runs out to a coast and stops is right and a fault forbidden from reaching one is not: the
// result of demanding every probe be land is a set that avoids the shore, which is the opposite mistake.
const beltLandShare = 0.5

// BuildBeltFaults places a fault set in the deformation zones around a planet's plate boundaries.
//
// land reports whether a world position is painted land, the same callback plates.Build takes. Offshore
// faults are real - the reference map has them all over the Mediterranean and the Arabian Sea - but the solve
// fixes every ocean cell at sea level, so a trace out there changes nothing and only clutters the diagnostic.
// They are therefore kept on land, and the density is measured against the part of each zone that *is* land
// so that the number an author sets keeps meaning what it says.
func BuildBeltFaults(p world.Planet, seed int64, cfg plates.Belt, bs []plates.Boundary,
	land func(xM, yM float64) bool) []FaultTrace {

	if !cfg.Wanted() || len(bs) == 0 {
		return nil
	}
	cfg = cfg.WithDefaults()

	s := noise.NewSource(seed, srcBeltFaults)
	// The bend lattice a trace's walk turns on. Its wavelength is tied to the zone rather than to the
	// planet's fault grain: a trace inside a belt should curve on the belt's own scale.
	bendCells := int(p.NoisePeriodM/(cfg.ZoneKm*1000) + 0.5)
	if bendCells < 1 {
		bendCells = 1
	}
	pl := &beltPlacer{
		p: p, s: s, bend: noise.NewLattice(bendCells*4, s), bendCells: bendCells, cfg: cfg,
		spread: cfg.Spread() * math.Pi / 180,
		conj:   cfg.ConjugateDeg * math.Pi / 180,
		land:   land,
	}
	refM := cfg.ReferenceCmYr / 100 // cm/yr to m/yr

	var out []FaultTrace
	for bi := range bs {
		b := &bs[bi]
		if len(b.V) < 2 {
			continue
		}
		seg := beltSegments(b, cfg, refM)
		if seg.zoneKm2 <= 0 {
			continue
		}

		// How much of this belt's zone a fault can actually be placed in, measured by proposing faults
		// exactly the way the placement loop below does and counting how many survive. The same draw and the
		// same test, so the acceptance rate the loop will see is the one the count is scaled by and the
		// density keeps meaning what it says.
		kept := 0
		for i := 0; i < beltLandProbes; i++ {
			if len(pl.propose(seg)) > 0 {
				kept++
			}
		}
		usable := float64(kept) / beltLandProbes
		if usable <= 0 {
			continue
		}

		want := cfg.Per1000Km2 * seg.zoneKm2 * usable / 1000
		n := int(want)
		// Stochastic rounding, so a short margin too small for one whole fault still gets one sometimes and
		// the density means what it says averaged over a planet rather than being floored to zero.
		if s.Float() < want-float64(n) {
			n++
		}

		for i := 0; i < n; i++ {
			for try := 0; try < beltPlaceTries; try++ {
				if set := pl.propose(seg); len(set) > 0 {
					out = append(out, set...)
					break
				}
			}
		}
	}
	return out
}

// beltPlacer is everything one fault's placement needs, carried together because proposing a fault and
// measuring how often a proposal succeeds have to be the same code. See propose.
type beltPlacer struct {
	p         world.Planet
	s         *noise.Source
	bend      *noise.Lattice
	bendCells int
	cfg       plates.Belt
	spread    float64
	conj      float64
	land      func(xM, yM float64) bool
}

// propose draws one fault in the zone, builds it, and returns the traces it is made of with any that ended
// up offshore removed. An empty result is a proposal the caller should redraw.
//
// The land test is applied to the **walked geometry**, not to the straight line the fault was proposed along,
// and that distinction is the whole reason this function exists. Between a proposal and a trace sit two
// things that move it: a fault over twelve kilometres is broken into en-echelon segments staggered across
// strike, and every segment is then walked with a perturbed heading. Testing the proposal let a trace be
// accepted on a headland and then stepped and walked out into open water, which is what the first run of
// this pass drew across two straits.
//
// Segments are filtered one at a time rather than the set being kept or dropped whole, because part of a
// fault continuing offshore while the rest of it is on land is the ordinary case at any coast.
func (pl *beltPlacer) propose(seg beltSeg) []FaultTrace {
	xM, yM, tangent, halfM, side := seg.sample(pl.s)

	a := tangent + (pl.s.Float()*2-1)*pl.spread
	if pl.cfg.Conjugate() > 0 && pl.s.Float() < pl.cfg.Conjugate() {
		// The second set, crossing the first. Which way it leans is drawn per fault, because a conjugate
		// pair is two directions and picking one of them globally would be the corduroy this exists to avoid.
		if pl.s.Float() < 0.5 {
			a += pl.conj
		} else {
			a -= pl.conj
		}
	}

	lengthM := (pl.cfg.LengthKm[0] + (pl.cfg.LengthKm[1]-pl.cfg.LengthKm[0])*pl.s.Float()) * 1000
	// A wide belt carries long faults. Scaled against the configured width so that the length range an
	// author sets is the one they get on a reference-rate collision.
	lengthM *= clampF(halfM/(pl.cfg.ZoneKm*1000), 0.4, 2.2)

	throw := pl.cfg.ThrowM[0] + (pl.cfg.ThrowM[1]-pl.cfg.ThrowM[0])*pl.s.Float()

	// Vergence, and this is the point of carrying the side at all. A thrust belt is doubly vergent: the
	// faults on each flank face outwards, away from the boundary and towards the foreland they are riding
	// over. So which block goes up is decided by which side of the line the fault sits on, rather than by
	// the coin flip a class fault set has to use for want of anything better.
	set := traceSet(pl.p, pl.bend, pl.bendCells, pl.s,
		xM, yM, a, lengthM, throw*seg.throwScale, side < 0, -1)

	if pl.land == nil {
		return set
	}
	out := set[:0]
	for _, f := range set {
		if traceLandShare(f, pl.land) >= beltLandShare {
			out = append(out, f)
		}
	}
	return out
}

// traceLandShare is how much of a built trace stands on painted land.
//
// Offshore faults are real - the reference map has them throughout the Mediterranean and the Arabian Sea -
// but the solve fixes every ocean cell at sea level, so a trace out there changes no height and does nothing
// but clutter the diagnostic the fault set is read from.
func traceLandShare(f FaultTrace, land func(xM, yM float64) bool) float64 {
	if len(f.PointsM) == 0 {
		return 0
	}
	on := 0
	for _, pt := range f.PointsM {
		if land(pt[0], pt[1]) {
			on++
		}
	}
	return float64(on) / float64(len(f.PointsM))
}

// beltSeg is one boundary prepared for sampling: cumulative length along it, the zone half-width at each
// vertex, and the zone's total area.
type beltSeg struct {
	b *plates.Boundary

	// cum[i] is the length along the polyline up to vertex i, so a uniform draw over cum[last] picks a point
	// uniformly along the *line* rather than uniformly among its vertices - which would over-sample wherever
	// the chain happened to be dense.
	cum  []float64
	half []float64 // zone half-width in metres at each vertex

	zoneKm2    float64
	throwScale float64
}

// beltSegments measures a boundary: how wide its zone is at every point, how much ground that is, and how
// much the closing rate should scale the throws in it.
func beltSegments(b *plates.Boundary, cfg plates.Belt, refM float64) beltSeg {
	seg := beltSeg{b: b, cum: make([]float64, len(b.V)), half: make([]float64, len(b.V))}

	base := cfg.ZoneKm * 1000
	closingTotal := 0.0
	for i, v := range b.V {
		w := beltWidth[v.Kind]
		rate := math.Abs(v.ClosingMYr)
		if v.Kind == plates.Transform {
			// A transform closes at nothing by definition, so its zone has to be scaled by how fast it is
			// *sliding* instead. Without this every transform margin in the world would have a zone of zero
			// and the San Andreas would be unfaulted.
			rate = math.Abs(v.SlipMYr)
		}
		// The square root, not the rate itself: doubling the convergence does not double the width of the
		// belt it deforms, and a linear scale makes the fastest margin swallow a continent.
		scale := math.Sqrt(clampF(rate/refM, 0.04, 6))
		seg.half[i] = base * w * scale
		closingTotal += rate
	}
	for i := 1; i < len(b.V); i++ {
		d := math.Hypot(b.V[i].XM-b.V[i-1].XM, b.V[i].YM-b.V[i-1].YM)
		seg.cum[i] = seg.cum[i-1] + d
		// The zone either side of this segment, as a trapezium on each flank.
		seg.zoneKm2 += d * (seg.half[i-1] + seg.half[i]) / 1e6
	}
	if n := len(b.V); n > 0 {
		mean := closingTotal / float64(n)
		seg.throwScale = clampF(mean/refM, beltThrowFloor, beltThrowCeil)
	}
	return seg
}

// sample draws one position in the zone: a point along the line, then an offset across it.
//
// It returns the world position, the belt's local strike there, the local zone half-width, and which side of
// the line the point fell on. The side is what vergence is read from, and the half-width is what a trace's
// length is scaled by.
func (seg beltSeg) sample(s *noise.Source) (xM, yM, strike, halfM, side float64) {
	total := seg.cum[len(seg.cum)-1]
	if total <= 0 {
		return seg.b.V[0].XM, seg.b.V[0].YM, 0, seg.half[0], 1
	}
	at := s.Float() * total
	// Walk to the segment holding it. Linear rather than a binary search on purpose: a boundary is a few
	// hundred vertices and this runs a few thousand times, so the search is not where the time goes and a
	// loop with no off-by-one in it is worth more here than the log.
	i := 1
	for i < len(seg.cum)-1 && seg.cum[i] < at {
		i++
	}
	t := 0.0
	if d := seg.cum[i] - seg.cum[i-1]; d > 0 {
		t = (at - seg.cum[i-1]) / d
	}

	a, b := seg.b.V[i-1], seg.b.V[i]
	lx := a.XM + (b.XM-a.XM)*t
	ly := a.YM + (b.YM-a.YM)*t
	nx := a.NX + (b.NX-a.NX)*t
	ny := a.NY + (b.NY-a.NY)*t
	if d := math.Hypot(nx, ny); d > 0 {
		nx, ny = nx/d, ny/d
	}
	halfM = seg.half[i-1] + (seg.half[i]-seg.half[i-1])*t

	// Across the line. The offset is biased towards zero so the swarm is dense at the margin and trails off,
	// and the side is drawn separately so both flanks are populated.
	side = 1
	if s.Float() < 0.5 {
		side = -1
	}
	off := halfM * math.Pow(s.Float(), beltFalloff) * side

	// The strike is the boundary's own tangent, which is the perpendicular of the normal. This single line is
	// the whole difference between a swarm that follows the Zagros round its arc and a set of traces pointing
	// wherever a noise lattice happened to say.
	strike = math.Atan2(-nx, ny)

	return lx + nx*off, ly + ny*off, strike, halfM, side
}

func clampF(v, lo, hi float64) float64 {
	if v < lo {
		return lo
	}
	if v > hi {
		return hi
	}
	return v
}