1/**
2 * The cymbal problem: a reciprocity solve for a listener near a rigid cymbal
3 * (problems/cymbal.py in the CUDA proof of concept). A point source at the
4 * listener's position L radiates in free space, and the only obstacle is a
5 * stationary, sound-hard cymbal:
6 *
7 * (Delta + k^2) u = -delta(x - L) outside the cymbal,
8 * du/dn = 0 on the cymbal,
9 *
10 * with the radiation condition. By reciprocity u(x) is also the pressure at the
11 * listener due to a unit source at x. We solve the exterior mode of the same
12 * impedance equation with lambda = 0, and show the field on three planes through
13 * a 1 m box around the cymbal. Lengths in meters, sound speed 343 m/s.
14 */
15import { CymbalShape, insideTest } from './cymbal.ts';
16import type { CVec, Slice } from './problem.ts';
18export const C_SOUND = 343;
19export const RADIUS = 0.2;
20export const LISTENER: [number, number, number] = [0.25, 0, 0.45];
21/** the evaluation box: [lo, hi] in x, y, z */
22export const BOX: [number, number][] = [
23 [-0.5, 0.5],
24 [-0.5, 0.5],
25 [-0.35, 0.65],
26];
27/** height of the horizontal slice, below the cymbal */
28export const SLICE_BELOW = -0.1;
30/**
31 * Point sources inside the cymbal whose field is the exact solution outside, for
32 * the self-test: positions (3, n) and strengths. A single source would be at most
33 * T/2 = 1 cm from the surface, and its data would vary on that scale, which the
34 * main problem's data does not. Instead the sources are spread 3 mm apart over
35 * the midsurface out to r = rMax, with strengths tapering smoothly to zero.
36 */
37export function selftestSources(shape: CymbalShape, rMax = 0.12, spacing = 0.003): { x0s: Float64Array; ws: Float64Array } {
38 const g: number[] = [];
39 for (let t = -rMax; t <= rMax + spacing / 2; t += spacing) g.push(t);
40 const pts: number[][] = [];
41 const ws: number[] = [];
42 for (const x of g) {
43 for (const y of g) {
44 const rho = x * x + y * y;
45 if (rho >= rMax * rMax) continue;
46 pts.push([x, y, shape.F(rho).v]);
47 ws.push(((1 - rho / rMax ** 2) ** 3 * spacing ** 2) / ((Math.PI * rMax ** 2) / 4));
48 }
49 }
50 const n = pts.length;
51 const x0s = new Float64Array(3 * n);
52 for (let j = 0; j < n; j++) for (let d = 0; d < 3; d++) x0s[d * n + j] = pts[j][d];
53 return { x0s, ws: Float64Array.from(ws) };
54}
56/** sum_j ws[j] G(x, x0s_j) at points x (3, n), or with nrm given its normal derivative. */
57export function sourceSum(x: Float64Array, x0s: Float64Array, ws: Float64Array, k: number, nrm?: Float64Array): CVec {
58 const n = x.length / 3;
59 const m = ws.length;
60 const out = new Float64Array(2 * n);
61 for (let i = 0; i < n; i++) {
62 let sr = 0;
63 let si = 0;
64 for (let j = 0; j < m; j++) {
65 const dx = x[i] - x0s[j];
66 const dy = x[n + i] - x0s[m + j];
67 const dz = x[2 * n + i] - x0s[2 * m + j];
68 const r = Math.hypot(dx, dy, dz);
69 const s = ws[j] / (4 * Math.PI * r);
70 let gr = s * Math.cos(k * r);
71 let gi = s * Math.sin(k * r);
72 if (nrm) {
73 const dn = (dx * nrm[i] + dy * nrm[n + i] + dz * nrm[2 * n + i]) / r;
74 // g (ik - 1/r) dn
75 const tr = (-gr / r - gi * k) * dn;
76 const ti = (gr * k - gi / r) * dn;
77 gr = tr;
78 gi = ti;
79 }
80 sr += gr;
81 si += gi;
82 }
83 out[2 * i] = sr;
84 out[2 * i + 1] = si;
85 }
86 return out;
87}
89/** n points around the cymbal, at least 5 mm from its surface (deterministic). */
90export function selftestPoints(shape: CymbalShape, n: number): Float64Array {
91 const near = insideTest(shape, 0.005);
92 const pts: number[][] = [];
93 let s = 12345;
94 const rnd = () => {
95 s = (Math.imul(s, 1664525) + 1013904223) >>> 0;
96 return s / 4294967296;
97 };
98 while (pts.length < n) {
99 const x = [-0.3 + 0.6 * rnd(), -0.3 + 0.6 * rnd(), -0.15 + 0.35 * rnd()];
100 if (!near(x)) pts.push(x);
101 }
102 const out = new Float64Array(3 * n);
103 for (let j = 0; j < n; j++) for (let d = 0; d < 3; d++) out[d * n + j] = pts[j][d];
104 return out;
105}
107/**
108 * Power radiated to infinity, k int |u_inf|^2 dOmega, for u = G(., x0) + S_k sigma +
109 * i alpha D_k phi1, with u ~ u_inf exp(ikr)/r. The far field pattern is computed
110 * with the smooth rule at the nodes, and the sphere integral with Gauss-Legendre
111 * in cos(theta) and the trapezoid rule in phi, exact for the band limit. For a
112 * rigid cymbal it equals the power emitted by the source, k/(4 pi) + Im v(x0).
113 */
114export function farFieldPower(
115 k: number,
116 x0: ArrayLike<number>,
117 xs: Float64Array,
118 nrm: Float64Array,
119 wts: Float64Array,
120 sigma: CVec,
121 phi1: CVec,
122 alpha: number,
123 leggauss: (n: number) => { x: Float64Array; w: Float64Array },
124): number {
125 const n = wts.length;
126 let rmax = Math.hypot(x0[0], x0[1], x0[2]);
127 for (let i = 0; i < n; i++) rmax = Math.max(rmax, Math.hypot(xs[i], xs[n + i], xs[2 * n + i]));
128 const L = Math.ceil(k * rmax) + 20;
129 const { x: ct, w: wct } = leggauss(L + 1);
130 const nphi = 2 * L + 2;
131 // qs = w sigma, qd = i alpha w phi1
132 const qs = new Float64Array(2 * n);
133 const qd = new Float64Array(2 * n);
134 for (let i = 0; i < n; i++) {
135 qs[2 * i] = wts[i] * sigma[2 * i];
136 qs[2 * i + 1] = wts[i] * sigma[2 * i + 1];
137 qd[2 * i] = -alpha * wts[i] * phi1[2 * i + 1];
138 qd[2 * i + 1] = alpha * wts[i] * phi1[2 * i];
139 }
140 let total = 0;
141 for (let a = 0; a <= L; a++) {
142 const st = Math.sqrt(1 - ct[a] ** 2);
143 for (let b = 0; b < nphi; b++) {
144 const ph = (2 * Math.PI * b) / nphi;
145 const d = [st * Math.cos(ph), st * Math.sin(ph), ct[a]];
146 const p0 = -k * (d[0] * x0[0] + d[1] * x0[1] + d[2] * x0[2]);
147 let ur = Math.cos(p0);
148 let ui = Math.sin(p0);
149 for (let i = 0; i < n; i++) {
150 const ph1 = -k * (d[0] * xs[i] + d[1] * xs[n + i] + d[2] * xs[2 * n + i]);
151 const er = Math.cos(ph1);
152 const ei = Math.sin(ph1);
153 const dn = d[0] * nrm[i] + d[1] * nrm[n + i] + d[2] * nrm[2 * n + i];
154 // e (qs + (-i k dn) qd)
155 const cr = qs[2 * i] + k * dn * qd[2 * i + 1];
156 const ci = qs[2 * i + 1] - k * dn * qd[2 * i];
157 ur += er * cr - ei * ci;
158 ui += er * ci + ei * cr;
159 }
160 total += ((wct[a] * 2 * Math.PI) / nphi) * (ur * ur + ui * ui);
161 }
162 }
163 return (k * total) / (16 * Math.PI * Math.PI);
164}
166/**
167 * Three planes through the box: xz through the listener, yz through the cymbal's
168 * axis, and xy below the cymbal. The mask is 1 outside the cymbal and away from
169 * the source.
170 */
171export function boxSlices(ngrid: number, x0: ArrayLike<number>, shape: CymbalShape): Slice[] {
172 const inside = insideTest(shape, 1e-3);
173 const defs: [Slice['name'], [number, number], number, number][] = [
174 ['xz', [0, 2], 1, x0[1]],
175 ['yz', [1, 2], 0, 0],
176 ['xy', [0, 1], 2, SLICE_BELOW],
177 ];
178 const out: Slice[] = [];
179 for (const [name, axes, fixed, value] of defs) {
180 const N = ngrid * ngrid;
181 const pts = new Float64Array(3 * N);
182 const mask = new Uint8Array(N);
183 const [ea, eb] = [BOX[axes[0]], BOX[axes[1]]];
184 const h = (ea[1] - ea[0]) / (ngrid - 1);
185 for (let row = 0; row < ngrid; row++) {
186 for (let col = 0; col < ngrid; col++) {
187 const q = row * ngrid + col;
188 const x = [0, 0, 0];
189 x[axes[0]] = ea[0] + ((ea[1] - ea[0]) * col) / (ngrid - 1);
190 x[axes[1]] = eb[0] + ((eb[1] - eb[0]) * row) / (ngrid - 1);
191 x[fixed] = value;
192 pts[q] = x[0];
193 pts[N + q] = x[1];
194 pts[2 * N + q] = x[2];
195 const atSource = Math.hypot(x[0] - x0[0], x[1] - x0[1], x[2] - x0[2]) < 1e-3 * h;
196 mask[q] = inside(x) || atSource ? 0 : 1;
197 }
198 }
199 out.push({ name, axes, fixed, value, extent: [[ea[0], ea[1]], [eb[0], eb[1]]], pts, mask });
200 }
201 return out;
202}