UnrealPrototyping/Tools/Terrain/internal/plates/plates.go

// Package plates is the tectonics a painted planet does not draw: which rigid pieces the lithosphere is in,
// how they move, and therefore where they are colliding.
//
// It exists because of one row in Docs/Terrain.md's pass table. Pass 1 writes `uplift` *and* `boundaries`,
// and pass 3 - faults - reads `boundaries`. D-53 dropped passes 1 to 4 on the painted path, because a
// painted template is already a statement about where the ranges are. What went out with them was the
// boundary set, and `uplift/painted_faults.go` substituted a noise grain field for it - a field that has
// never been told where a belt is. That substitution is visible in `map_uplift.png`: cyan traces striking
// across the bright belts at angles unrelated to them, the densest set sitting in a lowland, and several
// walking out over open ocean.
//
// The correction is not a better grain field. A range and its faults are not two things, one decorating the
// other: they are both consequences of the same convergence, and the line they are consequences of is the
// plate boundary. So the boundary is what gets built first, and the uplift and the faults are both read off
// it.
//
// **The plates live on the cylinder, not on a sphere.** A real plate moves by rotating about an Euler pole
// through the centre of the planet, and the velocity that produces varies along a boundary - which is the
// reason one margin is a head-on collision at one end and a strike-slip fault at the other. That variation
// is worth having; the sphere is not. Every other pass here measures distance in flat metres on a cylinder
// of fixed circumference with an 8 m cell that never varies (D-48), so a pass that believed in a sphere
// would be the only one whose distances disagreed with the solve's, and its velocities would converge at
// poles nothing else knows are there. The compromise keeps the property and drops the geometry: a plate's
// motion is a translation plus a rotation about a pole **in the map plane**,
//
//	v(x) = T + omega x (x - pole)
//
// which is the two-dimensional analogue and varies along a boundary for the same reason.
//
// **What the painting still owns.** Whether a plate is continental is read from the land mask rather than
// drawn from the seed: a plate covering the author's continent *is* a continental plate. That is the one
// place the painting feeds the model rather than competing with it, and it is what makes an ocean-continent
// margin land where an author would expect a subduction zone.
//
// **The tectonic grid is its own, and coarse.** The partition is rasterised at a few hundred metres rather
// than at the 8 m geology cell. A plate boundary belt is tens of kilometres wide and the finest thing read
// off the line is a fault trace, so a quarter-kilometre lattice is already finer than anything downstream
// can use, and it makes the whole pass a few million operations instead of a few hundred million. What
// leaves this package is polylines in **world metres**, which is the same form `uplift.FaultTrace` already
// travels in and for the same reason: a region filters the planet's set to what reaches its own frame, so a
// boundary crossing a region edge is one boundary and two decompositions agree.
package plates

import (
	"fmt"
	"math"

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

// Pass indices for this package's seeded streams. They sit above the detail passes' 40s so that adding one
// here cannot reshuffle any existing field.
const (
	srcSites  = 50
	srcMotion = 51
	srcWarp   = 52
)

// Config is what the manifest asks for. Every zero field takes a default from withDefaults.
type Config struct {
	// Layer and Legend are the painted tectonic layer: an image where a colour is a plate, and a legend
	// saying how each one moves. Naming them is what turns the plates from something the seed invents into
	// something an author draws, and it is the intended way to use this package - see paint.go.
	//
	// With no layer the plates come from Count and the seed, which is a Voronoi partition that knows nothing
	// about where the continents are. That mode's real job is Propose: it writes a first painting, which the
	// author then edits.
	Layer  string `json:"layer"`
	Legend string `json:"legend"`

	// Count is how many plates the lithosphere is in. Earth has seven or eight majors and a couple of dozen
	// minors; what matters here is that a boundary has to have room to be a mountain belt, so the useful
	// range on a hundred-kilometre planet is single digits.
	Count int `json:"count"`

	// SizeSpread is the ratio between the largest and smallest plate weight. The partition is a
	// multiplicatively weighted Voronoi, so a heavier site claims ground further away: 1 makes every plate
	// the same size, which is the one thing real plates never are.
	SizeSpread float64 `json:"size_spread"`

	// VelocityCmYr is how fast a plate moves, low to high. Earth runs 1 to 10; the number that matters
	// downstream is the *relative* speed across a boundary, which is a difference of two of these.
	VelocityCmYr [2]float64 `json:"velocity_cm_yr"`

	// SpinFraction is how much of a plate's speed is rotation about its own centre rather than translation.
	// Zero makes every margin uniform along its length, which is the defect the in-plane pole exists to
	// avoid; one makes the plate a pinwheel. A third of it is enough to turn a collision into a transform
	// over a few tens of kilometres.
	//
	// A pointer because zero is a real answer here and so is "say nothing". JSON cannot tell an absent
	// number from a zero one, and a plain float64 read them as the same thing - which is how a proposal came
	// out with every plate's spin at zero and every margin uniform, the one defect this field exists to
	// prevent. Absent takes the default; an explicit 0 means none.
	SpinFraction *float64 `json:"spin_fraction"`

	// WarpFraction is how far a boundary wanders from the straight Voronoi edge, as a fraction of the mean
	// plate spacing. Without it the partition is a polygon net and every margin is a ruled line.
	//
	// It is applied over two octaves, and that is not decoration either: one octave at the plate wavelength
	// gives a margin one long shallow bend, which at planet scale is still a ruled line with a kink in it.
	// The second octave at a third of the wavelength is what puts a promontory and a re-entrant into a
	// margin, and those are where a collision belt gets its along-strike segmentation from.
	WarpFraction float64 `json:"warp_fraction"`

	// ResolutionM is the tectonic grid's cell. See the package comment: coarse on purpose.
	ResolutionM float64 `json:"resolution_m"`

	// ContinentalFraction is the share of a plate's painted area that has to be land before it counts as
	// continental. Well below a half, because a continental plate carries a shelf and a passive margin as
	// well as its continent.
	ContinentalFraction float64 `json:"continental_fraction"`

	// ObliqueDeg is where a margin stops being convergent or divergent and becomes transform: the angle
	// between the relative velocity and the boundary normal, past which the strike-slip component is the
	// one in charge. 60 degrees means a margin stays convergent until the slip is over 1.7 times the
	// closing.
	ObliqueDeg float64 `json:"oblique_deg"`

	// Faults is the deformation zone around every margin: how wide it is and how densely it is broken. A
	// zero block means the margins carry no faults of their own, which is what every painted planet had
	// before it existed - the legend's per-class `faults` blocks are a separate, and now secondary, set.
	Faults Belt `json:"faults"`
}

// Default is the configuration a manifest that says nothing gets.
func Default() Config {
	return Config{
		Count:               7,
		SizeSpread:          1.7,
		VelocityCmYr:        [2]float64{1, 6},
		SpinFraction:        nil, // see Spin(); the default lives there so that an explicit 0 can mean none
		WarpFraction:        0.34,
		ResolutionM:         250,
		ContinentalFraction: 0.18,
		ObliqueDeg:          60,
	}
}

// DefaultSpinFraction is what a config that does not mention spin gets. It is not zero on purpose: a planet
// of plates that only drift has margins identical along their whole length, and the along-strike change from
// collision to transform is the most useful thing the model gives a fault set.
const DefaultSpinFraction = 0.35

// Spin is the configured spin fraction, or the default when the manifest said nothing. An explicit zero is
// honoured and means no rotation at all.
func (c Config) Spin() float64 {
	if c.SpinFraction == nil {
		return DefaultSpinFraction
	}
	if *c.SpinFraction < 0 {
		return 0
	}
	return *c.SpinFraction
}

func (c Config) withDefaults() Config {
	d := Default()
	if c.Count <= 0 {
		c.Count = d.Count
	}
	if c.SizeSpread < 1 {
		c.SizeSpread = d.SizeSpread
	}
	if c.VelocityCmYr[1] <= 0 {
		c.VelocityCmYr = d.VelocityCmYr
	}
	if c.VelocityCmYr[0] < 0 {
		c.VelocityCmYr[0] = 0
	}
	if c.WarpFraction < 0 {
		c.WarpFraction = d.WarpFraction
	}
	if c.ResolutionM <= 0 {
		c.ResolutionM = d.ResolutionM
	}
	if c.ContinentalFraction <= 0 {
		c.ContinentalFraction = d.ContinentalFraction
	}
	if c.ObliqueDeg <= 0 || c.ObliqueDeg >= 90 {
		c.ObliqueDeg = d.ObliqueDeg
	}
	return c
}

// Plate is one rigid piece of the lithosphere.
type Plate struct {
	ID int `json:"id"`

	// SiteXM, SiteYM is the Voronoi site in world metres, and Weight is what makes plates different sizes.
	SiteXM float64 `json:"site_x_m"`
	SiteYM float64 `json:"site_y_m"`
	Weight float64 `json:"weight"`

	// CentroidXM, CentroidYM is the plate's centre of area, measured the short way round the cylinder. It is
	// where the plate turns about, and it is the one position a painted plate has - a painting has no site.
	CentroidXM float64 `json:"centroid_x_m"`
	CentroidYM float64 `json:"centroid_y_m"`

	// The motion, in metres a year: a translation plus a rotation about a pole in the map plane. The pole is
	// always the centroid; it is stored rather than derived so that VelocityAt needs nothing but the plate.
	TransXM    float64 `json:"trans_x_m_yr"`
	TransYM    float64 `json:"trans_y_m_yr"`
	PoleXM     float64 `json:"pole_x_m"`
	PoleYM     float64 `json:"pole_y_m"`
	OmegaRadYr float64 `json:"omega_rad_yr"`

	// Continental is read from the painting rather than drawn from the seed: see the package comment.
	Continental  bool    `json:"continental"`
	LandFraction float64 `json:"land_fraction"`

	// AreaCells is the plate's size on the tectonic grid, which is what decides who overrides whom when two
	// oceanic plates converge.
	AreaCells int `json:"area_cells"`
}

// VelocityAt is the plate's motion at a world point, in metres a year.
//
// The lever arm is measured the short way round the cylinder. Without that a plate whose pole sits just east
// of the seam would spin the wrong way for every point just west of it, and the boundary running through the
// seam would be classified as convergent on one side and divergent on the other - the one bug this whole
// coordinate system exists to make impossible.
func (pl Plate) VelocityAt(p world.Planet, xM, yM float64) (vx, vy float64) {
	rx := wrapDelta(xM-pl.PoleXM, p.CircumferenceM())
	ry := yM - pl.PoleYM
	return pl.TransXM - pl.OmegaRadYr*ry, pl.TransYM + pl.OmegaRadYr*rx
}

// SpeedMYr is how fast the plate is going at its own site, which is the number worth printing.
func (pl Plate) SpeedMYr(p world.Planet) float64 {
	vx, vy := pl.VelocityAt(p, pl.SiteXM, pl.SiteYM)
	return math.Hypot(vx, vy)
}

// Model is a planet's tectonics: the plates, the grid they were rasterised on, and the boundaries between
// them.
type Model struct {
	P   world.Planet `json:"-"`
	Cfg Config       `json:"config"`

	Plates []Plate `json:"plates"`

	// The tectonic grid. GCellM is derived rather than taken: it is the circumference divided by a whole
	// number of columns, so the grid wraps exactly and a boundary crossing the seam is an ordinary one.
	GW     int     `json:"-"`
	GH     int     `json:"-"`
	GCellM float64 `json:"grid_cell_m"`

	// Cell is the plate id at every tectonic cell, row-major, X cyclic.
	Cell []int16 `json:"-"`

	// Boundaries is the whole planet's set, in world metres.
	Boundaries []Boundary `json:"boundaries"`
}

// GridXM and GridYM are the world position of a tectonic cell's centre. Y runs from the top of the polar
// pad, so a grid row and a planet row mean the same place.
func (m *Model) GridXM(gx int) float64 { return (float64(gx) + 0.5) * m.GCellM }
func (m *Model) GridYM(gy int) float64 { return m.P.YM(0) + (float64(gy)+0.5)*m.GCellM }

// GridIdx wraps X and clamps Y, the same way world.Planet.Idx does.
func (m *Model) GridIdx(gx, gy int) int {
	gx = ((gx % m.GW) + m.GW) % m.GW
	if gy < 0 {
		gy = 0
	} else if gy >= m.GH {
		gy = m.GH - 1
	}
	return gy*m.GW + gx
}

// PlateAt is which plate owns a world position.
func (m *Model) PlateAt(xM, yM float64) int {
	gx := int(math.Floor(xM / m.GCellM))
	gy := int(math.Floor((yM - m.P.YM(0)) / m.GCellM))
	return int(m.Cell[m.GridIdx(gx, gy)])
}

// Build draws a planet's plates and the boundaries between them, once, deterministically from the seed.
//
// land reports whether a world position is painted land. It is a callback rather than a raster so that this
// package knows nothing about templates: what it needs from the painting is one bit, and asking for it this
// way also means the caller decides how the land mask is sampled.
func Build(p world.Planet, seed int64, cfg Config, land func(xM, yM float64) bool) (*Model, error) {
	cfg = cfg.withDefaults()
	if cfg.Count < 2 {
		return nil, fmt.Errorf("a planet in %d plate(s) has no boundaries", cfg.Count)
	}

	circ := p.CircumferenceM()
	gw := int(circ/cfg.ResolutionM + 0.5)
	if gw < cfg.Count*4 {
		gw = cfg.Count * 4
	}
	gcell := circ / float64(gw)
	gh := int(float64(p.H)*p.CellM/gcell + 0.5)
	if gh < 2 {
		gh = 2
	}

	m := &Model{P: p, Cfg: cfg, GW: gw, GH: gh, GCellM: gcell}
	m.Plates = placeSites(p, seed, cfg, gh, gcell)
	giveMotion(p, seed, cfg, m.Plates)
	m.Cell = partition(m, seed)
	m.measure(land)
	m.Boundaries = m.buildBoundaries()
	return m, nil
}

// warpFineGain is how strong the second warp octave is against the first. Held well below a half: past that
// the displacement folds back on itself and the partition grows islands of one plate inside another, which
// is a nonsense the boundary tracer would faithfully chain into a ring.
const warpFineGain = 0.4

// spacingM is the mean distance between neighbouring sites: the length every other length in this package is
// a fraction of.
func spacingM(circ, heightM float64, count int) float64 {
	return math.Sqrt(circ * heightM / float64(count))
}

// placeSites scatters the plate centres, refusing any that lands on top of another.
//
// Rejection rather than relaxation. Lloyd's algorithm would give an even, hexagonal net, which is a worse
// answer than this one: plates are not even, and the interesting boundary geometry - a small plate wedged
// between two large ones, a long thin one - comes from exactly the irregularity relaxation removes. The
// minimum separation is only there to stop two sites coinciding, which produces a sliver no boundary tracer
// can chain.
func placeSites(p world.Planet, seed int64, cfg Config, gh int, gcell float64) []Plate {
	s := noise.NewSource(seed, srcSites)
	circ := p.CircumferenceM()
	heightM := float64(gh) * gcell
	top := p.YM(0)
	minSep := 0.5 * spacingM(circ, heightM, cfg.Count)

	out := make([]Plate, 0, cfg.Count)
	for len(out) < cfg.Count {
		for try := 0; ; try++ {
			x := s.Float() * circ
			y := top + s.Float()*heightM
			// After enough refusals the separation is the thing that is wrong, not the draw, so it is given
			// up rather than looped on for ever - a count near the area's limit can have no valid position
			// left at all.
			if try < 64 && tooClose(out, p, x, y, minSep) {
				continue
			}
			out = append(out, Plate{
				ID:     len(out),
				SiteXM: x,
				SiteYM: y,
				Weight: 1 + (cfg.SizeSpread-1)*s.Float(),
			})
			break
		}
	}
	return out
}

func tooClose(out []Plate, p world.Planet, x, y, minSep float64) bool {
	circ := p.CircumferenceM()
	for i := range out {
		dx := wrapDelta(x-out[i].SiteXM, circ)
		dy := y - out[i].SiteYM
		if math.Hypot(dx, dy) < minSep {
			return true
		}
	}
	return false
}

// giveMotion draws each plate's translation and its spin.
//
// **Every plate turns about its own centre of area**, which measure fills in once the partition exists. An
// earlier version put the pole a plate-width off to one side, on the reasoning that a pole at the middle
// would cancel symmetrically and leave the margins as uniform as a pure translation does. That reasoning is
// wrong: the relative velocity at a contact is
//
//	(T_a - T_b) + omega_a x (x - c_a) - omega_b x (x - c_b)
//
// which varies along the contact for any pole at all, and a pole at the centre puts the *largest* rotational
// contribution out at the margins, where it is wanted. The offset pole bought nothing and cost something
// real: it cannot be written into a painted legend, where an author says "this plate is also turning
// clockwise" and means about its own middle. A proposal therefore did not read back as the planet it was
// proposed from, which is what TestAProposalReadsBackAsTheSamePlanet caught.
func giveMotion(p world.Planet, seed int64, cfg Config, ps []Plate) {
	s := noise.NewSource(seed, srcMotion)
	circ := p.CircumferenceM()
	spacing := spacingM(circ, p.HeightM(), cfg.Count)
	for i := range ps {
		speed := s.Range(cfg.VelocityCmYr[0], cfg.VelocityCmYr[1]) / 100 // cm/yr to m/yr
		dir := s.Float() * 2 * math.Pi
		ps[i].TransXM = math.Cos(dir) * speed
		ps[i].TransYM = math.Sin(dir) * speed

		// The spin is set so that the rotational speed one spacing from the pole is SpinFraction of the
		// translation speed: the fraction means the same thing whatever the planet's size.
		sign := 1.0
		if s.Float() < 0.5 {
			sign = -1
		}
		if spacing > 0 {
			ps[i].OmegaRadYr = sign * cfg.Spin() * speed / spacing
		}
	}
}

// partition rasterises the plates onto the tectonic grid.
//
// A multiplicatively weighted Voronoi - nearest site by distance/weight - through a warped query point. The
// warp is what stops the result being a polygon net: it is sampled from a lattice in world coordinates, so
// two decompositions of the same planet warp the same point the same way, and its wavelength is deliberately
// long compared with the plate spacing, because a boundary that wiggled at a ten-kilometre wavelength would
// be a coastline rather than a plate margin.
func partition(m *Model, seed int64) []int16 {
	p := m.P
	circ := p.CircumferenceM()
	spacing := spacingM(circ, p.HeightM(), m.Cfg.Count)

	ws := noise.NewSource(seed, srcWarp)
	// Two octaves. Both cell counts are whole numbers of the noise period, because noise.Lattice.Sample wraps
	// modulo its own count and anything else is a discontinuity down one meridian.
	coarse := int(p.NoisePeriodM/spacing + 0.5)
	if coarse < 1 {
		coarse = 1
	}
	fine := coarse * 3
	wx := noise.NewLattice(coarse, ws)
	wy := noise.NewLattice(coarse, ws)
	fx := noise.NewLattice(fine, ws)
	fy := noise.NewLattice(fine, ws)
	amp := m.Cfg.WarpFraction * spacing

	out := make([]int16, m.GW*m.GH)
	for gy := 0; gy < m.GH; gy++ {
		yM := m.GridYM(gy)
		v := yM / p.NoisePeriodM * float64(coarse)
		fv := yM / p.NoisePeriodM * float64(fine)
		row := gy * m.GW
		for gx := 0; gx < m.GW; gx++ {
			xM := m.GridXM(gx)
			u := xM / p.NoisePeriodM * float64(coarse)
			fu := xM / p.NoisePeriodM * float64(fine)
			qx := xM + ((float64(wx.Sample(u, v))*2-1)+(float64(fx.Sample(fu, fv))*2-1)*warpFineGain)*amp
			qy := yM + ((float64(wy.Sample(u, v))*2-1)+(float64(fy.Sample(fu, fv))*2-1)*warpFineGain)*amp

			best, bestID := math.Inf(1), 0
			for i := range m.Plates {
				pl := &m.Plates[i]
				dx := wrapDelta(qx-pl.SiteXM, circ)
				dy := qy - pl.SiteYM
				d := math.Hypot(dx, dy) / pl.Weight
				if d < best {
					best, bestID = d, pl.ID
				}
			}
			out[row+gx] = int16(bestID)
		}
	}
	return out
}

// measure walks the partition once and fills in everything that can only be known after it exists: each
// plate's area, how much of it the author painted as land, and where its centre is.
//
// Painted rows only, for the land fraction. The polar pad is synthetic ocean that no class was ever painted
// on, so counting it would drag every plate that reaches a pole towards oceanic for a reason that is
// scaffolding rather than geography. The area and the centre do count the pad, because a plate really does
// extend over it.
//
// **The centre is a circular mean in X.** A plate painted or drawn across the meridian has cells at both ends
// of the raster, and an arithmetic mean of those columns puts its centre on the far side of the planet - and
// with it the pole the whole plate rotates about, which would make its velocity field nonsense and every
// margin around it wrong. This is the same rule as world.WrapX and Plate.VelocityAt's lever arm: the short
// way round is the only way round.
func (m *Model) measure(land func(xM, yM float64) bool) {
	n := len(m.Plates)
	landCells := make([]int, n)
	paintedCells := make([]int, n)
	sumSin := make([]float64, n)
	sumCos := make([]float64, n)
	sumY := make([]float64, n)

	circ := m.P.CircumferenceM()
	for gy := range m.GH {
		yM := m.GridYM(gy)
		painted := yM >= 0 && yM < m.P.HeightM()
		row := gy * m.GW
		for gx := range m.GW {
			id := int(m.Cell[row+gx])
			pl := &m.Plates[id]
			pl.AreaCells++

			xM := m.GridXM(gx)
			ang := 2 * math.Pi * xM / circ
			sumSin[id] += math.Sin(ang)
			sumCos[id] += math.Cos(ang)
			sumY[id] += yM

			if !painted {
				continue
			}
			paintedCells[id]++
			if land != nil && land(xM, yM) {
				landCells[id]++
			}
		}
	}

	for i := range m.Plates {
		pl := &m.Plates[i]
		if paintedCells[i] > 0 {
			pl.LandFraction = float64(landCells[i]) / float64(paintedCells[i])
		}
		pl.Continental = pl.LandFraction >= m.Cfg.ContinentalFraction
		if pl.AreaCells > 0 {
			pl.CentroidXM, pl.CentroidYM = centroid(sumSin[i], sumCos[i], sumY[i], pl.AreaCells, circ)
		}
		// Every plate turns about its own centre of area. See giveMotion for why it is not somewhere else.
		pl.PoleXM, pl.PoleYM = pl.CentroidXM, pl.CentroidYM
	}
}

// centroid turns the accumulated sums into a position, taking X the short way round the cylinder.
func centroid(sumSin, sumCos, sumY float64, cells int, circ float64) (xM, yM float64) {
	ang := math.Atan2(sumSin, sumCos)
	if ang < 0 {
		ang += 2 * math.Pi
	}
	return ang / (2 * math.Pi) * circ, sumY / float64(cells)
}

// wrapDelta brings a difference in X into -circ/2 .. +circ/2: the short way round the cylinder.
func wrapDelta(d, circ float64) float64 {
	if circ <= 0 {
		return d
	}
	d = math.Mod(d, circ)
	if d > circ/2 {
		d -= circ
	} else if d < -circ/2 {
		d += circ
	}
	return d
}