% A bar with free ends, like a glockenspiel or xylophone bar, held softly
% (like cords) at the nodes of its lowest mode, 0.224 of the length from
% each end. Its stiffness comes from bending triplets, so its overtones are
% not harmonic: for an ideal free bar they sit at 2.76 and 5.40 times the
% fundamental. The soft supports add low "bouncing" modes of the whole bar
% below the fundamental. (Pinning the support beads instead holds them
% rigidly, which moves the overtones a long way.)
%
% The bending stiffness is set from the target pitch: an ideal free bar
% has f1 = 22.37/(2*pi*L^2) * sqrt(EI/rhoA), and the triplets give
% EI = kb*a^3 and rhoA = m/a.

N = 41;            % beads
L = 0.2;           % length (m)
rhoA = 0.2;        % mass per length (kg/m), e.g. 25 x 3 mm aluminium
f1 = 880;          % target fundamental (Hz)
thick = 0.15;      % kb/k: bending relative to stretching stiffness

[pos, E, trip] = bp_chain(N, L);
a = L / (N - 1);
m = rhoA * a;
kb = m / a^4 * (2*pi*f1 * L^2 / 22.37)^2;

sys.springs = bp_springs(pos, E, kb / thick, 0);
supports = [bp_near(pos, [0.224*L 0 0], 0), bp_near(pos, [0.776*L 0 0], 0)];
[pos, cords, anchors] = bp_anchor(pos, supports, 500, 0.015);   % 500 N/m each
sys.pos = pos;
sys.mass = m;
sys.pinned = anchors;
sys.springs = [sys.springs; cords];
sys.bends = [trip, kb * ones(size(trip, 1), 1)];
sys = bp_decay(sys, [f1 2.5; 10*f1 0.5]);

hit.t = 0;
hit.beads = bp_near(pos, [L/2 0 0], 0);
hit.force = [0 0 -40];       % N
hit.duration = 5e-4;
sys.hits = hit;

sys.listen = struct('beads', N - 2, 'dir', [0 0 1]);
sys.duration = 2.5;