A Row of Springs Makes a Wave

// internal/field/spring.go

// Spring pulls a body back toward one place, and pulls harder the
// further the body has been taken from it. K is how much force one
// world pixel of displacement is worth.
type Spring struct {
	Rest Vec2
	K    float64
}

// Force is this spring's pull on one body: the displacement from rest,
// turned around and multiplied by K. A body sitting at rest is pulled
// by nothing.
func (s Spring) Force(b Body) Vec2 {
	return b.Pos.Sub(s.Rest).Scale(-s.K)
}

// Damper is a force that points against whatever a body is doing and
// grows with how fast it is doing it. C is how much force one world
// pixel per tick of speed is worth.
type Damper struct {
	C float64
}

// Force is the drag this damper puts on one body: its velocity, turned
// around and multiplied by C.
func (d Damper) Force(b Body) Vec2 {
	return b.Vel.Scale(-d.C)
}

// Stable reports whether a spring this stiff can be carried by a body of
// the given mass one whole tick at a time. Past this bound the update's
// own arithmetic runs away, whatever the spring is attached to.
func (s Spring) Stable(mass, damping float64) bool {
	return s.K > 0 && s.K < 4*mass-2*damping
}

Two types in twenty-odd lines, and neither has heard of the other. A body picks up both pushes the same way it picks up wind and weight, in whatever order the caller happens to write them, and the accumulator does not care which is which:

// cmd/waves/main.go — the whole of an oscillator's tick

func (o *oscillator) step() {
	o.body.ApplyForce(o.spring.Force(o.body))
	o.body.ApplyForce(o.damper.Force(o.body))
	o.body.Step()
}

That is the composability the accumulator was built for, cashed in a third time. A spring and a damper are two independent laws with two independent authors, meeting in a vector neither of them owns.

$ go run ./cmd/waves -mode one -ticks 8
waves: one body of mass 1.0 on a spring, k 0.10, damping 0.00
       rest at y 0.000, released from y 12.000, stable bound k < 4.00
  tick       force       accel         vel         pos
     1      -1.200      -1.200      -1.200      10.800
     2      -1.080      -1.080      -2.280       8.520
     3      -0.852      -0.852      -3.132       5.388
     4      -0.539      -0.539      -3.671       1.717
     5      -0.172      -0.172      -3.843      -2.125
     6       0.213       0.213      -3.630      -5.755
     7       0.576       0.576      -3.054      -8.810
     8       0.881       0.881      -2.173     -10.983
period: 19.786 ticks predicted, 19.786 measured over 20 swings

The first four ticks are the hand table, to the last decimal place. Tick 5 is the overshoot: the position goes negative, which means the body has crossed its rest height and is now above it, and the force turns positive on tick 6 without anything checking for the crossing. Nothing in Force tests a sign. Subtracting the rest place from the position gives a negative displacement on its own, and negating a negative gives a downward pull.

Two numbers on the last line agree that were computed in completely different ways. One comes from the table above, out of an angle and a division. The other comes from running the body for twenty full swings and timing the ticks where it crosses rest going up. When a prediction from the arithmetic and a measurement from the run agree to three decimals, the model of the thing and the thing are the same model.

One stiffness deserves a hand check because it comes out in whole numbers. At k = 1 the angle per tick is a sixth of a turn, so a swing takes exactly six ticks, and if the release is 12 the body only ever visits three heights:

$ go run ./cmd/waves -mode one -k 1 -ticks 7
waves: one body of mass 1.0 on a spring, k 1.00, damping 0.00
       rest at y 0.000, released from y 12.000, stable bound k < 4.00
  tick       force       accel         vel         pos
     1     -12.000     -12.000     -12.000       0.000
     2       0.000       0.000     -12.000     -12.000
     3      12.000      12.000       0.000     -12.000
     4      12.000      12.000      12.000       0.000
     5       0.000       0.000      12.000      12.000
     6     -12.000     -12.000       0.000      12.000
     7     -12.000     -12.000     -12.000       0.000
period: 6.000 ticks predicted, 6.000 measured over 20 swings

12, 0, −12, −12, 0, 12, and back to where it started on tick 7. Every entry is exact in binary, so this run repeats forever with no drift at all, and it is the cheapest possible check that a spring you have just written is wired up the right way round. If your version prints 24 on tick 1, the sign is inverted and you have built something that runs away from rest.

$ go run ./cmd/waves -mode damp -k 0.1 -c 0.05 -ticks 600
waves: one release from y 12.000, k 0.10, run twice
       free: no damper.  held: damping 0.05
  tick      free y  free swing      held y  held swing
   100      10.672      24.190       0.690       3.863
   200       8.121      24.164       0.014       0.302
   300       4.639      24.161      -0.002       0.024
   400       0.626      24.213      -0.000       0.002
   500      -3.458      24.254      -0.000       0.000
   600      -7.147      24.282      -0.000       0.000
after 600 ticks the free spring still swings 24.282 and the held one 1.006e-05
numbers: free 2b46f27575442e2d  held a2264cafdbf979df

Sixty seconds of world time. The free spring is still going 24 pixels top to bottom on the last tick, as far as it ever was, because nothing in the update takes energy out of it: a force goes in, comes back out on the return, and the books balance. The swing column wanders between 24.16 and 24.29 only because a swing lasting 19.786 ticks is never sampled at quite the same place twice, so the highest sample is rarely the highest point.

The held column is the same spring with one extra caller. Ten seconds in it has lost five sixths of its travel, twenty seconds in it is inside a third of a pixel, thirty seconds in it is inside a fortieth, and by the end it is a hundred-thousandth of a pixel from rest, which on a screen is a body that has stopped. The damper never decides to stop it. It only ever subtracts a fraction of whatever speed it finds, and a fraction of a fraction of a fraction goes to nothing on its own.

Both runs end with a hash over every position and velocity of all six hundred ticks, so the claim being made here is not "it looked settled". Run either one again and the sixteen characters come back the same, on any machine that adds floats the way this one does. A world that damps its water differently on Tuesday is a world whose replay is worthless.

// cmd/waves/main.go

// row is many oscillators side by side, each one primed with the state
// the reference oscillator held some whole number of ticks ago.
type row struct {
	osc     []oscillator
	off     []int
	hist    []state
	k, c, a float64
}

// newRow builds cols oscillators and hands column i the reference
// oscillator's state from off(i) ticks into its history. A column whose
// offset is bigger than the history is left standing at rest.
func newRow(k, c, amp float64, cols int, off func(int) int) *row {
	max := 0
	offs := make([]int, cols)
	for i := range offs {
		offs[i] = off(i)
		if offs[i] > max {
			max = offs[i]
		}
	}
	r := &row{osc: make([]oscillator, cols), off: offs, hist: history(k, c, amp, max), k: k, c: c, a: amp}
	for i := range r.osc {
		r.osc[i] = newOsc(k, c, 0)
		if o := offs[i]; o >= 0 {
			r.osc[i].put(r.hist[o])
		}
	}
	return r
}

func (r *row) step() {
	for i := range r.osc {
		r.osc[i].step()
	}
}

$ go run ./cmd/waves -mode strip -ticks 6
waves: the first eight of 32 springs, k 0.10, 1 ticks of offset per column
  tick     col 0     col 1     col 2     col 3     col 4     col 5     col 6     col 7
     0   -10.176   -11.742   -12.133   -11.311    -9.358    -6.470    -2.934     0.895
     1    -7.593   -10.176   -11.742   -12.133   -11.311    -9.358    -6.470    -2.934
     2    -4.250    -7.593   -10.176   -11.742   -12.133   -11.311    -9.358    -6.470
     3    -0.483    -4.250    -7.593   -10.176   -11.742   -12.133   -11.311    -9.358
     4     3.333    -0.483    -4.250    -7.593   -10.176   -11.742   -12.133   -11.311
     5     6.816     3.333    -0.483    -4.250    -7.593   -10.176   -11.742   -12.133
     6     9.617     6.816     3.333    -0.483    -4.250    -7.593   -10.176   -11.742

Put a finger on −10.176 in the top left and move it down and right. It is in column 0 at tick 0, column 1 at tick 1, column 2 at tick 2, and it keeps going to the edge of the pond. So does −12.133, the deepest point of that swing, two columns to its east. So does every other number on the grid. The surface pattern moves east at one column per tick, which on four-pixel columns is forty pixels a second.

Nothing is moving east. Read the table the other way, down a single column, and column 0 goes −10.176, −7.593, −4.250, −0.483, 3.333: that body is climbing, and climbing is all it will ever do apart from falling. It has no sideways velocity, no neighbour it talks to, no idea that a crest went past it. The travelling pattern is the offsets, and the offsets were fixed once when the row was built and have not changed since.

// cmd/waves/main.go

// drawPool paints the valley through the camera, then cuts the pond's
// waterline down to where the row of springs is standing.
func drawPool(b *render.Buffer, m *render.Tiles, cam render.Camera, r *row) {
	b.Fill(render.Void)
	vis := cam.Cells(m.T, m.Cols, m.Rows)
	for cy := vis.Y0; cy < vis.Y1; cy++ {
		for cx := vis.X0; cx < vis.X1; cx++ {
			sx, sy := cam.ToScreen(cx*m.T, cy*m.T)
			b.FillRect(render.Rect{X0: sx, Y0: sy, X1: sx + m.T, Y1: sy + m.T}, m.At(cx, cy))
		}
	}
	// The pond's still surface sits half way down its top row; every
	// column of it moves to wherever its own spring is.
	rest := PondY0*m.T + m.T/2
	for i := range r.osc {
		y := rest + int(math.Round(r.height(i)))
		x0 := PondX0*m.T + i*Column
		bank := m.At(PondX0+i*Column/m.T, PondY0-1)
		sx, _ := cam.ToScreen(x0, 0)
		_, top := cam.ToScreen(0, PondY0*m.T)
		_, sy := cam.ToScreen(0, y)
		b.FillRect(render.Rect{X0: sx, Y0: top, X1: sx + Column, Y1: sy}, bank)
		b.FillRect(render.Rect{X0: sx, Y0: sy, X1: sx + Column, Y1: sy + 2}, foam)
		b.FillRect(render.Rect{X0: sx, Y0: sy + 2, X1: sx + Column, Y1: sy + 4}, trough)
	}
}

$ go run ./cmd/waves -mode pool -amp 6 -ticks 60
waves: 32 springs across the pond, k 0.10, damping 0.00
       1 ticks of offset per column, 19.786 ticks a swing, 79.1 pixels a wave
       the pond is 128 pixels across and a crest crosses it at 4.0 pixels a tick
  tick     col 0     col 3     col 6     col 9         frame           numbers
     0    -5.088    -5.656    -1.467     3.956  098a306e6dfd  8599352473c7e8c2
    10     4.972     5.728     1.667    -3.796  dfb9333c871b  aa10974075cf75c3
    20    -4.851    -5.794    -1.864     3.633  1ff2d35b82ee  1f3db960dc05638d
    30     4.723     5.853     2.060    -3.465  ced8d5011214  6245986d203a1a8e
    40    -4.591    -5.905    -2.253     3.293  54493026ba47  6a52d733519bea87
    50     4.453     5.950     2.444    -3.118  4cf64777076b  aee2163803dd2159
    60    -4.309    -5.989    -2.631     2.939  129dad5edded  2e18c0563bd9c194
every column is the one recorded history, offset: true

The picture behind those hashes is the pond with a curve for a waterline: a bright two-pixel crest, a darker band under it, and dry bank above wherever the surface has dropped. This release is 6 pixels instead of the strip run's 12, so the waterline stays inside the pond's top row of tiles. A crest and a half fit across the water, they slide east at four pixels a tick, and the surface never stops, because no damper was asked for.

The last line is the interesting one. It is not a description; it is a question the program asks itself every run. Take the reference oscillator that was recorded before any column existed, run it forward to tick 60 plus this column's own offset, and compare the two numbers to the last bit. Thirty-two columns, sixty ticks of independent stepping, and every one of them is still exactly the recording, offset. That is what "nothing travels" means when a machine says it.

The two knobs on the surface do different jobs, and mixing them up wastes an afternoon. The offset per column sets the speed: one tick per column moves the crest four pixels a tick, three ticks per column slows it to 1.3. The stiffness sets how many ticks a swing takes, and the two together set the distance between crests. Turn k up and the water gets choppier without getting any faster:

$ go run ./cmd/waves -mode pool -k 0.6 -amp 6 -ticks 60
waves: 32 springs across the pond, k 0.60, damping 0.00
       1 ticks of offset per column, 7.899 ticks a swing, 31.6 pixels a wave
       the pond is 128 pixels across and a crest crosses it at 4.0 pixels a tick
  tick     col 0     col 3     col 6     col 9         frame           numbers
     0     6.488    -5.070     0.893     3.769  13ae7fce74fd  7aa734a9bf7e7ca8
    10    -0.145     4.566    -6.503     4.902  00471ef5dfd5  17dd84b18a6dcffa
    20    -6.459     4.158     0.405    -4.748  3a925b60161f  7558112f78f7f10c
    30     1.435    -5.397     6.422    -3.954  e5b676a9950e  270fc805e666b315
    40     6.173    -3.080    -1.688     5.538  194319d9cbbf  9d77e33ca4b47ab1
    50    -2.668     6.012    -6.085     2.849  991aea7e1cb8  572e10a8e16abf7f
    60    -5.640     1.880     2.903    -6.107  c0c919e2af9c  f37895636ec699e6
every column is the one recorded history, offset: true

31.6 pixels between crests instead of 79.1, so four of them fit across the pond and the waterline reads as chop instead of a swell. The crest speed on the third header line has not moved: still four pixels a tick, because the offsets did not change.

// cmd/waves/main.go — the offset is the only line that changed

	across := render.MapCols * render.TileSize / Column
	down := render.MapRows * render.TileSize / Column
	r := newRow(k, c, amp, across*down, func(i int) int {
		dx, dy := float64(i%across-hitX), float64(i/across-hitY)
		return -int(math.Round(math.Hypot(dx, dy) * perPatch))
	})

// cmd/waves/main.go — a patch joins in when the front reaches it

		for i := range r.osc {
			if -r.off[i] == t {
				r.osc[i].put(r.refAt(0)) // the front has just arrived here
				continue                 // and its first tick is the next one
			}
			r.osc[i].step()
		}

$ go run ./cmd/waves -mode ripple -c 0.05 -amp 8 -ticks 48 -shot ripple.png
waves: 23040 patches of ground, each on its own spring, k 0.10, damping 0.05
       struck at patch 136,48, the front crossing one patch per 1 tick
  tick     front   at blow   8 across  20 across         frame           numbers
     0         0     8.000      0.000      0.000  9fd0141a7752  c94fb67673c2016c
     8        32    -5.752      8.000      0.000  e73682eb7ccd  84eec5467cf771f4
    16        64     2.549     -5.752      0.000  f6e16110ffef  5d6368deda55e9ce
    24        96     0.334      2.549      1.489  6b8164cc405a  25072d7027efb206
    32       128    -2.148      0.334     -4.164  6aa808f3d9e8  a76db1ff5b1b32dd
    40       160     2.712     -2.148      4.699  a399b1059b1a  277c15d1f71af2b2
    48       192    -2.279      2.712     -3.656  943eac1e875f  6fdc06db0c9c9c76
every patch the front has reached is the one recorded history, offset: true
wrote ripple.png

The "at blow" and "8 across" columns are the strip table again, in a different disguise. The patch eight across holds 8.000 at tick 8, which is what the struck patch held at tick 0; it holds −5.752 at tick 16, which is what the struck patch held at tick 8. The column labelled twenty across is flat zeros until the front reaches it on tick 20, and then it starts the same story two seconds later than the middle.

The picture is rings. Each patch is drawn a little lighter or darker than the ground under it, by how far its spring has carried it, so a crest reads as a pale circle and a trough as a dark one, and they walk outward at four pixels a tick. The damper is what makes it a ripple instead of a permanent target painted on the field: the middle has almost given up by tick 48, while the outer ring is still going, because the outer ring started later and has had fewer ticks to lose.

// internal/field/spring_test.go

// TestARowIsOneHistoryOffset is the wave, written down as arithmetic.
// Column i is primed with the state the reference oscillator held i
// ticks into its own run, and from then on the two are stepped
// separately. If a wave is nothing but offsets, the two must agree bit
// for bit forever.
func TestARowIsOneHistoryOffset(t *testing.T) {
	const (
		k, c, amp = 0.1, 0.02, 12
		cols      = 16
		ticks     = 400
	)
	ref := newSwing(k, c, amp)
	hist := []struct{ Pos, Vel Vec2 }{{ref.b.Pos, ref.b.Vel}}
	for range ticks + cols {
		ref.step()
		hist = append(hist, struct{ Pos, Vel Vec2 }{ref.b.Pos, ref.b.Vel})
	}

	row := make([]*swing, cols)
	for i := range row {
		row[i] = newSwing(k, c, 0)
		row[i].b.Pos, row[i].b.Vel = hist[i].Pos, hist[i].Vel
	}
	for n := 1; n <= ticks; n++ {
		for i := range row {
			row[i].step()
			if row[i].b.Pos != hist[n+i].Pos || row[i].b.Vel != hist[n+i].Vel {
				t.Fatalf("at tick %d column %d holds %v %v, but the reference at tick %d held %v %v",
					n, i, row[i].b.Pos, row[i].b.Vel, n+i, hist[n+i].Pos, hist[n+i].Vel)
			}
		}
	}
}

$ go test -count=1 -v ./internal/field/
=== RUN   TestPushesSumThenDivide
--- PASS: TestPushesSumThenDivide (0.00s)
=== RUN   TestOnePushTwoMasses
--- PASS: TestOnePushTwoMasses (0.00s)
=== RUN   TestGravityIgnoresMass
--- PASS: TestGravityIgnoresMass (0.00s)
=== RUN   TestSofteningCapsThePull
--- PASS: TestSofteningCapsThePull (0.00s)
=== RUN   TestPullIsProportionalToDisplacement
--- PASS: TestPullIsProportionalToDisplacement (0.00s)
=== RUN   TestDamperOpposesTheBody
--- PASS: TestDamperOpposesTheBody (0.00s)
=== RUN   TestARowIsOneHistoryOffset
--- PASS: TestARowIsOneHistoryOffset (0.00s)
=== RUN   TestStiffnessPastFourTimesTheMassRunsAway
    spring_test.go:112: k 4.1 passed a million pixels from rest on tick 32
--- PASS: TestStiffnessPastFourTimesTheMassRunsAway (0.00s)
PASS
ok  	theworld/internal/field	0.002s

Four hundred ticks and sixteen columns is 6,400 comparisons of two floats each, and an equality test on floats is the right test here precisely because it is the brittle one. The two sides are not approximations of one truth; they are the same operations on the same bits in the same order, so a single reordering anywhere in Step would move one of them and not the other. The last test is the failure box kept: 3.9 stays on the field for two thousand ticks, 4.1 passes a million pixels on tick 32, and the log line records which tick it went.