1% The dulcimer's body: a rigid box under the string, with a round sound hole
2% in its top plate.
3%
4% The top plate lies at z = 0. The box hangs below it (down to z = -boxd) and
5% the string runs along x just above it, at height `gap`. Everything here is
6% metres and inverse seconds.
7%
8% The shell is a *mask*, not a material: `wall` is 0 inside the solid and 1
9% in air, and the solver's masked Laplacian drops any flux into a masked
10% cell, which is the rigid (Neumann) condition. So the walls reflect
11% everything and cost nothing — no timestep penalty, unlike a fast material.
12% What a rigid mask cannot do is absorb, so `absorb` paints a thin layer of
13% ordinary volume absorption over the shell's surfaces (the smoothed mask's
14% transition skin, 4*wall*(1-wall), which peaks exactly at the surface); that
15% is what lets the cavity's ring be tamed like real wood tames it.
16%
17% The masks are drawn with a sharp profile (about half a cell) on purpose. A
18% box is grid-aligned, so sharpness costs no staircase artefacts, and a plate
19% only two cells thick has to reach mask ≈ 0 at its centre or it leaks.
20function [c, sig, wall, lineprof, boardprof] = medium(x, y, z, h, Lx, Ly, Lz, c0, Ls, body, boxl, boxw, boxd, thick, holer, holex, gap, patchr, absorb)
21 s = 0.45 * h;
23 % The solid: the outer box minus its interior cavity.
24 ox = 0.5 * (1 - tanh((abs(x) - boxl/2) / s));
25 oy = 0.5 * (1 - tanh((abs(y) - boxw/2) / s));
26 oz = 0.5 * (1 - tanh((abs(z + boxd/2) - boxd/2) / s));
27 ix = 0.5 * (1 - tanh((abs(x) - (boxl/2 - thick)) / s));
28 iy = 0.5 * (1 - tanh((abs(y) - (boxw/2 - thick)) / s));
29 iz = 0.5 * (1 - tanh((abs(z + boxd/2) - (boxd/2 - thick)) / s));
30 shell = max(0, ox .* oy .* oz - ix .* iy .* iz);
32 % The sound hole: a cylinder of radius holer at (holex, 0), cut out of the
33 % top plate only.
34 r = sqrt((x - holex).^2 + y.^2);
35 hole = 0.5 * (1 - tanh((r - holer) / s));
36 top = 0.5 * (1 - tanh((abs(z + thick/2) - thick/2) / s));
37 shell = max(0, shell - hole .* top);
39 % `body` at 0 takes the box away entirely: a bare string in open air. The
40 % bridge patch stays where the plate would have been, so the bridge route
41 % keeps working as an abstract source even with nothing to push against.
42 shell = body * shell;
44 wall = 1 - shell;
46 % Uniform air. The walls being a mask is what keeps this uniform — and the
47 % timestep as large as air allows.
48 c = c0;
50 % Open boundary at the domain edge, plus the wall's surface absorption.
51 sig = sponge3(x, y, z, Lx, Ly, Lz, 0.15*Ly, 9000) + absorb * (4 * shell .* (1 - shell));
53 % Where the string radiates directly: a Gaussian tube two cells wide
54 % around the string's line (y = 0, z = gap), fading out at its ends.
55 % Normalized so its cross-section integrates to one — the tube is a stand-in
56 % for a line, and without the 1/(pi a^2) its strength would depend on the
57 % grid that draws it.
58 zs = gap;
59 a = 2*h;
60 tube = exp(-(y.^2 + (z - zs).^2) / a^2) / (pi * a^2);
61 span = 0.5 * (1 + tanh((Ls/2 - abs(x)) / (1.5*h)));
62 lineprof = tube .* span .* wall;
64 % Where the bridge force drives the air: a patch of the air just above the
65 % top plate, around the bridge end of the string (x = +Ls/2). This stands
66 % in for the top plate moving; the plate's own modes are not modelled.
67 % Normalized to integrate to one over its volume, like the tube.
68 w = 1.5*h;
69 pr = exp(-((x - Ls/2).^2 + y.^2) / patchr^2) / (pi * patchr^2);
70 skin = exp(-(z - w).^2 / w^2) / (sqrt(pi) * w);
71 boardprof = pr .* skin .* wall;
72end