1e29264cymsimweb: WebGPU companion of cymsim with YAML/form editor, live spectrogram and playbackJeremy Magland 1// Rayleigh-integral listener with the rim-diffraction back term, as in cymsim/radiation.py, but
2// with the retarded-time sum evaluated at the audio rate: the per-element acceleration is
3// decimated to the audio rate by a windowed-sinc low-pass spanning DECIMATION_SPAN audio samples
4// (applied in the time loop through per-element accumulator slots) and delayed by an integer
5// plus a fractional number of audio samples with linear interpolation. The per-element
6// coefficients, the delays and the decimation taps are computed here.
7import type { Listener } from "../config/schema.ts";
8import type { CymbalProfile } from "./geometry.ts";
9import type { Mesh } from "./mesh.ts";
11/** Audio samples spanned by the decimation kernel (accumulator slots per element). */
12export const DECIMATION_SPAN = 8;
14/**
15 * Decimation taps: taps[w * m + kk] is the weight with which the sub-step of phase kk (0..m-1)
16 * within audio sample j0 contributes to audio sample j0 + w - (W/2 - 1), for slot w = 0..W-1.
17 * The kernel is a Blackman-windowed sinc with cutoff 0.45 fs, normalised to unit DC gain.
18 */
19export function decimationTaps(m: number, W = DECIMATION_SPAN): Float64Array {
20 const taps = new Float64Array(W * m), half = (W * m) / 2, fc = 0.45 / m;
21 let sum = 0;
22 for (let w = 0; w < W; w++) for (let kk = 0; kk < m; kk++) {
23 const d = kk + (W / 2 - 1 - w) * m; // sub-step offset from the sample centre
24 const x = Math.PI * fc * d, sinc = d === 0 ? 1 : Math.sin(2 * x) / (2 * x);
25 const u = (d + half) / (2 * half), win = 0.42 - 0.5 * Math.cos(2 * Math.PI * u) + 0.08 * Math.cos(4 * Math.PI * u);
26 taps[w * m + kk] = sinc * Math.max(win, 0);
27 sum += taps[w * m + kk];
28 }
29 for (let i = 0; i < taps.length; i++) taps[i] /= sum;
30 return taps;
31}
33export interface Radiation {
34 coefF: Float64Array; coefB: Float64Array; // rho A / (2 pi r), front and back
35 i0F: Uint32Array; frF: Float64Array; // delay in audio samples: integer and fraction
36 i0B: Uint32Array; frB: Float64Array;
37 ringLength: number; // history ring length in audio samples
38 taps: Float64Array; // (DECIMATION_SPAN, substeps) decimation weights
39 rFront: Float64Array; rBack: Float64Array;
40}
42export function buildRadiation(mesh: Mesh, profile: CymbalProfile, z0Node: Float64Array, listener: Listener, fs: number, substeps: number, nRim = 720): Radiation {
43 const E = mesh.e, [Lx, Ly, Lz] = listener.position;
44 const cx = new Float64Array(E), cy = new Float64Array(E), cz = new Float64Array(E);
45 for (let e = 0; e < E; e++) {
46 cx[e] = mesh.centroid[2 * e]; cy[e] = mesh.centroid[2 * e + 1];
47 cz[e] = (z0Node[mesh.tri[3 * e]] + z0Node[mesh.tri[3 * e + 1]] + z0Node[mesh.tri[3 * e + 2]]) / 3;
48 }
49 let R = 0;
50 for (let n = 0; n < mesh.n; n++) R = Math.max(R, Math.hypot(mesh.xy[2 * n], mesh.xy[2 * n + 1]));
51 const rimX = new Float64Array(nRim), rimY = new Float64Array(nRim), rimZ = new Float64Array(nRim), d2 = new Float64Array(nRim);
52 for (let j = 0; j < nRim; j++) {
53 const phi = (2 * Math.PI * j) / nRim;
54 rimX[j] = R * Math.cos(phi); rimY[j] = R * Math.sin(phi); rimZ[j] = profile.z0(rimX[j], rimY[j]);
55 d2[j] = Math.hypot(rimX[j] - Lx, rimY[j] - Ly, rimZ[j] - Lz);
56 }
57 const rFront = new Float64Array(E), rBack = new Float64Array(E);
58 const coefF = new Float64Array(E), coefB = new Float64Array(E), i0F = new Uint32Array(E), frF = new Float64Array(E);
59 const i0B = new Uint32Array(E), frB = new Float64Array(E);
60 const { c_air: c, rho_air: rho, back_diffraction_gain: g } = listener;
61 let maxI0 = 0;
62 for (let e = 0; e < E; e++) {
63 rFront[e] = Math.hypot(Lx - cx[e], Ly - cy[e], Lz - cz[e]);
64 let best = Infinity;
65 for (let j = 0; j < nRim; j++) best = Math.min(best, Math.hypot(cx[e] - rimX[j], cy[e] - rimY[j], cz[e] - rimZ[j]) + d2[j]);
66 rBack[e] = best;
67 coefF[e] = (rho * mesh.area[e]) / (2 * Math.PI * rFront[e]);
68 coefB[e] = (g * rho * mesh.area[e]) / (2 * Math.PI * rBack[e]);
69 const dF = (rFront[e] / c) * fs, dB = (rBack[e] / c) * fs;
70 i0F[e] = Math.floor(dF); frF[e] = dF - i0F[e]; i0B[e] = Math.floor(dB); frB[e] = dB - i0B[e];
71 maxI0 = Math.max(maxI0, i0F[e], i0B[e]);
72 }
73 return { coefF, coefB, i0F, frF, i0B, frB, ringLength: maxI0 + 4, taps: decimationTaps(substeps), rFront, rBack };
74}