The simulator

Point-scatterer and Fourier-wavefield backends, scenes, heightfields, and objects.

A GPU-accelerated synthetic aperture sonar data simulator. Produces POSSM-compatible HDF5 files that the beamform_tdbp pipeline consumes directly.

What You Get Out

run_simulation(...) writes a single HDF5 at output_path with the layout expected by beamformer/sas/data_loader.py:

/ReceiverInformation (attrs: SampleRateHZ, CenterFrequencyHZ) and ElementPositionsM/{x,y,z}
/TransmitterInformation/Transmitter1/ElementPositionsM/{x,y,z}
/SimEnvironment attr waterSpeedMetersPerSecond
/Pings/<N> (1-based) with Data (complex baseband, shape (n_channels, n_samples)), sensorPosition, timeSeconds, and {roll,pitch,heading}Degrees

The baseband data is directly usable by the beamformer; no post-hoc conversion required.

Two Backends

run_simulation(scene, motion, cfg, out, method='ps' | 'fw', ...) dispatches to one of:

Point-Scatter (`method='ps'`, default)

Per-ping CUDA pipeline:

Upload scatterer positions + reflectivities + heightfield + normals to GPU once (invariant across pings).
For each ping, compute TX/RX positions in NED, optionally include intra-ping platform velocity and body-frame angular velocity for per-channel motion during TOF.
gpu_raytrace produces delays, amplitudes, visibility (ray-heightfield LOS shadows), and grazing angle per scatterer per channel.
gpu_echo accumulates complex echoes into a (n_channels, n_samples) time series via an atomicAdd CUDA kernel convolved with the LFM chirp.

Accurate but slow: roughly a minute per ping at 1000 scatterers/m^2. Requires CuPy. Supports shadows. shadow_enabled=False disables the LOS check.

Fourier Wavefield (`method='fw'`)

Implements Sanford et al. 2024 (references/sanford2024.pdf) in the paper's four steps:

A. Geometric rendering: the seafloor heightfield is rendered once from nadir (world grid, Sanford Eq. 11); each aspect (spacing fw_aspect_spacing_deg, default 1 deg since 2026-09-10; the regression ladder pins 0.5) rotates the normals (Eq. 4) and runs a per-aspect line-of-sight occlusion pass sheared along y tan(phi) with the true per-pixel grazing angle. Placed CAD meshes are rendered per aspect with a tilted orthographic camera along the sonar ray, first hit per pixel with the facet normal (simulator/fw_mesh.py, FW_MESH_FACETS), so facing walls, self-occlusion and layover come from the geometry; the mesh's surface samples are then not injected as points. Procedural boxes, cylinders, wedges and truncated cones register their own closed mesh and take the same path. Whatever is still injected (cables, manual points, porous meshes) is occluded per aspect against the render's horizon (FW_INJECT_OCCLUDE).
B. Scattering intensity: Lambert from the normal maps with the same BSS calibration as PS, projected onto the slant-range plane; mesh facets add (TS / A_mesh) cos(angle) per camera pixel.
C. Wavefield: random phase screen, 2-D FFT, Stolt mapping to (k_q, ω), system window, inverse FFT. The result is ONE monostatic, stop-and-hop, straight-track wavefield e(q, t). The Stolt gather interpolates along k_r (Lanczos-5), so the range FFT is padded to keep the far edge of the swath inside the kernel's passband (_compute_fw_pad_dims, docs/fw_calibration_tests.md fix 10); a 100 m swath at 18.75 mm needs 16384 range bins.
D. Data generation: per ping and per receiver: gather e at the phase-centre along-track coordinate (FFT upsample fw_upsample_factor, default 8), then apply one exact delay map dτ(ping, ch, t) holding the bistatic excess, sway/heave with range-dependent grazing, attitude and (unless simulation.stop_and_hop) intra-pulse drift (fw_engine.compute_step_d_delays). There is no post-processing after the HDF5 is written; the stop_and_hop attribute is the config value, as for PS.

Not modelled: the angle dependence of the bistatic excess and the off-broadside (Doppler) part of drift / wide-beam sway; both are applied at their broadside value. Validation ladder and tolerances: docs/fw_calibration_tests.md. Design: docs/superpowers/specs/2026-09-09-fw-step-d-exact-design.md.

fw_pixel_size <= lambda/4 (~1.25 mm at 300 kHz).

Both backends write identical HDF5 layout.

Coordinates and the Altitude Rule

NED (x=North / along-track, y=East / cross-track, z=Down). The heightfield stores positive-up elevations; when converted to NED scatterer positions, z = -elevation. Platform altitude above the seafloor is -pos_z; P_0 = [0, 0, -altitude[0]].

run_simulation(..., max_range=...) is a slant range in meters. By convention, SAS platform altitude is roughly 10% of max ground range. E.g. ground range 60 m → altitude 6–8 m → slant range ~60 m.

Design-Rule Warnings and Preflights (2026-09-05)

simulator/scenes/run.py::instantiate runs three checks before any heavy work, on every path (CLI, cloud worker, studio, GUI), not only in the GUI validation layer:

Design rules (simulator/design_rules.py, the formulas of docs/reference/sas_sonar.md): R2 PRF vs range ambiguity, R6 altitude band, R7 along-track advance per ping ≤ L/2, R8 channel spacing vs λ/2, R17 TX azimuth width ≥ channel spacing (ghost rule), R18 baseband Nyquist, R19 elevation null. Printed as [sim-prep] design rule ... lines; never blocks a run. APERTURE_DESIGN_RULE_WARNINGS=0 silences them.
Host-RAM preflight (mem_preflight.host_memory_plan): predicts the whole run's peak (build, simulation, beamforming) and compares it with MemAvailable less a 3 GB reserve (APERTURE_RAM_RESERVE_GB). A run that fits runs in RAM. One that does not goes out of core (simulator/spill.py, 2026-09-30): every full-length scatterer array (positions, reflectivities, cell and zone indices, the beam-tile sort's keys and sorted copy) becomes a memory map over an unlinked file in APERTURE_SPILL_DIR (default ~/.cache/aperturelab/spill; keep it on a local Linux disk, not a Windows drive under WSL). RAM then holds only the grids, the per-thread block work, the object scatterers and the CUDA / beamformer base, a few GB at any density, and the run needs ~42 B per scatterer of free disk. The output is the same as the in-RAM run (the same operations on the same blocks; test_spill.py). Measured on the benchmark scene at 20,000/m^2: 3.7 GB anonymous peak against 11.3 GB in RAM, 176 s against 161 s. MemoryError only when neither fits. APERTURE_SPILL=0 refuses instead of spilling, =1 always spills; APERTURE_SKIP_RAM_PREFLIGHT=1 bypasses the check.
GPU preflight (mem_preflight.preflight_gpu_ram, PS backend, just before the first device upload): the resident elevation + normals grids (32 B/cell) plus the echo accumulators are checked against the device total (raise) and what is free right now (warn). Same bypass flag.
Object placement contract (run.validate_objects, also run by the studio's validate_scene): every key on an object dict must be a parameter of that kind's placer or one of OBJECT_METADATA_KEYS (kind, scour, annotation, name, id, _scatterer_range); PLACERS in run.py is the single kind -> placer -> renames table. A bad key is a one-line objects[i] (kind 'name'): unexpected key ... error at prep, never a TypeError inside placement (2026-09-05: the studio's name identity key killed a run that way). Guard tests: simulator/tests/test_object_placement_contract.py tie the studio's documented object keys and the GUI's dropped keys to the signatures.
Object burial (scour.burial_fraction): with burial_is_height_fraction: false (the default) it is the fraction of the equilibrium scour-pit depth the object settles into, so it needs enabled: true and does nothing on mud or clay, neither of which scour in this version. With the flag true it is the fraction of the object's own height that is buried, measured from the bed after the object's own pit is dug, with or without scour, and the object's only sink; the pit is still dug around it but does not add to it. Full rules, the seating reference and the label wiring: realistic_simulation_notes.md, the seating section.

Related loader rule: a bottom.zone_grid_path that does not exist beside the YAML is a FileNotFoundError; the inline zone_grid next to it is only a 1x1 stub, and the old silent fallback to it built zone-less scenes. APERTURE_ALLOW_MISSING_ZONE_GRID=1 restores the fallback (with a warning) for deliberate use.

Scatterer Density

Scene-level parameter (Scene(scatterer_density=...), units scatterers/m^2). Controls speckle fidelity and GPU memory:

Density	Use
`200`	Debug. Debug at 200 first.
`2000`	Standard production runs.
`5000`	High-fidelity speckle / publication figures.
`33000`	Physics-accurate but memory-intensive.

Memory scales linearly: a 3000 m^2 flat_seafloor at 2000/m^2 yields ~6 M scatterers (~3 GB at 32 channels). Iterate at 200/m^2; only bump density when you know the setup is right.

Workflow

from simulator import Scene, PlatformMotion, SonarConfig, run_simulation

cfg = SonarConfig()                               # 300 kHz / 75 kS/s / 60 kHz BW defaults
motion = PlatformMotion.straight_line(
    speed=1.5, altitude=8.0, n_pings=128, ping_rate=3.3,
    sway_amplitude=0.15, sway_period=10.0,        # optional 6-DOF perturbations
    roll_amplitude=0.01,
)

scene = Scene(
    x_min=-5, x_max=65, y_min=0, y_max=70,
    cell_size=0.05, bottom_type='sand',
    scatterer_density=200.0,                       # debug first!
)
scene.set_bottom_zone_map(                         # optional zoned seafloor
    [['sand'], ['rock']],
    rms_heights={'sand': 0.001, 'rock': 0.1},
    seed=42,
)
scene.add_procedural_rocks(                        # optional rock population
    zone=(1, 0), density=2.0, d_min=0.2, d_max=1.0, seed=543,
)
scene.add_box_object(                              # objects raise the heightfield
    x_center=20, y_center=20,
    x_extent=2.0, y_extent=1.0, height=0.5, ts_db=-15.0, seed=43,
)
scene.build(seed=42)                               # must be called before run_simulation

run_simulation(scene, motion, cfg,
               output_path='sim_output/scene.h5',
               method='ps', max_range=60.0,
               shadow_enabled=True, verbose=True)

The Heightfield

2D elevation grid (n_x, n_y) with cell_size. Owns:

elevation[n_x, n_y]: positive-up heights (m). Objects raise local cells to max(existing, object_top).
bottom_type (uniform) or bottom_type_grid (per-cell index into _zone_bottom_types) when SeafloorBuilder runs a zone map.
object_scatterers: explicit list of (x, y, z_ned, complex_refl) tuples added by place_*.
Surface normals via compute_normals() with a 60 deg max-slope clamp (prevents single-cell object edges from producing near-horizontal normals that silence seafloor scatterers).

Primitives on the heightfield (all exposed via Scene.add_*):

place_box_object(x_center, y_center, x_extent, y_extent, height, ts_db)
place_cylinder_object(x_center, y_center, length, diameter, orientation='along_track'|'cross_track', ts_db, heading_deg=None, proud=False): horizontal cylinder. The axis sits on the datum (one radius stands proud) unless heading_deg is given, which turns the stamp, scatterers and FW facets together to any heading; proud=True with it rests the body on the bed, one diameter high
place_rockan_mine(x_center, y_center, heading, ts_db): tapered wedge (1.02 x 0.50–0.80 x 0.385 m, from simulator/object_sizes.py)
place_manta_mine(x_center, y_center, ts_db): truncated cone (0.98 m base x 0.44 m tall, from simulator/object_sizes.py; cylinder size presets for mine-like bodies live in the same module, see scenes-framework.md, Published object sizes)
place_mesh_object(x_center, y_center, mesh_path, scale, rotation_z, ts_db) : loads STL/OBJ/PLY via trimesh, rasterizes into heightfield via downward ray shooting, then samples surface scatterers uniformly
place_superellipsoid_rock(a, b, c, n_shape, ts_db, ...) ; |x/a|^n + |y/b|^n + |z/c|^n = 1. n_shape<2 pinched, >2 boxy
place_procedural_rocks(x_min, x_max, y_min, y_max, density, d_min, d_max, alpha) : Poisson-placed superellipsoids, log-normal axis ratios, size-dependent n_shape (large rocks angular, small rocks rounded)
place_faceted_rock(x_center, y_center, a, b, c, n_shape, planes, ts_db, ...): convex superellipsoid clipped by half spaces (planes rows nx, ny, nz, d in the body frame); six vertical sides plus a top make a basalt column. sink_m on both rock placers buries the body: the footprint shrinks, buried scatterers are dropped, the FW mesh is clipped at the bed.
add_point_scatterer(x, y, z, ts_db, n_sub=10, jitter_m=0.01) ; Rayleigh sub-scatterers inside a small jitter disc (avoids the single-pixel coherent halo)

Large-scale bathymetry (2026-09-21)

Everything else on the heightfield is small-scale relief: the fractal background tops out around 8 cm rms, ripples at a few centimetres. The bottom.bathymetry operator list adds the scene-scale relief: the regional slope and bedform field a real survey line crosses: up to bottom.bathymetry_clip_m (default 5 m).

procedural_features.compose_bathymetry(ops, n_x, n_y, cell_size, clip_m=5.0) sums the operators, clips the sum to +/- clip_m and makes it zero-mean (clip and demean iterate to a fixed point, so the returned field satisfies both bounds at once). The result is added to hf.elevation as the first relief stage in scenes/run.py::instantiate(), before the zone warp and every other operator, so ripples, pockmarks, trawl scars, outcrops, rocks and objects all settle on the bathymetry. Because the field is zero-mean the bottom changes shape without changing depth, so the altitude rule (altitude ~ 10 % of max ground range) still holds. The realised span is stored as hf.bathymetry_range_m ((0.0, 0.0) when no operators are configured).

`kind`	Function	Parameters	What it models
`rolling`	`fractal_background`	`rms_m`, `correlation_length_m`, `spectral_exponent`, `seed`	Goff-Jordan `1/k^beta` rolling relief at the scene scale. Same generator as the small-scale background, driven an order of magnitude harder and with an outer scale of tens of metres.
`tilt`	`bathymetry_tilt`	`along_deg`, `across_deg`	Planar regional slope, zero at the grid centre (so it adds no mean depth). `across_deg > 0` deepens the bed toward the far range.
`sand_waves`	`bathymetry_sand_waves`	`wavelength_m`, `amplitude_m`, `direction_deg`, `asymmetry`, `envelope_fraction`, `seed`	A migrating sand-wave train: profile `sin(phi) + asymmetry * sin(2 phi) / 2` with `phi = 2 pi (x cos(theta) + y sin(theta)) / wavelength`, rescaled to a peak of `amplitude_m`. `asymmetry` steepens one flank and flattens the other (lee / stoss) without skewing the crest heights.

envelope_fraction < 1 turns the wave train into a sand-wave field covering that fraction of the scene: smoothed noise thresholded at the 1 - envelope_fraction quantile (correlation length 2 wavelengths), passed through zone_blend.inward_feather with the same sigma. The feather is inward, so the amplitude is exactly zero at the field boundary and rises from there: a step is impossible: and the outline is organic, never rectangular (realistic_simulation_notes.md, sections 2 and 3). A field that runs off the edge of the scene is not feathered there: inward_feather counts pixels beyond the array as inside, so the bedforms continue past the scene instead of fading at its border.

The envelope's correlation length is hundreds of pixels on a production grid, so the noise is synthesised on a decimated grid (working sigma ~4 px) and bilinearly zoomed back; the threshold and the feather still run at full resolution, so the mask edge is pixel-exact. Deterministic in seed, but the realisation differs from a full-resolution smooth.

save_dem_png(hf, path, meta=None) in scenes/run.py writes the resulting heightfield as a 16-bit DEM: v = round((z - z_min) / (z_max - z_min) * 65535), with the decode constants in the PNG's aperturelab tEXt chunk (beamformer/sas/imaging.py::read_png_meta). Rows are flipped like every other image product, so the file is in the house orientation (along-track UP, range RIGHT, row 0 = x_min at the image bottom); a reader that wants the heightfield's array order applies [::-1, :] after decoding. It is the data counterpart of save_heightmap_png, a normalised 8-bit hillshade preview in the same orientation that cannot be decoded back to metres.

The SeafloorBuilder

Driven by Scene.set_bottom_zone_map(grid, rms_heights, seed) and Scene.add_procedural_rocks(...). During scene.build():

Divide the heightfield XY extent into equal-area rectangular zones, one per cell of grid (rows index along-track).
Synthesize ONE global FFT power-law roughness realisation and composite it per pixel through feathered per-zone amplitude maps (no per-zone FFTs, so no seams). Each type is in one of two spectrum modes (see "Roughness spectrum" below): absolute, the surface is W(k) = b k^-γ with b = spectral_strength and no rescale (built-in types), or rms, the legacy realisation normalised to rms_height_m (custom types by default, and any type named in rms_heights).
Cells already occupied by placed objects are protected via protect_mask so object tops stay smooth.
Rocks auto-populate by bottom_type.default_rock_density unless disabled; explicit add_procedural_rocks overrides apply last. An override carrying boulder=BoulderParams (a preset or any knob on the YAML rock group) is sampled by simulator/boulder_field.py instead: lithology, Zingg shape family, Thomas clusters or an edge band, burial and optional scour. See boulder_fields.md.

Optional Tang ripple overlays (scene.set_tang_ripples(ripple_wavelength, rms_height, direction_deg, skewness_b)) layer on top: skewness_b>0 uses the Tang eq. (5) exponential transform for sharp peaks / wide troughs (b=18 m^-1 matches Traykovski data).

Ripples from wave forcing

A custom type's ripple: block may carry a forcing: sub-block (hs_m, tp_s, depth_m, d50_mm, wave_dir_deg, r_reset, r_washout) and an age_days. simulator/ripple_forcing.py resolves them at prep time, in simulator/scenes/run.py's custom-type loop, so CLI, GUI, Studio and cloud all see the same result:

State from the Penko & Kearney ripple reset parameter Λ = Hs² / (4 h d50) (TREX13 thresholds 50 / 150, site-tuned): relict below r_reset, active between, washout at or above r_washout.
Active: wavelength, RMS height and crest direction are replaced by the Nelson et al. 2013 equilibrium for that sea state (λ/A = 1/[0.72 + 2.0e-3 Δ (1 − e^{−(1.57e-4 Δ)^1.15})], Δ = A/d50, η/λ = 0.12 λ^−0.056; rms = η/(2√2); direction_deg = wave_dir_deg). 0.5 m / 6 s waves in 7.5 m over 0.23 mm sand give 0.26 m ripples, as observed at TREX13.
Relict: the explicit geometry is kept (it is the author's record of what the last event left). Washout: no overlay, the type shows its base roughness only; design rule R23 warns.
Age: bioturbation decay exp(−D k0² t), D = 5e-9 m²/s (NSEA), scales the RMS height and, as a heuristic for crest rounding, the skewness. A 25 cm ripple e-folds in ~4 days, a 1.5 m one in ~130.

The run log prints one [sim-prep] ripple '<type>': <state> ... line per forced or aged block. References: references/penko_kearney_ripple_reset_icce.pdf, references/kearney_penko_2022_nsea.pdf; the Nelson coefficients were taken from secondary sources (Coastal Wiki, NSEA) and should be checked against the 2013 paper when it is to hand.

Bottom Types

Defined in simulator/bottom.py. BSS at 20 deg grazing calibrates the amplitude; spectral params drive FFT macro-roughness.

name	mean_reflectivity	roughness_std (m)	bss_db_at_20deg	γ (spectral_exponent)	b (spectral_strength)	rms_height_m	default_rock_density
`sand`	0.30	0.002	-28.0	3.25	6.5e-5	0.02	0.0
`rippled_sand`	0.30	0.002	-28.0	3.25	6.5e-5	0.037	0.0
`mud`	0.15	0.001	-32.0	3.0	1.74e-7	0.005	0.0
`rock`	0.60	0.010	-15.0	2.0	4.82e-4	0.15	15.0
`gravel`	0.45	0.005	-22.0	2.5	1.08e-4	0.08	3.0
`silt`	0.20	0.001	-38.0	3.0	6.26e-8	0.003	0.0
`shell_hash`	0.38	0.006	-24.0	2.8	1.99e-5	0.025	0.0
`patchy_sand`	0.30	0.002	-28.0	3.25	6.5e-5	0.02	0.0

Angular law (2026-09-29): sand, rippled_sand, shell_hash and patchy_sand scatter by perturbation theory (angular_law: pt, below); mud, silt, gravel and rock stay Lambert. patchy_sand is sand with the sub-resolution patch texture (patch_shape 3, patch_size_m 0.015), tuned to K α = 4 for the Fidelity reference instrument at 10,000/m^2. Custom types default to Lambert and no texture unless they set these fields.

All built-ins are spectrum_mode: absolute (since 2026-09-22), so rms_height_m is the legacy target, used only for custom types cloned from them. Expected statistics on a 100 m x 100 m scene at a 5 cm cell (range slope = rms of the np.gradient slope that feeds the scatterer normals):

type	rms height	range slope	source of b
`sand`, `rippled_sand`	12.7 cm	3.5 deg	Lyons et al. 2022 Table I (γ 3.25, measured 3.3 to 4.9 deg), b set for 3.5 deg
`shell_hash`	4.4 cm	3.5 deg	sand's slope at its own γ
`mud`	5 mm	0.25 deg	reproduces the historical rms_height_m
`silt`	3 mm	0.15 deg	same
`gravel`	8 cm	12.6 deg	same
`rock`	15 cm	46.6 deg	same

Roughness spectrum

W(k) = b k^-γ is the Jackson & Richardson (2007) 2-D spectrum used by Lyons, Olson & Hansen (2022, JASA 152, 1363, Eq. 9): k is the 2-D wavenumber magnitude in rad/m and ∫W d²k is the height variance, so b is in m^(4-γ). The rms height is Eq. 12 and the total rms slope Eq. 11, both integrals between an outer scale k_L (scene size) and an inner scale k_c (cell). In absolute mode nothing is rescaled, so a larger scene has more rms height (h ∝ L^((γ-2)/2), Eq. 13) and a finer cell more slope (s² ∝ k_c^(4-γ)).

FFT normalisation (seafloor_builder._absolute_scale): complex noise n with E|n|² = 2, amplitude A = sqrt(b) (2π)^(1-γ/2) sqrt(N) / dx · f^(-γ/2) (f in cycles/m, N = n_x n_y), z = Re(ifft2(n A)). Then Var z = Σ W(k_m) Δk_x Δk_y, the discrete Eq. 12. seafloor_builder.powerlaw_moments gives the expected rms height and the rms slope per component as np.gradient measures it (central difference transfer sin²(k_x dx)/dx², Nyquist box π/dx); simulator/tests/test_seafloor_spectrum.py checks realisations against it. Two conventions to keep apart: Lyons' Eq. 11 is the TOTAL slope (both components) with a sharp cut at k_c = 2π/dx, and one component through a central difference on the grid carries 5.4x less variance at γ = 3.25. Table I's b = 2e-5 read literally is therefore 1.9 deg of range slope here (4.5 deg in the paper's own convention at dx = 5 cm).

shell_hash (2026-09-18) is sand armoured with shell fragments. No APL-UW entry exists for it; Jackson and Richardson (2007, ch. 13) put shell-covered sand a few dB above clean sand and below gravel at 100 to 400 kHz, hence -24 dB between sand and gravel, with a grainier per-scatterer roughness on a bed as flat as sand. It is usually painted as patches (bottom.patches) rather than as a zone.

roughness_std is sub-wavelength per-scatterer z-jitter; spectral_strength (absolute mode) or rms_height_m (rms mode) drives the FFT macro-roughness from SeafloorBuilder. bss_db_at_20deg=None means uncalibrated; Scene will refuse to build unless you pass allow_uncalibrated_bottom=True.

Angular law: sin^n or perturbation theory (2026-09-29)

A material scatters with sin^n(grazing) in power, n = angular_exponent (Lambert, 2), unless it sets angular_law: pt, which the built-in sand-like types do. Then the law is the first-order perturbation cross section of a fluid sediment, as the APL-UW handbook gives it (simulator/scattering_laws.py; APL-UW TR 9407, 1994, Eqs. 50 to 52, references/apluw_tr9407.pdf; Jackson and Richardson 2007, ch. 13):

sigma(theta) ~ sin^4(theta) |Y(theta)|^2 (cos^2 theta + 1/400)^(-γ/2),
Y = [(a_ρ - 1)^2 cos^2 + a_ρ^2 - κ^2] / [a_ρ sin + P]^2,
P = sqrt(κ^2 - cos^2),  κ = (1 + i δ) / a_ν,

with density_ratio a_ρ, sound_speed_ratio a_ν, loss_parameter δ and the material's spectral_exponent γ. Defaults are fine sand from Lyons, Olson and Hansen (2022, Table I, NW Elba): 1.87, 1.12, 0.01, 3.25. For a power-law spectrum the shape does not depend on frequency or b, so the BSS at ref_grazing_deg still sets the level. Relative to Lambert it is 5 to 12 dB lower at 2 to 5 deg grazing (about sin^4), peaks at the critical angle (26.8 deg for sand) and dips after it. Not modelled: the Kirchhoff term near normal incidence (the law is held at its 75 deg value above 75 deg) and volume scattering (insignificant for sand, Lyons 2022 Fig. 4). The raytrace reads it as a negative code in the per-cell exponent grid plus a 1025-point table in sin(grazing) (Heightfield.angular_law); scenes that use only sin^n are bit-identical to before. Measured on the Fidelity fixture: on a flat bed the image level against range follows the law to 0.17 dB, and on rough sand SI at long range rises far above the Lambert value (about 6 against 3.5 at 95 m for a 4 deg rms slope), as Lyons, Olson and Hansen report for real sand. The FW backend ignores angular_law (Lambert for every material). The validation fixtures pin their sand to Lambert (sand_lambert custom type), because their predictions are written for that law.

Seafloor surface sampled at each scatterer (2026-09-29)

Each seafloor scatterer takes the bilinearly interpolated height of the heightfield at its own (x, y), and the raytrace a bilinearly interpolated normal. Before, both were the value of the scatterer's 5 cm cell, which built every slope as 5 cm flat steps; on steep relief the image carried a texture at the cell spacing (on 3 cm rms ripples a spectral peak 500x the background at 30 m) and the far-range statistics were off by 3 to 6%. PS_SURFACE_INTERP=0 restores the nearest-cell surface, to reproduce a render made before this date bit for bit. Flat beds are unchanged.

Sub-resolution patch texture (2026-09-29)

With only independent point scatterers the speckle is a Poisson field: SI = 1 + 2/(ρ A_eff), so heavy tails come only from low density. A material with patch_shape > 0 multiplies each scatterer's power by a gamma variate of mean 1 and variance 1/patch_shape, constant over patches of mean side patch_size_m (Voronoi cells of a jittered lattice), the patch model of Abraham and Lyons (2002) (simulator/patch_texture.py). Patches much larger than the resolution cell give K speckle with α = patch_shape; patches smaller than it give α of about patch_shape times the patches per cell, so α then grows as resolution coarsens. The exact SI for a measured PSF is patch_texture.predicted_si. Mean backscatter is unchanged, and the texture comes from a hash of position, not from the scatterer RNG, so a scene without it is bit-identical. Across zone edges the texture variance is feathered like the other per-zone quantities. The built-in patchy_sand uses it; the sample-image dataset generator gives three quarters of its scenes' sand a texture with patch_shape log-uniform in 1.2 to 8, which spans the α of 2.3 to 7.9 measured at sea. Not in the FW backend.

APL-UW TR 9407 sediments (2026-10-01)

The handbook's 23 named bottoms (Table 2, Rough Rock to Clay) are built-in types named tr9407_<key> (tr9407_medium_sand, tr9407_clay, ...; keys in simulator/tr9407_bottom.SEDIMENTS). The sediments with a grain size are regenerated from the handbook's grain-size equations (Eqs. 2 to 10) and checked against the printed table; the three rock rows are copied. (MASTODON's copy of the table has Cobble ν = 2.50; the handbook prints 1.80, which is used here.) Each type takes from the handbook:

level: the backscattering strength at 20 deg, at 100 kHz, the top of the handbook's 10 to 100 kHz range. The types carry one level and do not extrapolate to SAS frequencies.
angular law: the same curve, tabulated (angular_law: tr9407, tr9407_sediment: <key>), so muds and clays follow their volume-scattering curve rather than Lambert. FW stays Lambert.
roughness: the handbook's γ (3.25) and w₂ as b = w₂ 10^(2γ-8), in absolute mode (the same 2-D spectrum convention). Rough Rock's w₂ = 0.207 cm⁴ gives about 1 m rms relief on a 100 m scene.

The curve is the handbook's full model. The interface term blends the Kirchhoff term (near normal incidence) with the composite-roughness term (perturbation theory averaged over the large-scale slope, with shadowing). Once the large-scale rms slope passes about 7 deg it switches to the empirical large-roughness term for rock and gravel. The volume term is added on top. It reproduces the handbook's Table 3 to its 0.1 dB rounding: 210 printed values for rough rock, rock, cobble, sandy gravel, and coarse and medium sand, at 1 to 90 deg and 10, 30 and 100 kHz (simulator/tests/test_tr9407_bottom.py). Two symbols are illegible in the handbook scan, and that agreement fixes both:

the Kirchhoff constant's exponent is 1/2 + 1/(2α), as printed in Mourad and Jackson (1989, Eq. 22, references/mourad_jackson1989.pdf). It also agrees with the exact Kirchhoff integral of Jackson et al. (1986, Eq. 38, references/jackson1986.pdf). The other reading puts the 90 deg peak 3.6 dB low.
the large-roughness term (Eq. 53) starts σ₁ sinᵐ[θπ / (1 + 0.81 θc²/θ²)], with θπ = 180 deg. The sine's argument rises to 90 deg well above the critical angle and collapses at low grazing.

Levels at 20 deg (100 kHz) run from -6.2 dB (Rough Rock) through -12.0 (Rock), -24.2 (Medium Sand) and -27.2 (Fine Sand, next to sand's -28) to -32.0 (Clay).

Radar (kind='rf')

The point-scatter backend runs radar scenes with the same kernels: Thorp is gated off, stop-and-hop is the default, the antenna can be depressed (depression_deg) and materials carry an angular exponent (angular_exponent, ref_grazing_deg) read per cell by the raytracer. Reference systems, the calibration convention and the validation ladder: sar_radar.md.

Real bathymetry grids (`kind: grid`)

Besides the procedural relief operators, a scene can fly over a real survey grid. tools/bathy_prepare.py reads a north-up GeoTIFF (float elevations, positive up, square pixels, a projected CRS such as UTM) and writes a prepared .npz:

conda run -n py312 python tools/bathy_prepare.py --src grid.tif --out grid_lp3.npz --lowpass 3

The scene then names it in its bathymetry list:

bottom:
  bathymetry:
    - {kind: grid, path: grid_lp3.npz, origin_e: 629591.8, origin_n: 3921379.7,
       heading_deg: 296.6, side: starboard}

(origin_e, origin_n) is the map point under the track at scene x = 0, heading_deg the along-track direction clockwise from north. Behavior (simulator/terrain_grid.py, run._apply_bathymetry):

Bottom following. Heights are taken relative to the bed under the track, smoothed over 20 m along-track, so the vehicle keeps its configured altitude above the floor it is flying.
Not clipped. The grid's relief is added as surveyed, outside bathymetry_clip_m; procedural operators in the same list still compose, scale and clip as before and add on top.
Starboard only. side: port is refused for now.
No pitch. The vehicle flies level, so the bed under the track should climb no more than about 3 degrees along-track; tools/bathy_plan_lines.py refuses lines that do, and lines that touch missing data.
Output in the grid's projection. With drc_format: geotiff, the DRC is written in the grid's EPSG with a rotated transform, so swaths on any heading land on the grid in a GIS.
Smooth the grid first. Survey grids carry processing artifacts (a fine cross-hatch, facet creases between cells) a few centimeters high that render as stripes of shadow at sonar grazing angles. A Gaussian low-pass of 3 m standard deviation (--lowpass 3) removes them, and with them most relief shorter than about 10 m; the material's own roughness supplies the centimeter scale.

tools/bathy_plan_lines.py tiles a survey line into consecutive 80 m swaths, one scene each.

Platform Motion

PlatformMotion.straight_line(speed, altitude, n_pings, ping_rate, heading=0, ...) generates a straight track. Sinusoidal perturbation pairs (sway_amplitude+sway_period, heave_*, roll_*, pitch_*, yaw_*) inject 6-DOF wiggle, all zero by default. Angles are radians. motion.altitude returns -pos_z; motion.time is arange(n_pings)/ping_rate.

motion.kind: spline_track builds a curved track instead, through simulator/track.py: a centripetal Catmull-Rom curve through track.control_points, a per-point speed law (piecewise-linear in arc length), and a kinematic attitude derivation with a first-order vehicle response per axis (track.dynamics). The FW backend requires a straight aperture, so its control points must be colinear at one altitude (rule R12d, appcore/validation/rules.py); the PS backend has no such restriction.

Noise

simulation.noise_level adds white complex Gaussian receiver noise of that RMS, relative to the echo scale (the dataset generators use it).

simulation.ambient_noise: true (2026-10-01) adds ambient sea noise from the APL-UW TR 9407 handbook (Section II.D; simulator/ambient_noise.py) on top of it:

surface noise from wind (ambient_noise_wind_mps, 10 m wind; air/sea temperature difference ambient_noise_air_sea_dt_c) and rain (ambient_noise_rain_mmph, capped at 10 mm/h as the handbook advises), radiated as a dipole and received through each element's own beam pattern, direct path, with Thorp absorption over the platform depth (ambient_noise_platform_depth_m; 0 = none). For an omnidirectional hydrophone this is the handbook's 2π A E₃(0.23 α D);
thermal noise, -15 + 20 log f(kHz) dB, less the element's directivity index and efficiency (ambient_noise_thermal_efficiency_db).

The noise is shaped to that spectrum across the band of the complex samples and is independent between channels. Echoes are simulated relative to a unit source at 1 m, so the projector's source level ambient_noise_source_level_db (default 210 dB re 1 µPa at 1 m) sets the SNR, which then follows the sonar equation. The handbook gives wind and rain for 1 to 100 kHz, so above that they are power-law extrapolations. At SAS frequencies thermal noise dominates below about 13 m/s of wind. The element pattern is baffled (no response behind the array), which leaves the echoes unchanged because the seafloor is only on the sonar side. Not modelled: multipath noise, near-surface bubble loss, shipping, snapping shrimp. Ignored for kind: rf.

SonarConfig Defaults

Match POSSM scene066: fc=300 kHz, fs=75 kHz, BW=60 kHz, pulse_length=4 ms, sos=1500, 36 channels at 33.33 mm spacing, 3 cm square TX/RX elements. create_baseband_chirp() returns the LFM waveform used for matched filtering; beamwidth is computed on demand via tx_beamwidth_rad / rx_beamwidth_rad (never hardcode).

Point-Scatter Performance (GPU path)

Beam-wedge tiling of a streamed cloud (2026-09-13)

When the cloud is too large to stay GPU-resident, the prefilter used to stream the whole cloud through PCIe on every ping (at 80k/m^2 on a 71 x 100 m scene: 580 M scatterers, 11.6 GB per ping, about 3 s of each 4.3 s ping, GPU at 30 to 40 %). simulator/cloud_window.py now sorts the cloud once into 2 m (x, y) tiles (stable radix sort on int16 keys, RAM-gated) and, per ping, streams and scans only the tiles the TX beam wedge can reach: a conservative per-tile test derived from the prefilter's own keep rule, with the tile's circumscribed radius and the cloud's z extent folded in, so the survivor set is unchanged (property tests over random attitudes, both sonar sides and elevated scatterers; byte-identical HDF5 and DRC output on a straight scene and on the arcing testing_fw_linear_track scene at 5000/m^2). Because the wedge is a thin fan around broadside, the selected share is similar at every heading: 19 % of the cloud per ping on the straight scene (6.3 to 15.4 pings/s, streaming forced), 33 % on the 90-degree arc (3.75 to 6.13 pings/s). The log reports the tiling time, ping 1's selection and the run's average streamed share. PS_X_WINDOW=0 disables it; the sort is skipped with a log line when free RAM cannot cover its transient (22 bytes per scatterer).

Resident clouds tile on the GPU (2026-09-13, later that day). A cloud that will stay device-resident (<= 25 % of free VRAM) is uploaded as is and sorted there by cloud_window.tile_cloud_gpu: same tile grid, same stable permutation as the host path (pinned by test_gpu_tiling_reproduces_the_host_permutation_exactly), so the HDF5 is unchanged, and there is no second host copy of the cloud. The host copy had pushed the build-to-sim peak RSS 0.9 GB over the preflight budget at 96 M scatterers (test_build_and_sim_handoff_peak_rss_budget). Streamed clouds still tile on the host, page-locked, because the sorted copy is what the per-ping copies read from.

Parallel tile sort (2026-09-13, evening). The host tiling is a counting sort run block-wise on the shared thread pool (simulator/cpu_pool.py, PS_BUILD_THREADS, 1 = serial): each block of 4 M rows computes its tile keys and a per-tile count, an exclusive prefix sum over blocks gives every block the destination of its rows for every tile (tile start + that tile's rows in earlier blocks), and each block stable-sorts its own keys and scatters its rows straight into the output. Blocks are contiguous index ranges, so the result is the same permutation as the global np.argsort(keys, kind="stable") it replaces (test_parallel_tile_sort_matches_the_global_stable_sort, serial and threaded, odd block sizes). The page-locked output is allocated on the pool first so the locking overlaps the key passes. 100 M rows with page-locked output (WSL): 21.6 s serial -> 5.0 s threaded, i.e. about 104 s -> 30 s for the 580 M cloud.

Single-pass, double-buffered streaming (2026-09-13, same day)

The streamed prefilter used to count survivors in one pass and fill them in a second, so every selected chunk crossed PCIe twice per ping, and each copy was synchronous. It is now one pass: two staging slots on the GPU, host-to-device copies issued on a dedicated non-blocking stream so chunk k+1 lands while chunk k's keep mask and compaction run, and survivors appended to persistent output buffers that grow on demand (no per-ping allocation, no concatenate). The tiled cloud is allocated page-locked when it is at most 40 % of free RAM (PS_PIN_CLOUD=0 keeps it pageable), so .set() streams straight from it with no host memcpy; a pageable cloud still bounces through two chunk-sized pinned buffers. Output order is chunk order, identical to before: the straight scene's HDF5 and DRC are bit-identical to the original two-pass unwindowed run. Streaming forced at 5000/m^2: straight 15.4 to 25.6 pings/s, the 90-degree arc 6.1 to 13.0 pings/s. PS_TIMING=1 prints the prefilter and whole-ping milliseconds with every progress line; on the arc scene the prefilter is now about 40 ms of an 85 ms ping.

Reference: the File -> New Scene default (60 m x 80 m image, HISAS-class 36-channel array, 150 pings) at 10,000 scatterers/m^2 = 57.6M scatterers, RTX 4500 Ada. Sim stage 226 s -> 19 s (per ping 1.56 s -> 0.11 s) on 2026-09-02; the beamformed image correlates 0.9998 with the old path (differences are fp32 summation-order speckle, ~53 dB down).

Per-ping stages, in the order they run, with the switch that restores the previous implementation. All are ON by default.

Stage	What it does now	Off switch
Cloud residency	The scatterer cloud is uploaded once and stays on the GPU when it needs <= 25 % of free VRAM; larger clouds stream from host chunk-wise as before (`simulate._resident_or_pinned_cloud`).	`PS_RESIDENT_CLOUD=0`
TX-beam prefilter	One mask kernel + one compaction over the resident cloud (`gpu_engine._prefilter_resident`) instead of a two-pass CuPy chain.	`PS_PREFILTER_KERNEL=0`
Line of sight	`check_los_hard_mip_kernel`: the ray march skips any segment (<= 8 cells long) whose highest point clears the 3x3-dilated max elevation of its coarse cell; the fine samples that do get tested are exactly the legacy ones, so visibility is bit-identical.	`PS_LOS_MIP=0`
Raytrace	Fused CUDA kernel (`fused_raytrace_kernel`, one thread per channel x scatterer) writes delays + complex64 amplitudes directly; replaces the CuPy elementwise chain (~10 (n_ch, n_sc, 3) temporaries per chunk).	`PS_FUSED_RAYTRACE=0`
Cull + sort	`raytrace_ping_fused_culled`: a per-scatterer score kernel yields max-channel amplitude and mean delay; the cull mask and delay argsort run on those two vectors; survivors' inputs are compacted and only then does the fused kernel write the (n_ch, n_keep) arrays, already sorted. Same keep rule as `_compact_active_scatterers`.	`PS_FUSED_CULL=0`
Echo (FFT path)	`accumulate_deltas_c64_kernel` reads the complex64 amplitudes directly and uses `sincospi` on the fractional carrier cycles (186 dB vs the split real/imag + `cos`/`sin` kernel).	`PS_DELTA_C64=0`
Raytrace + echo in one pass	`raytrace_echo_kernel`: one thread per delay-sorted survivor (indexed through the cull's index list, no compaction), looping over channels, scattering each pair's delta straight from registers into the FFT delta planes. Runs of equal bin address within a warp are summed by a segmented shuffle scan before ONE fp64 atomicAdd pair. The carrier phase is float32: fc·τ split with an FMA two-product so the fractional cycle is exact to ~6e-8 (fp64 `sincospi` on the raw ~3e4-cycle argument ran at 1/64 rate and was 3/4 of the kernel). Chunks are sized for this path's ~64 B/scatterer working set, so the default scene is one chunk per ping and the cull runs against the true per-ping peak. 146 dB vs the two-kernel path at equal chunking.	`PS_FUSED_ECHO=0`

None of this depends on the amplitude cull being enabled: with echo_amplitude_cull = 0 (Preferences -> Run) the fused pass keeps every visible scatterer and skips the score kernel. Only the TX-beam prefilter stays tied to the cull, since it is itself a cull. On the 4125i scene at 20k/m^2 (145M scatterers, 103 pings) the sim stage is 90 s at cull 0, 44 s at 1e-4, 36 s at 1e-3; before the decoupling cull 0 took 5 min. The default is 1e-4 since 2026-09-03 (simulate.run_simulation, the GUI preference run/echo_amplitude_cull, and therefore browser/AWS runs): 1e-3 culled the dim far-range ripple-trough and shadow-edge scatterers, so ripple fields beyond ~75 m broke into along-track dashes and point targets grew sidelobe crosses; 1e-4 is visually identical to culling off.

The ping loop also trims CuPy's memory pool whenever its cached total passes 60 % of the free VRAM measured at sim start (PS_POOL_TRIM_FRACTION). The fused path's per-ping transients change size with the survivor count, so cached blocks stop matching and the pool only grows; on the 145M-scatterer 4125i scene it reached 46 GB "total" on a 24 GB card. WSL's driver spills that to host memory rather than failing, and mid-track pings ran 2-10x slower (worse from the GUI, which holds VRAM of its own). With the trim: 67 s -> 36 s for that scene.

Second round, same scene at 20,000/m^2 (115M scatterers, 2.3 GB resident): per ping 215 ms -> 98 ms, sim 34.6 s -> 18.7 s, wall 71 s -> 56 s. What remains per ping (~98 ms): echo/raytrace kernel 30 ms, LOS 21 ms, score 14 ms, prefilter 12 ms, argsort 6 ms, the rest in compaction/finalize. LOS coarse factor is 16 (measured fastest, bit-identical at every factor). Outside the sim, instantiate() is ~10 s (the single-zone scatterer draw now writes through slices instead of a 57.6M-element fancy index; the zone-feather gaussian is 2.3 s of it) and the two PNG products ~7 s (zlib level 1; speckle does not compress and the row filter dominates).

Profiling recipe: simulator/perf/bench_echo.py (CUDA-event timers on a synthetic scene) or wrap gpu_engine / gpu_raytrace entry points with event timers around generate_ping_data_gpu on a real YAML.

Object placement hands arrays to the store (2026-09-13)

Every place_* method used to feed ObjectScattererStore.append one (x, y, z_ned, refl) tuple per scatterer. On the 9obj scene at 80k/m^2 that loop ran 39 M times and was 77 s of a 641 s run (the third-largest stage after echo simulation and the tiling sort). The placers now call ObjectScattererStore.extend_arrays(xs, ys, z_ned, refl) with the arrays they already had (z_ned may be a scalar); the store keeps the same float64 / complex128 chunks and the same order, so the object cloud is bit-identical to the tuple path (checked for every placer, including the Manta top disc whose per-point cos/sin became one vectorised call). The legacy append face stays for the texgraph placement drain and tests. Parsed mesh geometry is also cached per (path, mtime, size) as read-only arrays, so a scene that places the same OBJ many times parses it once. Measured on the 9obj scene (WSL): placement 77 s -> 20 s; what remains is trimesh's surface sampler and OBJ parsing of distinct files.

Mesh objects prepared in parallel (2026-09-13)

place_mesh_object is now two steps. prepare_mesh_object is pure: it loads the mesh, transforms it, rasterises its top surface over the footprint window and draws the surface scatterers from the object's own seed, reading only grid geometry and density. apply_mesh_prep writes the heightfield stamp (or feather), the protect and footprint masks, the facet record and the scatterers. instantiate() prepares every mesh object at once with prepare_mesh_objects on the shared thread pool (PS_BUILD_THREADS, 1 = inline), then applies them in YAML order between the other kinds, so overlapping stamps (max against the current bed) and per-object scatterer ranges are exactly what the one-by-one placer produced (test_prepared_mesh_placement_matches_ sequential_placement; the full 9obj placement pass was also compared against the committed code: scatterers, elevation, footprint, facet records and ranges identical). Distinct mesh files are parsed once, serially, before the threads start: trimesh's OBJ parse holds the GIL, and letting 16 threads contend for it made a low-density scene slower than serial. 9obj scene (WSL): placement 18.7 s serial -> 11.1 s at 80k/m^2, 4.3 -> 3.7 s at 500/m^2. What remains is mostly that serial parse (about 9 s for the nine distinct OBJ files); a process pool for the parse alone is the obvious next step if it matters.

Texture-graph nodes evaluated in parallel (2026-09-13)

evaluate_graph used to walk the topological order one node at a time. Nodes are pure functions of the context, their connected inputs and a seed derived from (base_seed, node_id), so independent nodes now run concurrently: a node is submitted to the shared thread pool (PS_BUILD_THREADS, 1 = the serial loop) as soon as every node it links from has finished, results are stored on the calling thread, and a producer is released only once every consumer has finished. The material palette keeps its last-writer-in-topological-order rule. Same node code on the same inputs, so the height, reflectivity, gains, material ids, palette and masks are byte-identical to the serial order (simulator/tests/texgraph/test_parallel_eval.py, plus the real 9obj graph applied to its heightfield compared against the committed evaluator). The 9obj graph is three generators (tang_ripple 4.6 s, musgrave 2.7 s, voronoi 2.6 s in WSL) and cheap operators, so the parallel time is the slowest generator: 10.5 s -> 7.1 s (WSL). The ripple node is a serial 10-iteration Hilbert/FFT loop; multi-threaded FFTs inside it are the next step if that stage matters.

Threaded scatterer sampling (2026-09-13)

Heightfield.generate_scatterers draws one PCG64 stream in a fixed order (uniform x, uniform y, then per zone normal / rayleigh / uniform phase), and cached references pin the cloud bit for bit to that stream, so the draws cannot be parallelised. They are only about a third of the sampler though: at 580 M scatterers roughly 30 s of draws against 55 s of cos/sin, cell indexing, grid gathers and float32 casts. Those now run on a bounded thread pool (_BlockRunner, PS_BUILD_THREADS, default min(16, cores), 1 = the serial reference path) while the calling thread keeps drawing the next block; numpy releases the GIL for all of it. Every draw still happens on the calling thread in the same order with the same block size, and each block gets the same operations it always did, so the output is byte-identical to the serial and to the historical unblocked path (test_threaded_build_is_bit_identical compares all three on a multi-zone scene with both texture gains and on a single-zone scene). The x and y draws are now blocked too, which removes the last full-length float64 transient (4.6 GB at 580 M). Measured on the 9obj scene at 40k/m^2 (290 M, WSL): scene.build 41.0 s serial -> 21.3 s threaded.

Calibration Harness

calibrate_simulator.py validates simulator fidelity against measurement targets. It builds three scenarios: level (flat seafloor, checks speckle CV), PSF (single calibration target, checks point-spread function), and ghosts (targets at multiple ranges, checks ghost suppression): runs the simulator + beamformer, and writes sim_output/calibration/calibration_report.txt with CV, dynamic range, contrast, and sidelobe metrics. Design spec: docs/superpowers/specs/2026-04-15-simulator-calibration-design.md.

Quick Presets

simulator/presets.py ships ready-to-run (scene, motion, cfg) tuples: flat_seafloor(), point_targets(), calibration_target(range_m=50), rippled_seafloor(), complex_scene(). All call scene.build() before returning.

← Signal processing Scene framework →