Porting a tide engine to Swift with the original as the test oracle

TL;DR — Porting numeric code by hand and eyeballing the output is how you
ship plausible-looking wrong answers. Instead, treat the original
implementation as a test oracle: scripts that run the reference and dump
golden fixtures, one per algorithm layer, each ported layer gated on matching
the oracle to 1e-6 before the next begins. And when the port passes, you're
still not done — the last gate is matching reality, not the reference.
Repo: slackwater-engine.

We needed a pure-Swift harmonic tide/current prediction engine — tide math is
deterministic astronomy, so given a station's harmonic constituents you can
compute heights and high/low turns for any minute, years ahead, with zero
network. The best open reference implementation is
Neaps
(@neaps/tide-predictor, JavaScript, MIT). This post is about how we ported
it — the oracle-driven method — because the method transfers to any port of
numeric code between languages.

The problem: numeric ports fail silently

The Neaps algorithm is ~600 lines of real math: astronomical mean-longitude
polynomials, Doodson coefficients, Schureman/IHO node-correction formulas
full of terms like this one:

private func _nupp(_ N: Double, _ i: Double, _ omega: Double) -> Double {
    let I = d2r * _I(N, i, omega), nu = d2r * _nu(N, i, omega)
    let tan2 = pow(sin(I), 2) * sin(2 * nu) / (pow(sin(I), 2) * cos(2 * nu) + 0.0727)
    return r2d * 0.5 * atan(tan2)
}

Flip a sign in one of those, or mix degrees and radians once, and the engine
still produces beautiful sinusoidal tide curves — highs and lows in roughly
the right places, just quietly off by centimetres or minutes in ways that
drift with the 18.6-year lunar nodal cycle. Eyeballing plots against a tide
table catches none of it. A port validated by "looks right" is a liability
you discover at a shallow anchorage.

Wire the original up as an oracle

The fix is to make the reference implementation generate the tests. A small
Node script imports Neaps, evaluates each layer of the algorithm at fixed
instants, and writes the results as JSON fixtures the Swift tests load:

// tools/gen-golden.mjs — Neaps is the oracle. TideEngine (Swift) must match.
import { astro, constituents, createTidePredictor } from '@neaps/tide-predictor';

// Spread of UTC timestamps across the 18.6-year nodal cycle to exercise corrections.
const TIMES = [
  '2000-01-01T00:00:00Z', '2010-06-15T12:00:00Z', '2020-03-21T06:30:00Z',
  '2026-07-12T00:00:00Z', '2026-07-12T18:45:00Z', '2031-11-02T09:15:00Z',
  '2035-09-23T00:00:00Z', '2040-12-31T23:00:00Z',
];

write('astronomy.json', {
  note: 'Neaps astro(date) unwrapped to plain degrees; TideEngine must match within 1e-6.',
  times: TIMES,
  values: TIMES.map((t) => ({ time: iso(t), astro: unwrap(astro(new Date(t))) })),
});

Two properties make this an oracle rather than just "some tests":

1. Fixtures are generated, never hand-written. node tools/gen-golden.mjs regenerates every fixture from @neaps/tide-predictor@0.10.0. There is no opportunity to transcribe an expected value wrong, and bumping the oracle version re-derives the whole suite.
2. There's a fixture per layer, not just end-to-end. An end-to-end height comparison tells you that you diverged; a per-layer fixture tells you where. Every divergence is caught at the layer that caused it.

Port layer by layer, gate each layer

The algorithm decomposes into five layers, each consuming the one below.
Each got the same treatment: generate the fixture → write the failing Swift
test → port the math → gate on green before starting the next layer.

Layer	Golden check	Tolerance
Astronomy	15 quantities (mean longitudes + node angles) × 8 instants across the nodal cycle — 120 values	1e-6°
Node corrections	IHO f/u for 17 base constituents × 3 instants	1e-6
Constituents	V₀ + compound f/u for ~39 constituents (incl. M4, MS4, 2MK3)	1e-6
Timeline prediction	48 h height series, mixed-tide station, 10-min step	1e-6 m
Extremes	high/low count, kind, time, height	60 s / 0.02 m

Because both sides implement the same algorithm, the interesting tolerance is
floating-point agreement, not "close enough for tides." The prediction test
asserts both — the domain tolerance per point, and a divergence tripwire on
the max error:

var maxErr = 0.0
for (got, expected) in zip(points, fixture.points) {
    let err = abs(got.height - expected.height)
    maxErr = max(maxErr, err)
    #expect(err < 0.02, "height at \(expected.time): got \(got.height), want \(expected.height)")
}
// The engines share the algorithm, so agreement should be far tighter than 2cm.
#expect(maxErr < 1e-6, "max height error \(maxErr) exceeds 1e-6 — engines diverging")

If a future change drifts the port to, say, 5 mm of error — still a perfectly
plausible tide — the tripwire fires long before the domain tolerance would.

The per-layer gates earn their keep on the fiddly semantic mismatches between
the languages, because each one surfaces in the one layer it belongs to:

• Timeline construction. Neaps floors the start and ceils the end to the step interval. Get that wrong and every sample is shifted — the astronomy and constituent layers stay green, and the prediction fixture pinpoints it.
• Array identity vs value semantics. Neaps recomputes node corrections every 24 h and callers detect the refresh with a JS !== identity check. Swift arrays are value types — there is no identity. The port adds an explicit generation counter that callers compare instead. The extremes fixture is what proves the refresh logic actually matches.
• What not to port. Neaps ships both Schureman and IHO node-correction schemes; the predictor only ever uses IHO. The oracle made that safe to verify — no fixture ever exercised Schureman, so it stayed unported.

Don't port the catalog — generate it

The single biggest porting risk wasn't the math — it was the data. Neaps
bundles ~5,000 lines of constituent definitions: 394 constituents with
speeds, Doodson coefficients, and compound-member decompositions, plus an
IHO Annex-B name parser that resolves compound names at runtime. Hand-porting
5,000 lines of numeric literals is transcription, and transcription is where
silent errors live.

So the catalog isn't code in the Swift port at all. A second codegen script
walks Neaps' constituent map and emits it as a bundled JSON resource — with
the name parsing already resolved:

// tools/gen-catalog.mjs — the IHO Annex-B name decomposition + sign resolution
// runs HERE, in Neaps, at build time; Swift consumes the resolved members
// and never needs the parser.
for (const key of Object.keys(constituents)) {
  const c = constituents[key];
  if (key !== c.name) aliases[key] = c.name;   // NOAA "NU2" → canonical "nu2"
  seen.set(c.name, {
    name: c.name,
    speed: c.speed,
    coefficients: c.coefficients ?? null,      // Doodson, when fundamental
    members: c.members?.map((x) => ({ name: x.constituent.name, factor: x.factor })),
  });
}

A compound constituent in the emitted catalog.json (52 KB, 394 constituents,
83 name aliases) looks like this:

{"name": "2MK3", "speed": 42.9271398,
 "coefficients": [3, -1, 0, 0, 0, 0, 1],
 "members": [{"name": "M2", "factor": 2}, {"name": "K1", "factor": -1}]}

And the entire Swift side of "compound constituent support" collapses to two
recursions over pre-resolved members:

/// Node correction (f, u°). IHO fundamental if the constituent has one, else
/// combined from members: f = Π f_memberᵃᵇˢ⁽ᶠᵃᶜᵗᵒʳ⁾, u = Σ factor·u_member.
func correction(_ name: String, _ a: Astro) -> (f: Double, u: Double) {
    let key = canonical(name)
    if let fundamental = ihoCorrection(key, a) { return fundamental }
    guard let e = byName[key] else { return (1, 0) }
    var f = 1.0, u = 0.0
    for m in e.members ?? [] {
        let c = correction(m.name, a)
        u += m.factor * c.u
        f *= pow(c.f, abs(m.factor))
    }
    return (f, u)
}

The aliases matter more than they look: real station data is inconsistent
about constituent names (NOAA publishes NU2, MM, RHO; the catalog's
canonical names are nu2, Mm, rho1), and the engine silently ignores
constituents it can't resolve. Without alias resolution, published harmonic
constants would predict — just with pieces of the tide missing.

Net result: the whole engine is 577 lines of Swift against Neaps'
~6,300-line bundle, because the 5,000 data lines crossed over as data.

Golden checks → swift test → CI

The first iterations ran as a throwaway command-line checker. Once the layers
were green, the checks converted to swift-testing
suites — fixtures under Tests/TideEngineTests/Fixtures/ — so the whole
oracle runs as plain swift test, and CI is 18 lines:

jobs:
  test:
    name: swift test
    runs-on: macos-15
    steps:
      - uses: actions/checkout@v4
      - name: Run tests
        run: swift test

That conversion is worth calling out as its own step: an oracle that lives in
a side script rots; an oracle that is the test suite re-verifies the entire
port on every push, forever.

Matching the reference isn't the finish line

At this point the Swift engine reproduces Neaps to floating-point agreement.
That proves the port is faithful. It proves nothing about whether the
algorithm predicts actual water — if the reference had a bug, we ported the
bug at 1e-6 fidelity.

So the final gate compares against the tide authority itself. A third codegen
script takes published harmonic constants for Friday Harbor, WA
(NOAA 9449880, via @neaps/tide-database) and fetches NOAA's own official
high/low predictions from the CO-OPS API for the same window:

const fh = stations.find((s) => s.id === 'noaa/9449880');
const offset = fh.datums.MSL - fh.datums.MLLW; // shift MSL-relative harmonics to chart datum

const url = `https://api.tidesandcurrents.noaa.gov/api/prod/datagetter?...`
  + `&station=9449880&product=predictions&datum=MLLW&interval=hilo&units=metric&time_zone=gmt`;

The test predicts extremes from the constants and matches each official
high/low to the nearest computed one of the same kind — with deliberately
looser tolerances than the oracle tests, because now the differences are
real (datum rounding, node-correction epoch, constituent-set differences
between NOAA's internal engine and ours):

#expect(timeErr < 15, "\(off.kind) at \(off.time): time off by \(timeErr) min")
#expect(heightErr < 0.15, "\(off.kind) at \(off.time): height off \(heightErr) m")

Result across 12 highs/lows over three days: max time error 7.9 minutes,
max height error 3.5 cm — well inside tide-table tolerance. The ~8-minute
residual is expected and instructive: NOAA's engine uses a different
node-correction epoch and constituent set, so exact agreement with the
authority was never on the table. Two tolerance regimes, two different
claims: 1e-6 says "the port is faithful"; ±15 min/±0.15 m says "the
predictions are real."

The same validation pass confirmed the data path for Canadian waters:
CHS doesn't publish bundled harmonics the way NOAA does, but its IWLS API
serves tide predictions directly — and current predictions for 35+ BC
stations — so Canadian coverage rides the authority's own numbers online
rather than redistributed constituents.

Trim the API before anyone depends on it

Porting layer-by-layer leaves scaffolding: everything was public so the
interim checker could poke at it, and some ported code turned out dead (an
astronomical-speeds derivation the engine never uses — constituent speeds
ship pre-computed in the catalog). The last commit of the phase is a
deliberate audit: delete the dead code, demote the internals, and ship the
smallest surface that does the job —

public struct Station {
    public init(constituents: [HarmonicConstituent], offset: Double = 0)
    public func heights(from: Date, to: Date, step: TimeInterval = 600) -> [TidePoint]
    public func extremes(from: Date, to: Date) -> [TideExtreme]
}

Four public types, two methods. Astro, the catalog, the node-correction
table — all internal. Every symbol you expose during the port is a symbol
you can't rename after someone imports it; the oracle makes the trim safe,
because swift test proves the demotions broke nothing.

Three transferable lessons

1. An oracle turns a risky port into a mechanical one. Generate golden vectors from the reference per layer, and every divergence is caught at the layer that caused it — sign flips, degree/radian mixups, and value-vs-reference semantic mismatches all surfaced as a single red layer, never as "the tides look slightly wrong."
2. Generate data catalogs from the source; only port algorithms. The 5,000-line constituent catalog crossed languages as a JSON resource with the hard parts (name parsing, compound decomposition) pre-resolved at codegen time. Data you generate can't be mistranscribed, and regenerating beats re-porting when upstream moves.
3. Matching the reference isn't the finish line — matching reality is. A faithful port of a wrong reference is still wrong. Close the loop against the authority's published numbers, with honest (looser) tolerances and an explanation for the residual.

The engine is MIT, extracted as the open-source core for offline tide work on
our own boats: slackwater-engine.
If your home waters disagree with it, the fixtures show you exactly which
layer to blame — send a fix.

Related: Why Claude won't transcribe your PDF — and what to do instead · Why generic weather MCPs fail for marine navigation (use NDBC buoys)