% The dulcimer's body: a rigid box under the string, with a round sound hole
% in its top plate.
%
% The top plate lies at z = 0. The box hangs below it (down to z = -boxd) and
% the string runs along x just above it, at height `gap`. Everything here is
% metres and inverse seconds.
%
% The shell is a *mask*, not a material: `wall` is 0 inside the solid and 1
% in air, and the solver's masked Laplacian drops any flux into a masked
% cell, which is the rigid (Neumann) condition. So the walls reflect
% everything and cost nothing — no timestep penalty, unlike a fast material.
% What a rigid mask cannot do is absorb, so `absorb` paints a thin layer of
% ordinary volume absorption over the shell's surfaces (the smoothed mask's
% transition skin, 4*wall*(1-wall), which peaks exactly at the surface); that
% is what lets the cavity's ring be tamed like real wood tames it.
%
% The masks are drawn with a sharp profile (about half a cell) on purpose. A
% box is grid-aligned, so sharpness costs no staircase artefacts, and a plate
% only two cells thick has to reach mask ≈ 0 at its centre or it leaks.
function [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)
  s = 0.45 * h;

  % The solid: the outer box minus its interior cavity.
  ox = 0.5 * (1 - tanh((abs(x) - boxl/2) / s));
  oy = 0.5 * (1 - tanh((abs(y) - boxw/2) / s));
  oz = 0.5 * (1 - tanh((abs(z + boxd/2) - boxd/2) / s));
  ix = 0.5 * (1 - tanh((abs(x) - (boxl/2 - thick)) / s));
  iy = 0.5 * (1 - tanh((abs(y) - (boxw/2 - thick)) / s));
  iz = 0.5 * (1 - tanh((abs(z + boxd/2) - (boxd/2 - thick)) / s));
  shell = max(0, ox .* oy .* oz - ix .* iy .* iz);

  % The sound hole: a cylinder of radius holer at (holex, 0), cut out of the
  % top plate only.
  r = sqrt((x - holex).^2 + y.^2);
  hole = 0.5 * (1 - tanh((r - holer) / s));
  top = 0.5 * (1 - tanh((abs(z + thick/2) - thick/2) / s));
  shell = max(0, shell - hole .* top);

  % `body` at 0 takes the box away entirely: a bare string in open air. The
  % bridge patch stays where the plate would have been, so the bridge route
  % keeps working as an abstract source even with nothing to push against.
  shell = body * shell;

  wall = 1 - shell;

  % Uniform air. The walls being a mask is what keeps this uniform — and the
  % timestep as large as air allows.
  c = c0;

  % Open boundary at the domain edge, plus the wall's surface absorption.
  sig = sponge3(x, y, z, Lx, Ly, Lz, 0.15*Ly, 9000) + absorb * (4 * shell .* (1 - shell));

  % Where the string radiates directly: a Gaussian tube two cells wide
  % around the string's line (y = 0, z = gap), fading out at its ends.
  % Normalized so its cross-section integrates to one — the tube is a stand-in
  % for a line, and without the 1/(pi a^2) its strength would depend on the
  % grid that draws it.
  zs = gap;
  a = 2*h;
  tube = exp(-(y.^2 + (z - zs).^2) / a^2) / (pi * a^2);
  span = 0.5 * (1 + tanh((Ls/2 - abs(x)) / (1.5*h)));
  lineprof = tube .* span .* wall;

  % Where the bridge force drives the air: a patch of the air just above the
  % top plate, around the bridge end of the string (x = +Ls/2). This stands
  % in for the top plate moving; the plate's own modes are not modelled.
  % Normalized to integrate to one over its volume, like the tube.
  w = 1.5*h;
  pr = exp(-((x - Ls/2).^2 + y.^2) / patchr^2) / (pi * patchr^2);
  skin = exp(-(z - w).^2 / w^2) / (sqrt(pi) * w);
  boardprof = pr .* skin .* wall;
end