982850fWoodblock strike simulator: Python reference model and interactive web appJeremy Magland 1// Radiation: retarded free-field pressure at the listener point.
2// Port of radiate() and check_setup() in woodblock.py. The top and bottom
3// faces move together, so each strip is a local dipole; the surface velocity
4// sum is computed once per strip and reused for both faces with opposite sign.
6import type { WoodblockParams } from '../params';
7import type { ModalBasisData } from './basis';
8import { elementShapes } from './basis';
10export interface RadiateResult {
11 p: Float64Array; // pascals at the listener, at fsHi
12 warnings: string[];
13}
15export function radiate(
16 cfg: WoodblockParams,
17 Qdd: Float64Array[],
18 Xdd: Float64Array,
19 basis: ModalBasisData,
20 fsHi: number,
21 M: number,
22 onProgress?: (frac: number) => void,
23): RadiateResult {
24 const rho0 = cfg.air_density;
25 const c = cfg.speed_of_sound;
26 const N = basis.n;
27 const T = Xdd.length;
28 const L = cfg.length_m;
29 const W = cfg.width_m;
30 const halfH = 0.5 * cfg.block_height_m;
31 const patch = (W * L) / M;
33 const { phiElem, sElem } = elementShapes(basis, M);
35 const lx = cfg.listener_x_m;
36 const ly = cfg.listener_y_m;
37 const lz = cfg.listener_z_m;
39 const distTop = new Float64Array(M);
40 const distBot = new Float64Array(M);
41 let maxDist = 0;
42 for (let j = 0; j < M; j++) {
43 const x = (sElem[j] - 0.5) * L;
44 const dx = x - lx;
45 const dy = 0.0 - ly;
46 const dzT = halfH - lz;
47 const dzB = -halfH - lz;
48 distTop[j] = Math.sqrt(dx * dx + dy * dy + dzT * dzT);
49 distBot[j] = Math.sqrt(dx * dx + dy * dy + dzB * dzB);
50 maxDist = Math.max(maxDist, distTop[j], distBot[j]);
51 }
52 const sx = cfg.slit_x_m - lx;
53 const sy = cfg.slit_y_m - ly;
54 const sz = cfg.slit_z_m - lz;
55 const slitDist = Math.sqrt(sx * sx + sy * sy + sz * sz);
56 maxDist = Math.max(maxDist, slitDist);
58 const maxDelay = Math.ceil((maxDist / c) * fsHi) + 2;
59 const out = new Float64Array(T + maxDelay);
60 const wdd = new Float64Array(T);
62 for (let j = 0; j < M; j++) {
63 // Surface acceleration of strip j, shared by the +z and -z sources
64 wdd.fill(0);
65 for (let n = 0; n < N; n++) {
66 const phi = phiElem[n][j];
67 if (phi === 0) continue;
68 const q = Qdd[n];
69 for (let t = 0; t < T; t++) wdd[t] += phi * q[t];
70 }
71 const ampTop = (patch * rho0) / (4.0 * Math.PI * distTop[j]);
72 const ampBot = -(patch * rho0) / (4.0 * Math.PI * distBot[j]);
73 const kTop = Math.round((distTop[j] / c) * fsHi);
74 const kBot = Math.round((distBot[j] / c) * fsHi);
75 for (let t = 0; t < T; t++) {
76 out[kTop + t] += ampTop * wdd[t];
77 out[kBot + t] += ampBot * wdd[t];
78 }
79 if (onProgress) onProgress((j + 1) / M);
80 }
82 // Slit monopole
83 const ampSlit = (cfg.slit_area_m2 * rho0) / (4.0 * Math.PI * slitDist);
84 const kSlit = Math.round((slitDist / c) * fsHi);
85 for (let t = 0; t < T; t++) out[kSlit + t] += ampSlit * Xdd[t];
87 // Setup warnings (check_setup in the reference)
88 const warnings: string[] = [];
89 const spacing = L / M;
90 const fMax = basis.f[basis.n - 1];
91 if (spacing > (0.25 * c) / fMax) {
92 const need = Math.ceil(L / ((0.25 * c) / fMax));
93 warnings.push(
94 `element spacing is coarse at ${fMax.toFixed(0)} Hz; consider ${need} radiating elements`,
95 );
96 }
97 const d = Math.sqrt(lx * lx + ly * ly + lz * lz);
98 if (d < 0.5 * L) {
99 warnings.push(
100 `listener is ${(d * 100).toFixed(0)} cm from the block centre, inside the source region; point-source elements are crude here`,
101 );
102 }
104 return { p: out.subarray(0, T).slice(), warnings };
105}