function [friction, dashpots, law] = bp_decay(pos, mass, springs, T60)
% BP_DECAY Friction and dashpots that give chosen decay times.
%   [friction, dashpots] = bp_decay(pos, mass, springs, T60) returns values
%   for sys.friction and sys.dashpots such that, at small amplitude, each
%   transverse mode of a flat net decays at the amplitude rate
%       sigma(f) = sigma0 + sigma1 * (2*pi*f)^2,
%   passing through the decay times in T60 = [f1 T1; f2 T2] (Hz, seconds to
%   fall by 60 dB). With a single row [f1 T1], every mode gets the same
%   decay time. T1 must be at least T2 when f1 < f2 (higher modes cannot
%   ring longer under this law).
%
%   Each bead gets friction 2*sigma0*m, and each spring gets a dashpot
%   2*sigma1*g in parallel, where g is the spring's tension divided by its
%   length at rest. The damping is then 2*sigma0*M + 2*sigma1*K for the
%   transverse motion, so the modes stay uncoupled and the law is exact for
%   them. Larger motion, and the in-plane modes, are damped by the same
%   friction and dashpots but not by this law.
%
%   law = [sigma0 sigma1] is the third output.
n = size(pos, 1);
s = log(1000) ./ T60(:, 2);           % amplitude decay rates (1/s)
w2 = (2 * pi * T60(:, 1)).^2;
if size(T60, 1) == 1
    sigma0 = s(1);
    sigma1 = 0;
elseif size(T60, 1) == 2
    sigma1 = (s(2) - s(1)) / (w2(2) - w2(1));
    sigma0 = s(1) - sigma1 * w2(1);
else
    error('bp_decay: T60 must be [f1 T1] or [f1 T1; f2 T2]');
end
if sigma1 < 0
    error('bp_decay: the higher frequency must not ring longer than the lower one');
end
if sigma0 < 0
    error('bp_decay: decay at the higher frequency is too fast for sigma0 + sigma1*w^2; raise its T60 or lower the other');
end
friction = 2 * sigma0 * mass(:) .* ones(n, 1);
i = springs(:, 1);
j = springs(:, 2);
L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
e = L - springs(:, 4);
f = springs(:, 3) .* e;
if size(springs, 2) > 4
    f = f + springs(:, 5) .* e.^3;
end
dashpots = [i, j, 2 * sigma1 * f ./ L];
law = [sigma0, sigma1];
end