1function [friction, dashpots, law] = bp_decay(pos, mass, springs, T60)
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% Each bead gets friction 2*sigma0*m, and each spring gets a dashpot
13% 2*sigma1*g in parallel, where g is the spring's tension divided by its
14% length at rest. The damping is then 2*sigma0*M + 2*sigma1*K for the
15% transverse motion, so the modes stay uncoupled and the law is exact for
16% them. Larger motion, and the in-plane modes, are damped by the same
17% friction and dashpots but not by this law.
18%
19% law = [sigma0 sigma1] is the third output.
20n = size(pos, 1);
21s = log(1000) ./ T60(:, 2); % amplitude decay rates (1/s)
22w2 = (2 * pi * T60(:, 1)).^2;
23if size(T60, 1) == 1
24 sigma0 = s(1);
25 sigma1 = 0;
26elseif size(T60, 1) == 2
27 sigma1 = (s(2) - s(1)) / (w2(2) - w2(1));
28 sigma0 = s(1) - sigma1 * w2(1);
29else
30 error('bp_decay: T60 must be [f1 T1] or [f1 T1; f2 T2]');
31end
32if sigma1 < 0
33 error('bp_decay: the higher frequency must not ring longer than the lower one');
34end
35if sigma0 < 0
36 error('bp_decay: decay at the higher frequency is too fast for sigma0 + sigma1*w^2; raise its T60 or lower the other');
37end
38friction = 2 * sigma0 * mass(:) .* ones(n, 1);
39i = springs(:, 1);
40j = springs(:, 2);
41L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
42e = L - springs(:, 4);
43f = springs(:, 3) .* e;
44if size(springs, 2) > 4
45 f = f + springs(:, 5) .* e.^3;
46end
47dashpots = [i, j, 2 * sigma1 * f ./ L];
48law = [sigma0, sigma1];
49end