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