359 lines
15 KiB
Go
359 lines
15 KiB
Go
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
|
|
}
|