% Drum head: a triangular net of beads under tension, pinned at the rim,
% struck by a soft mallet. Units are SI. The net stands in for a membrane
% with areal density sigma and tension T; its fundamental is near
% 2.405/(2*pi*R)*sqrt(T/sigma).

R = 0.15;          % radius (m)
a = 0.01;          % bead spacing (m)
sigma = 0.25;      % areal density (kg/m^2)
T = 2000;          % membrane tension (N/m)
strain = 0.02;     % spring stretch at rest. Smaller is closer to a real
                   % membrane and more nonlinear, but stiffer in-plane, so
                   % the time step shrinks.

[pos, tri, rim] = bp_hexdisk(R, a);
E = bp_edges(tri);
f0 = T * a / sqrt(3);        % spring tension that gives membrane tension T
k = f0 / (strain * a);       % spring stiffness

sys.pos = pos;
sys.mass = sigma * (sqrt(3)/2) * a^2;    % lattice area per bead
sys.pinned = rim;
sys.springs = bp_springs(pos, E, k, f0);
% Damping, set by decay times (T60: seconds to fall 60 dB) at two
% frequencies: low modes ring about 2.5 s, modes near 2 kHz about 0.75 s.
% bp_decay sets sys.friction and sys.dashpots from them; the page plots
% the decay time each mode gets.
sys = bp_decay(sys, [200 2.5; 2000 0.75]);

% A soft mallet: a 1 ms raised-cosine push on the beads within 3 cm.
hit.t = 0;
hit.beads = bp_near(pos, [0.05 0 0], 0.03);
hit.force = [0 0 -0.4];      % peak force per bead (N)
hit.duration = 1e-3;
sys.hits = hit;
% For more strikes, make sys.hits a struct array, e.g.
%   sys.hits(2) = struct('t', 0.5, 'beads', bp_near(pos, [0 0.06 0], 0.02), ...
%                        'force', [0 0 -0.6], 'duration', 1e-3);
% (every hit in the array needs the same fields). A hit may instead give an
% instant velocity: struct('t', 0, 'beads', idx, 'velocity', [0 0 -3]).

sys.listen = struct('beads', bp_near(pos, [-0.04 0.03 0], 0), 'dir', [0 0 1]);
sys.duration = 2;

fprintf('%d beads (%d pinned), %d springs\n', size(pos, 1), numel(rim), size(E, 1));
fprintf('expected fundamental: %.1f Hz\n', 2.405 / (2*pi*R) * sqrt(T / sigma));