function S = bp_springs(pos, E, k, f0, k3)
% BP_SPRINGS Springs along edges E with stiffness k and tension f0 at rest.
%   S = bp_springs(pos, E, k, f0) returns rows [i j k L0] for sys.springs,
%   choosing each rest length L0 so the spring pulls with force f0 at the
%   bead positions pos: L0 = L - f0/k. k and f0 may be scalars or one per edge.
%   S = bp_springs(pos, E, k, f0, k3) appends a cubic stiffness column.
L = vecnorm(pos(E(:,1), :) - pos(E(:,2), :), 2, 2);
m = size(E, 1);
k = k(:) .* ones(m, 1);
f0 = f0(:) .* ones(m, 1);
L0 = L - f0 ./ k;
if any(L0 <= 0)
    error('bp_springs: tension f0 is too large for stiffness k (rest length would be <= 0)');
end
S = [E, k, L0];
if nargin > 4
    S = [S, k3(:) .* ones(m, 1)];
end
end