Wind Is a Function of Place

// internal/field/lay.go

// Lay is a wind written down cell by cell: one force for every cell of
// the world's grid, and an answer for any position in between. It is the
// first thing in this world that is a function of place rather than a
// number: hand it where a body is standing and it hands back what the
// air there is doing to it.
type Lay struct {
	Cols, Rows int
	T          float64 // cell size, in world pixels

	cell []Vec2
}

// NewLay writes down one force per cell by asking at for each of them.
// The grid is walked in row order, so two runs build the same table.
func NewLay(cols, rows int, t float64, at func(cx, cy int) Vec2) *Lay {
	l := &Lay{Cols: cols, Rows: rows, T: t, cell: make([]Vec2, cols*rows)}
	for cy := 0; cy < rows; cy++ {
		for cx := 0; cx < cols; cx++ {
			l.cell[cy*cols+cx] = at(cx, cy)
		}
	}
	return l
}

// Cell is the force written down for one cell. Coordinates off the grid
// read as the nearest edge cell, so a body drifting past the rim keeps
// getting an answer instead of an index panic.
func (l *Lay) Cell(cx, cy int) Vec2 {
	cx = clamp(cx, 0, l.Cols-1)
	cy = clamp(cy, 0, l.Rows-1)
	return l.cell[cy*l.Cols+cx]
}

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

// At is the force at any position on the field: the four cell centres
// nearest that position, each weighted by how near it is. The four are
// added in one fixed order, so the answer is the same on every run.
func (l *Lay) At(p Vec2) Vec2 {
	u := p.X/l.T - 0.5
	v := p.Y/l.T - 0.5
	i, j := int(math.Floor(u)), int(math.Floor(v))
	fx, fy := u-float64(i), v-float64(j)

	tl := l.Cell(i, j).Scale((1 - fx) * (1 - fy))
	tr := l.Cell(i+1, j).Scale(fx * (1 - fy))
	bl := l.Cell(i, j+1).Scale((1 - fx) * fy)
	br := l.Cell(i+1, j+1).Scale(fx * fy)
	return tl.Add(tr).Add(bl).Add(br)
}

The two lines that subtract a half are the only awkward part. A cell's force belongs at its centre, and cell 2's centre is at pixel 40, not 32, so dividing the position by the cell size gives the wrong origin by half a cell. Subtract the half and the whole number part names the cell whose centre is up and to the left, while the fraction says how far past that centre the body has got. Those two fractions are the weights: at half a cell past in both directions all four weights are 0.25, and directly on a centre one weight is 1 and the rest are 0.

Four scaled vectors are then added in the order the code lists them, top row before bottom row, which is the same rule the accumulator lives by. A sum whose order is decided anywhere else is a sum the seed does not control.

// cmd/wind/main.go

// base is the wind out on open soil: half a world pixel per tick per
// tick, blowing east.
var base = field.Vec2{X: 0.5, Y: 0}

// exposure is how much of the base wind reaches one kind of ground: open
// water has nothing standing in the way, soil has the valley floor's
// scrub, and a rock cell is the rim itself.
func exposure(t sim.Terrain) float64 {
	switch t {
	case sim.Water:
		return 1.6
	case sim.Soil:
		return 1.0
	}
	return 0.2
}

// windAt is the force written down for one cell: the base wind scaled by
// what the ground there lets through, plus a spill away from every rock
// neighbour, because air that cannot go through a wall goes around it.
func windAt(g *sim.Grid, cx, cy int) field.Vec2 {
	here, err := g.At(sim.Coord{X: cx, Y: cy})
	if err != nil {
		return field.Vec2{}
	}
	w := base.Scale(exposure(here))
	for _, d := range []sim.Coord{{X: -1}, {X: 1}, {Y: -1}, {Y: 1}} {
		n := sim.Coord{X: cx + d.X, Y: cy + d.Y}
		t, err := g.At(n)
		if err != nil || t == sim.Rock {
			w = w.Add(field.Vec2{X: -float64(d.X), Y: -float64(d.Y)}.Scale(0.15))
		}
	}
	return w
}

// lay builds the whole wind from a world: one force per cell, sampled at
// any position by the four cells nearest it.
func lay(g *sim.Grid) *field.Lay {
	return field.NewLay(g.W, g.H, Tile, func(cx, cy int) field.Vec2 {
		return windAt(g, cx, cy)
	})
}

$ go run ./cmd/wind -mode lay
wind: a 12x8 valley from seed 5, base 0.500,0.000
      ground                 wind, in pixels per tick per tick
      ############   0.10 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.10
      #..........#   0.25 0.67 0.52 0.52 0.52 0.52 0.52 0.52 0.52 0.52 0.38 0.05
      #..........#   0.25 0.65 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.35 0.05
      #.~~~~~....#   0.25 0.65 0.80 0.80 0.80 0.80 0.80 0.50 0.50 0.50 0.35 0.05
      #.~~~~~~...#   0.25 0.65 0.80 0.80 0.80 0.80 0.80 0.80 0.50 0.50 0.35 0.05
      #...~~~~...#   0.25 0.65 0.50 0.50 0.80 0.80 0.80 0.80 0.50 0.50 0.35 0.05
      #...~......#   0.25 0.67 0.52 0.52 0.81 0.52 0.52 0.52 0.52 0.52 0.38 0.05
      ############   0.10 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.18 0.10

      one transect east along row 3, through the pond and out the far shore
      position       cell     the field's answer
      24.000,56.000   1,3 soil   0.650,0.000
      32.000,56.000   2,3 water  0.725,0.000
      40.000,56.000   2,3 water  0.800,0.000
      104.000,56.000  6,3 water  0.800,0.000
      112.000,56.000  7,3 soil   0.650,0.000
      120.000,56.000  7,3 soil   0.500,0.000

      the four cells nearest 32.000,32.000, and what each is worth
      cell 1,1  centre   24.0,24.0  wind 0.650,0.150  weight 0.25
      cell 2,1  centre   40.0,24.0  wind 0.500,0.150  weight 0.25
      cell 1,2  centre   24.0,40.0  wind 0.650,0.000  weight 0.25
      cell 2,2  centre   40.0,40.0  wind 0.500,0.000  weight 0.25
      the field's answer at 32.000,32.000: 0.575,0.075

The pond stands out at 0.80 against the floor's 0.50, and the column beside the eastern rim reads 0.05, the quietest air in the valley: the wall's spill points west and takes 0.45 out of the 0.50 arriving there. Along the top and bottom the spill points inward, so row 1 reads 0.52 where row 2 reads 0.50: the extra 0.02 is all a sideways 0.15 adds to the length of a vector already 0.50 long. The table prints lengths, and the vectors under it point in different directions.

The transect is the blend at work. At pixel 24 the body is on cell 1's centre and gets cell 1's 0.650 exactly. At pixel 40 it is on cell 2's centre and gets 0.800. At pixel 32, halfway between, it gets 0.725, and the cell column says it is standing in cell 2 while half of what it feels comes from cell 1. The pond's edge is a cliff in the ground and a ramp in the wind, which is what a leaf crossing the shore should feel.

The last block is a position on no centre at all. Four cells, each a quarter, two of them carrying a downward 0.150 spilled off the rim above, and the answer 0.575, 0.075 is a vector that appears in none of them. That number is the point of the whole arrangement: the field has more places in it than it has cells.

// cmd/wind/main.go — the tick

for i := range bodies {
	f := l.At(bodies[i].Pos)
	bodies[i].ApplyForce(f)
	bodies[i].Step()
}

$ go run ./cmd/wind -mode valley -ticks 60 -drag=false
wind: 16 bodies on the seed 5 valley, base 0.500,0.000
      no drag at all
  tick  on the field  fastest    mean         frame           numbers
     0            16    0.000   0.000  6f2c67e383c5  a0dada96a2cd0116
    20            12   18.905  10.027  701b90a34739  54dfb866ccb7ac01
    40             4   16.905  10.989  9a26cdc5ef49  663c70af9eac8bc8
    60             0   14.974  10.225  89930306f1ad  ab562de6eb15dca8

Every one of them is gone inside six seconds, and the field is 192 pixels wide. The fastest body was crossing at nearly 19 pixels a tick, which is more than a cell every tick. Nothing here is broken. A force applied over and over with nothing opposing it builds speed without limit, exactly as it did when the only force was gravity, and a wind that never stops pushing is a wind that eventually throws everything it touches off the edge of the world. Real air prevents that by pushing back on whatever is moving through it.

// internal/field/drag.go

// Drag is what a body pushing through air or water gets back: a force
// against the velocity, growing in proportion to it. B is how thick the
// stuff is; nothing here scales with the body's size, because nothing in
// this world has one yet.
type Drag struct{ B float64 }

// Force is the drag on one body this tick.
func (d Drag) Force(b Body) Vec2 { return b.Vel.Scale(-d.B) }

// Terminal is the speed at which this drag exactly cancels a force of
// the given size. Above it the drag wins and the body slows; below it
// the force wins and the body speeds up.
func (d Drag) Terminal(force float64) float64 { return force / d.B }

$ go run ./cmd/wind -mode drift -ticks 120
wind: a steady 0.500,0.000 on masses 0.5 and 4.0, drag B 0.25
      predicted ceiling for both: 2.000 pixels per tick
  tick   v(m=0.5)   drag(m=0.5)     v(m=4)     drag(m=4)
     1     1.0000       -0.0000     0.1250       -0.0000
     2     1.5000       -0.2500     0.2422       -0.0312
     3     1.7500       -0.3750     0.3521       -0.0605
     4     1.8750       -0.4375     0.4550       -0.0880
     5     1.9375       -0.4688     0.5516       -0.1138
    10     1.9980       -0.4990     0.9511       -0.2203
    20     2.0000       -0.5000     1.4499       -0.3533
    40     2.0000       -0.5000     1.8487       -0.4596
    80     2.0000       -0.5000     1.9886       -0.4969
   120     2.0000       -0.5000     1.9991       -0.4998
      after 120 ticks: 2.000000 and 1.999134
      they are 1.000 and 1.000 of the way to the ceiling

Two bodies eight times apart in mass, one steady wind, and both of them stop speeding up at 2 pixels a tick. The light one is inside a thousandth of it by tick 10 and the heavy one needs the whole 120 to get that close, but the number they approach is the same, and the drag column says why: at 2.0 the drag reads 0.5000 against a wind of 0.5000, so the accumulator receives two forces that add to nothing and the velocity has no reason to change again.

That is a different run from one with wind alone. With no drag, the same wind keeps pushing the lighter body ahead faster forever. Here the ratio closes. Weight scales with mass and drag does not, so drag is the force that eventually makes every body agree.

The -0.0000 in the first row is not a misprint. The drag is computed from the velocity at the start of the tick, and at the start of tick 1 both bodies were standing still, so the force is zero scaled by minus a quarter. Floating point keeps a negative zero and prints it.

// cmd/wind/main.go

// fall drops the same two bodies through still air, where the force that
// drives them is their own weight and so is not the same for both.
func fall(ticks int) {
	leaf := field.NewBody(field.Vec2{}, 0.5)
	stone := field.NewBody(field.Vec2{}, 4)
	...
	for t := 1; t <= ticks; t++ {
		for _, b := range []*field.Body{&leaf, &stone} {
			b.ApplyForce(b.Weight(gravity))
			b.ApplyForce(air.Force(*b))
			b.Step()
		}
	}
}

$ go run ./cmd/wind -mode fall -ticks 120
wind: dropped in still air, gravity 0.000,0.500, drag B 0.25
      mass 0.5: weight 0.250, so the ceiling is 1.000 pixels per tick
      mass 4.0: weight 2.000, so the ceiling is 8.000 pixels per tick
  tick   v(leaf)   fell    v(stone)   fell
     1    0.5000    0.5      0.5000    0.5
     2    0.7500    1.2      0.9688    1.5
     3    0.8750    2.1      1.4082    2.9
     4    0.9375    3.1      1.8202    4.7
     8    0.9961    7.0      3.2262   15.6
    16    1.0000   15.0      5.1514   50.7
    32    1.0000   31.0      6.9857  151.2
    64    1.0000   63.0      7.8714  393.9
   120    1.0000  119.0      7.9965  840.1
      the leaf is at 1.00000 of its ceiling, the stone at 0.99957

On tick 1 they fall at the same rate, because both start at rest and drag on something at rest is nothing. From tick 2 they separate, and after twelve seconds the stone has fallen 840 pixels and the leaf 119. With no air in the calculation the two fall in step; with air, the leaf falls slowly because the air pushes back hard for the leaf's mass. Here is that sentence as a column of numbers. The stone's weight is eight times the leaf's and the air treats their velocities identically, so the stone needs eight times the speed before the drag can match it.

// cmd/wind/main.go

var (
	air  = field.Drag{B: 0.25} // how thick the air is
	pond = field.Drag{B: 1.25} // and the water, five times over
)

// dragAt is which stuff a body at p is moving through. It is the second
// field in this program and the plainer kind: one value per cell, read
// from the cell the body stands in and not blended with its neighbours.
func dragAt(g *sim.Grid, wet bool, p field.Vec2) field.Drag {
	if !wet {
		return air
	}
	c := sim.Coord{X: int(p.X) / Tile, Y: int(p.Y) / Tile}
	if t, err := g.At(c); err == nil && t == sim.Water {
		return pond
	}
	return air
}

Two fields now, and they are not the same kind of thing. The wind is a vector per cell and gets blended, because a leaf feels a little more air as it approaches open water. The drag is one number per cell and gets read whole, because the leaf is either on the water or off it and there is no halfway. Which one a field needs is a question about the quantity, not about the code: a blend is right when the thing being sampled really does vary smoothly, and wrong when it names a state.

// internal/field/drag.go

// Force is the drag on one body this tick. The last two lines are the
// cap: a drag may bring a body to rest inside one tick and may not push
// it the other way, whatever B and the mass work out to.
func (d Drag) Force(b Body) Vec2 {
	f := b.Vel.Scale(-d.B)
	if stop := b.Mass * b.Vel.Len(); f.Len() > stop {
		f = b.Vel.Unit().Scale(-stop)
	}
	return f
}

// Bare is the same drag with the cap taken off: the arithmetic straight
// out of the definition, kept so the two can be run against each other.
func (d Drag) Bare(b Body) Vec2 { return b.Vel.Scale(-d.B) }

$ go run ./cmd/wind -mode valley -ticks 10 -wet -watch 4 -from 1 -to 10
wind: 16 bodies on the seed 5 valley, base 0.500,0.000
      air B 0.25, pond B 1.25, capped
watching body 4, mass 0.5
  tick  on the field  fastest    mean         frame           numbers
     0            16    0.000   0.000  6f2c67e383c5  a0dada96a2cd0116
t=  1  B 1.25  wind 0.742,0.000  drag -0.000,-0.000  vel 1.483,0.000  pos 35.267,66.402
t=  2  B 1.25  wind 0.756,0.000  drag -0.742,-0.000  vel 1.511,0.000  pos 36.778,66.402
t=  3  B 1.25  wind 0.770,0.000  drag -0.756,-0.000  vel 1.540,0.000  pos 38.317,66.402
t=  4  B 1.25  wind 0.784,0.000  drag -0.770,-0.000  vel 1.568,0.000  pos 39.886,66.402
t=  5  B 1.25  wind 0.799,0.000  drag -0.784,-0.000  vel 1.598,0.000  pos 41.484,66.402
t=  6  B 1.25  wind 0.800,0.000  drag -0.799,-0.000  vel 1.600,0.000  pos 43.084,66.402
t=  7  B 1.25  wind 0.800,0.000  drag -0.800,-0.000  vel 1.600,0.000  pos 44.684,66.402
t=  8  B 1.25  wind 0.800,0.000  drag -0.800,-0.000  vel 1.600,0.000  pos 46.284,66.402
t=  9  B 1.25  wind 0.800,0.000  drag -0.800,-0.000  vel 1.600,0.000  pos 47.884,66.402
t= 10  B 1.25  wind 0.800,0.000  drag -0.800,-0.000  vel 1.600,0.000  pos 49.484,66.402

The same body, the same seed, the same pond, and it crosses the water at a steady 1.600 pixels a tick. Every drag reading from tick 2 onward is the previous tick's wind with its sign flipped, which is the cap doing precisely what it says: it removes the speed the wind put on and nothing further.

The cap is not free, and it is important to be exact about what it cost. Where it binds, this body is no longer obeying the drag law at all. The law's ceiling over the pond would be 0.800 divided by 1.25, which is 0.64 pixels a tick. What the run settles at is 1.600, which is 0.800 divided by the mass of 0.5. A cap that binds every tick has quietly replaced the physics with itself. It is a guard against nonsense, not a model, and a coefficient that keeps it engaged permanently is a coefficient chosen wrongly.