% Chaos that keeps its pitch: a small drum head bulged into a shallow dome
% by a steady load (like air pressure underneath), struck hard.
%
% A flat net only becomes chaotic when struck hard enough to bend its pitch
% by about a semitone, and then it sounds harsh. Curvature changes that:
% in a dome, two modes can drive a third directly (as in gongs and
% cymbals), so energy moves between modes at much smaller amplitude. Here
% the motion is chaotic (turn on "perturbed twin") while the main partial
% moves by only about 20 cents; what changes from strike to strike is the
% shimmer of the partials, not the pitch. Try force 5 (regular) and 15
% (chaotic), and compare with load = 0 (flat).

R = 0.06;          % radius (m)
a = 0.01;          % bead spacing (m)
sigma = 1;         % areal density (kg/m^2)
T = 2000;          % membrane tension (N/m)
strain = 0.005;    % spring stretch at rest

[pos, tri, rim] = bp_hexdisk(R, a);
E = bp_edges(tri);
f0 = T * a / sqrt(3);
k = f0 / (strain * a);

sys.pos = pos;
sys.mass = sigma * (sqrt(3)/2) * a^2;
sys.pinned = rim;
sys.springs = bp_springs(pos, E, k, f0);
[sys.friction, sys.dashpots] = bp_decay(pos, sys.mass, sys.springs, [200 2.5; 2000 1.5]);

% Steady upward force on every bead (N). The solver first finds the domed
% shape where springs and load balance (sys.relax, on by default with a
% load); the sag is reported with the run.
sys.load = [0 0 2];

hit.beads = bp_near(pos, [0.3*R 0.1*R 0], 1.5*a);
hit.force = [0 0 -15];       % peak force per bead (N)
hit.duration = 2e-3;
sys.hits = hit;

sys.listen = struct('beads', bp_near(pos, [-0.3*R 0.2*R 0], 0), 'dir', [0 0 1]);
sys.duration = 2.5;