function sys = bp_decay(sys, T60)
% BP_DECAY Set the damping from chosen decay times.
%   sys = bp_decay(sys, T60) sets sys.friction and sys.dashpots, and adds a
%   damping column to sys.bends and sys.hinges if present, such that at
%   small amplitude each mode 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). Call it after the elements are set.
%
%   Each bead gets friction 2*sigma0*m; each spring a dashpot in parallel
%   with isotropic part 2*sigma1*g and axial part 2*sigma1*(keff - g), where
%   g is the spring's tension divided by its length at rest and keff its
%   stiffness there (a spring's stiffness at rest is (keff - g) d d' + g I);
%   each bending triplet the damping 2*sigma1*kb, and each hinge 2*sigma1*kh.
%   The damping matrix at rest is then 2*sigma0*M + 2*sigma1*K for all the
%   motion, so the modes stay uncoupled and the law holds exactly at small
%   amplitude, for flat and curved nets alike.
%
%   The law [sigma0 sigma1] is stored in sys.decayLaw (not used by the solver).
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
pos = sys.pos;
n = size(pos, 1);
sys.friction = 2 * sigma0 * sys.mass(:) .* ones(n, 1);
if isfield(sys, 'springs') && ~isempty(sys.springs)
    S = sys.springs;
    i = S(:, 1);
    j = S(:, 2);
    L = vecnorm(pos(i, :) - pos(j, :), 2, 2);
    e = L - S(:, 4);
    f = S(:, 3) .* e;
    keff = S(:, 3);
    if size(S, 2) > 4
        f = f + S(:, 5) .* e.^3;
        keff = keff + 3 * S(:, 5) .* e.^2;
    end
    g = f ./ L;
    sys.dashpots = [i, j, 2 * sigma1 * g, 2 * sigma1 * (keff - g)];
else
    sys.dashpots = zeros(0, 4);
end
if isfield(sys, 'bends') && ~isempty(sys.bends)
    sys.bends = [sys.bends(:, 1:4), 2 * sigma1 * sys.bends(:, 4)];
end
if isfield(sys, 'hinges') && ~isempty(sys.hinges)
    sys.hinges = [sys.hinges(:, 1:5), 2 * sigma1 * sys.hinges(:, 5)];
end
sys.decayLaw = [sigma0, sigma1];
end