1/**
2 * Scenes: what the medium is.
3 *
4 * A scene returns two fields of position, the sound speed c (m/s) and the
5 * absorption sig (1/s). Both are ordinary fields, which is what lets one
6 * function describe both the scatterer and the open boundary: the absorbing
7 * layer around the outside is the statement that the medium swallows sound out
8 * there, and a scene is free to put absorption inside the domain too, making a
9 * lossy scatterer.
10 *
11 * These are built on the CPU, once, whenever a parameter changes. At 128^3
12 * that is two million evaluations and takes well under a second; at 192^3 it
13 * is seven million and is noticeable. The flat sibling writes its scenes as
14 * MATLAB and evaluates them through numbl's interpreter, which is the more
15 * interesting arrangement and the obvious thing to bring over here later.
16 *
17 * Interfaces are smoothed over about a cell. A jump between two neighbouring
18 * cells is not resolved by the grid: it scatters the grid's own staircase
19 * rather than the shape that was asked for.
20 */
21import { coord, type Grid } from './grid.ts';
22import { C_AIR, SPONGE_FRAC, SPONGE_MAX } from './units.ts';
24export interface SceneParam {
25 key: string;
26 label: string;
27 min: number;
28 max: number;
29 step: number;
30 value: number;
31 unit?: string;
32}
34export interface Medium {
35 c: Float32Array<ArrayBuffer>;
36 sig: Float32Array<ArrayBuffer>;
37 cmin: number;
38 cmax: number;
39}
41export interface Scene {
42 label: string;
43 blurb: string;
44 params: SceneParam[];
45 /** Source settings this scene wants, applied when it is selected. */
46 source?: Record<string, number>;
47 /** Sound speed at (x, y, z), in m/s, and any absorption of its own. */
48 medium(x: number, y: number, z: number, h: number, v: Record<string, number>): [number, number];
49}
51/**
52 * Absorption profile for an open boundary: zero in the interior, ramping up
53 * quadratically over a layer of width w inside each face of the cube and
54 * reaching SPONGE_MAX at the wall. This is what makes the finite grid stand in
55 * for an unbounded medium.
56 *
57 * The ramp is gradual on purpose. An absorbing layer is itself an impedance
58 * mismatch, so a sudden one reflects; spreading it over a couple of
59 * wavelengths keeps that small. It is not a perfectly matched layer, and at
60 * grazing incidence it does leak.
61 */
62function sponge(x: number, y: number, z: number, L: number): number {
63 const w = SPONGE_FRAC * L;
64 const dx = Math.max(0, w - (L / 2 - Math.abs(x))) / w;
65 const dy = Math.max(0, w - (L / 2 - Math.abs(y))) / w;
66 const dz = Math.max(0, w - (L / 2 - Math.abs(z))) / w;
67 const d = Math.max(dx, Math.max(dy, dz));
68 return SPONGE_MAX * d * d;
69}
71/** A smoothed indicator: 1 well inside the surface, 0 well outside. */
72const inside = (signedDistance: number, h: number): number =>
73 0.5 * (1 - Math.tanh(signedDistance / (1.5 * h)));
75/** Deterministic value noise on a coarse lattice, trilinearly interpolated. */
76function makeNoise(m: number, seed: number): (u: number, v: number, w: number) => number {
77 let s = (seed * 0x9e3779b9) >>> 0;
78 const rand = () => {
79 s = (s + 0x6d2b79f5) >>> 0;
80 let t = s;
81 t = Math.imul(t ^ (t >>> 15), t | 1);
82 t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
83 return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
84 };
85 const g = new Float64Array(m * m * m);
86 for (let i = 0; i < g.length; i++) g[i] = 2 * rand() - 1;
87 const at = (i: number, j: number, k: number) => {
88 const w = (a: number) => ((a % m) + m) % m;
89 return g[w(i) + m * w(j) + m * m * w(k)];
90 };
91 // u, v, w in [0, 1) over the whole domain.
92 return (u, v, w) => {
93 const fx = u * m;
94 const fy = v * m;
95 const fz = w * m;
96 const i = Math.floor(fx);
97 const j = Math.floor(fy);
98 const k = Math.floor(fz);
99 const sx = fx - i;
100 const sy = fy - j;
101 const sz = fz - k;
102 // Smoothstep on each axis so the field is continuous in its derivative,
103 // which keeps the medium resolvable by the grid.
104 const ex = sx * sx * (3 - 2 * sx);
105 const ey = sy * sy * (3 - 2 * sy);
106 const ez = sz * sz * (3 - 2 * sz);
107 const lerp = (a: number, b: number, t: number) => a + (b - a) * t;
108 const c00 = lerp(at(i, j, k), at(i + 1, j, k), ex);
109 const c10 = lerp(at(i, j + 1, k), at(i + 1, j + 1, k), ex);
110 const c01 = lerp(at(i, j, k + 1), at(i + 1, j, k + 1), ex);
111 const c11 = lerp(at(i, j + 1, k + 1), at(i + 1, j + 1, k + 1), ex);
112 return lerp(lerp(c00, c10, ey), lerp(c01, c11, ey), ez);
113 };
114}
116/** Rebuilt whenever the random medium's seed or scale changes. */
117let noiseCache: { key: string; f: (u: number, v: number, w: number) => number } | null = null;
119export const scenes: Record<string, Scene> = {
120 sphere: {
121 label: 'Sphere',
122 blurb:
123 'One spherical scatterer in a uniform background: the reference case, and the ' +
124 'three-dimensional problem with a classical series solution. Speeds above 1 behave ' +
125 'nearly rigid, below 1 nearly pressure-release, and exactly 1 is no scatterer at all.',
126 params: [
127 { key: 'R', label: 'radius', min: 0.05, max: 0.5, step: 0.01, value: 0.25, unit: 'm' },
128 { key: 'cin', label: 'speed ratio', min: 0.1, max: 4, step: 0.05, value: 2.5 },
129 { key: 'absorb', label: 'absorption', min: 0, max: 6000, step: 100, value: 0, unit: '1/s' },
130 ],
131 medium(x, y, z, h, v) {
132 const r = Math.sqrt(x * x + y * y + z * z);
133 const q = inside(r - v.R, h);
134 return [C_AIR * (1 + (v.cin - 1) * q), v.absorb * q];
135 },
136 },
138 pair: {
139 label: 'Two spheres',
140 blurb:
141 'Two identical scatterers side by side. Each one reradiates what it receives, ' +
142 'including what it receives from the other, so the pattern behind them is not the sum ' +
143 'of two single-sphere patterns.',
144 params: [
145 { key: 'R', label: 'radius', min: 0.05, max: 0.35, step: 0.01, value: 0.16, unit: 'm' },
146 { key: 'sep', label: 'separation', min: 0.15, max: 1.0, step: 0.01, value: 0.5, unit: 'm' },
147 { key: 'cin', label: 'speed ratio', min: 0.1, max: 4, step: 0.05, value: 3 },
148 ],
149 medium(x, y, z, h, v) {
150 const d = v.sep / 2;
151 const r1 = Math.sqrt(x * x + (y - d) * (y - d) + z * z);
152 const r2 = Math.sqrt(x * x + (y + d) * (y + d) + z * z);
153 const q = Math.max(inside(r1 - v.R, h), inside(r2 - v.R, h));
154 return [C_AIR * (1 + (v.cin - 1) * q), 0];
155 },
156 },
158 aperture: {
159 label: 'Aperture',
160 blurb:
161 'A screen with a circular hole in it. Behind the hole is the three-dimensional ' +
162 'diffraction pattern, which a plane cut through the flat problem cannot show: the ' +
163 'Airy-like rings are a property of the circle, not of the slit.',
164 params: [
165 { key: 'a', label: 'hole radius', min: 0.04, max: 0.5, step: 0.01, value: 0.15, unit: 'm' },
166 { key: 'th', label: 'screen thickness', min: 0.02, max: 0.2, step: 0.01, value: 0.06, unit: 'm' },
167 { key: 'cw', label: 'screen speed', min: 0.05, max: 1, step: 0.05, value: 0.2 },
168 ],
169 medium(x, y, z, h, v) {
170 // The screen is slow rather than fast. Reflection at an interface goes
171 // as |c2 - c1|/(c2 + c1) at constant density, so c = 0.2 reflects about
172 // as much as c = 5 would, but the timestep is set by the fastest speed
173 // anywhere on the grid, so a slow screen is free and a fast one taxes
174 // every step of the whole run. What a slow screen costs instead is
175 // resolution inside itself, which its own absorption then swallows.
176 const slab = inside(Math.abs(x) - v.th / 2, h);
177 const rho = Math.sqrt(y * y + z * z);
178 const hole = inside(rho - v.a, h);
179 const q = slab * (1 - hole);
180 return [C_AIR * (1 + (v.cw - 1) * q), 4000 * q];
181 },
182 },
184 random: {
185 label: 'Random medium',
186 blurb:
187 'Weak random structure everywhere, with no scatterer in particular. A pulse through ' +
188 'it arrives on time and then keeps arriving: multiple scattering turns the tail into a ' +
189 'coda, which is what a real inhomogeneous medium does.',
190 params: [
191 { key: 'amp', label: 'contrast', min: 0, max: 0.6, step: 0.01, value: 0.25 },
192 { key: 'scale', label: 'blob size', min: 0.05, max: 0.6, step: 0.01, value: 0.18, unit: 'm' },
193 { key: 'seed', label: 'seed', min: 1, max: 20, step: 1, value: 1 },
194 ],
195 medium(x, y, z, _h, v) {
196 // `buildMedium` has rebuilt the lattice for this scale and seed before
197 // calling us, so it is always here.
198 const f = noiseCache!.f;
199 const L = DOM.L;
200 return [C_AIR * (1 + v.amp * f((x + L / 2) / L, (y + L / 2) / L, (z + L / 2) / L)), 0];
201 },
202 },
203};
205/** The domain the scenes are currently being evaluated over. Set by `build`;
206 * the random medium needs it to place its lattice. */
207const DOM = { L: 2 };
209export const sceneKeys = Object.keys(scenes);
211/** Evaluate a scene over the grid, adding the absorbing layer. */
212export function buildMedium(scene: Scene, grid: Grid, v: Record<string, number>): Medium {
213 DOM.L = grid.L;
214 if (scene === scenes.random) {
215 const m = Math.max(2, Math.round(grid.L / Math.max(v.scale, 1e-3)));
216 const key = `${m}:${v.seed}`;
217 if (!noiseCache || noiseCache.key !== key) noiseCache = { key, f: makeNoise(m, v.seed) };
218 }
220 const { n, npts, L, h } = grid;
221 const c = new Float32Array(npts);
222 const sig = new Float32Array(npts);
223 let cmin = Infinity;
224 let cmax = -Infinity;
225 for (let iz = 0; iz < n; iz++) {
226 const z = coord(iz, grid);
227 for (let iy = 0; iy < n; iy++) {
228 const y = coord(iy, grid);
229 const row = n * iy + n * n * iz;
230 for (let ix = 0; ix < n; ix++) {
231 const x = coord(ix, grid);
232 const [cv, sv] = scene.medium(x, y, z, h, v);
233 const k = ix + row;
234 c[k] = cv;
235 sig[k] = sv + sponge(x, y, z, L);
236 if (cv < cmin) cmin = cv;
237 if (cv > cmax) cmax = cv;
238 }
239 }
240 }
241 return { c, sig, cmin, cmax };
242}