1% Drum head: a triangular net of beads under tension, pinned at the rim,
2% struck by a soft mallet. Units are SI. The net stands in for a membrane
3% with areal density sigma and tension T; its fundamental is near
4% 2.405/(2*pi*R)*sqrt(T/sigma).
6R = 0.15; % radius (m)
7a = 0.01; % bead spacing (m)
8sigma = 0.25; % areal density (kg/m^2)
9T = 2000; % membrane tension (N/m)
10strain = 0.02; % spring stretch at rest. Smaller is closer to a real
11 % membrane and more nonlinear, but stiffer in-plane, so
12 % the time step shrinks.
14[pos, tri, rim] = bp_hexdisk(R, a);
15E = bp_edges(tri);
16f0 = T * a / sqrt(3); % spring tension that gives membrane tension T
17k = f0 / (strain * a); % spring stiffness
19sys.pos = pos;
20sys.mass = sigma * (sqrt(3)/2) * a^2; % lattice area per bead
21sys.pinned = rim;
22sys.springs = bp_springs(pos, E, k, f0);
23% Damping, set by decay times (T60: seconds to fall 60 dB) at two
24% frequencies: low modes ring about 2.5 s, modes near 2 kHz about 0.75 s.
25% bp_decay sets sys.friction and sys.dashpots from them; the page plots
26% the decay time each mode gets.
27sys = bp_decay(sys, [200 2.5; 2000 0.75]);
29% A soft mallet: a 1 ms raised-cosine push on the beads within 3 cm.
30hit.t = 0;
31hit.beads = bp_near(pos, [0.05 0 0], 0.03);
32hit.force = [0 0 -0.4]; % peak force per bead (N)
33hit.duration = 1e-3;
34sys.hits = hit;
35% For more strikes, make sys.hits a struct array, e.g.
36% sys.hits(2) = struct('t', 0.5, 'beads', bp_near(pos, [0 0.06 0], 0.02), ...
37% 'force', [0 0 -0.6], 'duration', 1e-3);
38% (every hit in the array needs the same fields). A hit may instead give an
39% instant velocity: struct('t', 0, 'beads', idx, 'velocity', [0 0 -3]).
41sys.listen = struct('beads', bp_near(pos, [-0.04 0.03 0], 0), 'dir', [0 0 1]);
42sys.duration = 2;
44fprintf('%d beads (%d pinned), %d springs\n', size(pos, 1), numel(rim), size(E, 1));
45fprintf('expected fundamental: %.1f Hz\n', 2.405 / (2*pi*R) * sqrt(T / sigma));