% A domed cymbal: a shallow bow with a raised bell, free at the rim and
% clamped at the top of the bell, struck sharply near the edge. Units are SI.
%
% Its bending stiffness comes from hinges across the interior edges of the
% triangle mesh (sys.hinges). A hinge resists changes of the angle between
% its two triangles from their angle at the rest positions, so the curved
% shape below rests without stress. The springs carry no load at rest and
% only resist stretching. Curvature couples bending to stretching already at
% small amplitude, which a flat plate lacks: set Hbow = 0 and Hbell = 0 to
% compare with the same plate flat.
%
% The linear modes of a curved net need all three coordinates of every free
% bead (here about 1100), which takes some seconds; set sys.modes = false to
% skip them.

R = 0.2;           % radius (m)
a = 0.02;          % bead spacing (m)
sigma = 8.9;       % mass per area (kg/m^2), about 1 mm of bronze
Hbow = 0.02;       % height of the bow (m)
Hbell = 0.02;      % extra height of the bell (m)
Rbell = 0.06;      % radius of the bell (m)
kh = 230;          % hinge stiffness (N m per rad^2); sets the pitch
k = 8e6;           % spring stiffness (N/m); resists stretching

[pos, tri] = bp_hexdisk(R, a);
r = sqrt(sum(pos(:, 1:2).^2, 2));
pos(:, 3) = Hbow * (1 - (r / R).^2) + Hbell * max(0, 1 - (r / Rbell).^2).^2;
E = bp_edges(tri);
hh = bp_hinges(tri);

sys.pos = pos;
sys.mass = sigma * (sqrt(3)/2) * a^2;
sys.pinned = bp_near(pos, [0 0 max(pos(:, 3))], 1.05 * a);   % the mount
sys.springs = bp_springs(pos, E, k, 0);                       % no tension at rest
sys.hinges = [hh, kh * ones(size(hh, 1), 1)];
sys = bp_decay(sys, [100 2.5; 5000 1.2]);

% A stick on the edge: a sharp 0.3 ms push, straight down, on the beads
% within 3 cm of a point near the rim.
hit.t = 0;
hit.beads = bp_near(pos, [0.17 0 0], 0.03);
hit.force = [0 0 -150];      % peak force per bead (N)
hit.duration = 3e-4;
sys.hits = hit;

sys.listen = struct('beads', bp_near(pos, [-0.1 0.06 0], 0.01), 'dir', [0 0 1]);
sys.duration = 2.5;

fprintf('%d beads, %d springs, %d hinges, %d struck\n', ...
    size(pos, 1), size(E, 1), size(hh, 1), numel(hit.beads));