1function [friction, dashpots, bends, law] = bp_decay(pos, mass, springs, T60, bends)
2% BP_DECAY Friction and dashpots that give chosen decay times.
3% [friction, dashpots] = bp_decay(pos, mass, springs, T60) returns values
4% for sys.friction and sys.dashpots such that, at small amplitude, each
5% transverse mode of a flat net decays at the amplitude rate
6% sigma(f) = sigma0 + sigma1 * (2*pi*f)^2,
7% passing through the decay times in T60 = [f1 T1; f2 T2] (Hz, seconds to
8% fall by 60 dB). With a single row [f1 T1], every mode gets the same
9% decay time. T1 must be at least T2 when f1 < f2 (higher modes cannot
10% ring longer under this law).
11%
12% [friction, dashpots, bends] = bp_decay(pos, mass, springs, T60, bends)
13% also damps the bending triplets (rows [i j k kb]): the returned bends
14% get a fifth column cb, for sys.bends. Pass the bends whenever the
15% system has them, or plate modes will not follow the law.
16%
17% Each bead gets friction 2*sigma0*m, each spring a dashpot 2*sigma1*g in
18% parallel (g is the spring's tension divided by its length at rest), and
19% each bending triplet the damping cb = 2*sigma1*kb. The damping is then
20% 2*sigma0*M + 2*sigma1*K for the transverse motion, so the modes stay
21% uncoupled and the law is exact for them. Larger motion, and the in-plane
22% modes, are damped by the same elements but not by this law.
23%
24% law = [sigma0 sigma1] is the fourth output.
25n = size(pos, 1);
26s = log(1000) ./ T60(:, 2); % amplitude decay rates (1/s)
27w2 = (2 * pi * T60(:, 1)).^2;
28if size(T60, 1) == 1
29 sigma0 = s(1);
30 sigma1 = 0;
31elseif size(T60, 1) == 2
32 sigma1 = (s(2) - s(1)) / (w2(2) - w2(1));
33 sigma0 = s(1) - sigma1 * w2(1);
34else
35 error('bp_decay: T60 must be [f1 T1] or [f1 T1; f2 T2]');
36end
37if sigma1 < 0
38 error('bp_decay: the higher frequency must not ring longer than the lower one');
39end
40if sigma0 < 0
41 error('bp_decay: decay at the higher frequency is too fast for sigma0 + sigma1*w^2; raise its T60 or lower the other');
42end
43friction = 2 * sigma0 * mass(:) .* ones(n, 1);
44i = springs(:, 1);
45j = springs(:, 2);
46L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
47e = L - springs(:, 4);
48f = springs(:, 3) .* e;
49if size(springs, 2) > 4
50 f = f + springs(:, 5) .* e.^3;
51end
52dashpots = [i, j, 2 * sigma1 * f ./ L];
53if nargin > 4
54 bends = [bends(:, 1:4), 2 * sigma1 * bends(:, 4)];
55else
56 bends = zeros(0, 5);
57end
58law = [sigma0, sigma1];
59end