1function [pos, springs, anchors] = bp_anchor(pos, beads, kz, h)
2% BP_ANCHOR Hold beads softly, like a cord or a felt, rather than pinning them.
3% [pos, springs, anchors] = bp_anchor(pos, beads, kz, h) adds four fixed
4% anchor beads around each listed bead, at distance h (default 0.02 m)
5% along +x, -x, +y and -y in the bead's own plane, and returns
6% pos the positions with the anchors appended
7% springs rows [i j k L0] for the new springs (append to sys.springs)
8% anchors indices of the new beads (append to sys.pinned)
9% The springs are pre-tensioned so that each held bead feels a stiffness
10% kz (N/m) against motion out of the plane, and 2*kz against motion in
11% the plane. The forces balance at rest, and the anchors lie in the
12% object's plane, so a flat object keeps its linear-mode analysis.
13if nargin < 4
14 h = 0.02;
15end
16beads = beads(:);
17nb = numel(beads);
18n0 = size(pos, 1);
19dirs = [1 0 0; -1 0 0; 0 1 0; 0 -1 0];
20f0 = kz * h / 4; % four springs of tension f0 give kz out of plane
21newpos = zeros(4 * nb, 3);
22springs = zeros(4 * nb, 4);
23for q = 1:nb
24 for d = 1:4
25 r = 4 * (q - 1) + d;
26 newpos(r, :) = pos(beads(q), :) + h * dirs(d, :);
27 springs(r, :) = [beads(q), n0 + r, kz, h - f0 / kz];
28 end
29end
30pos = [pos; newpos];
31anchors = (n0 + 1:n0 + 4 * nb)';
32end