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