function [pos, tri, rim] = bp_hexdisk(R, a)
% BP_HEXDISK Triangular lattice of beads with spacing a, clipped to a disk of radius R.
%   [pos, tri, rim] = bp_hexdisk(R, a)
%   pos  n-by-3 bead positions in the plane z = 0, centered at the origin
%   tri  triangles of the lattice (rows of three bead indices)
%   rim  beads with fewer than six lattice neighbors (the edge); pin these
%        for a drum head. Every other bead then sits at an equilibrium
%        under uniform spring tension.
m = ceil(R / a) + 2;
[I, J] = ndgrid(-2*m:2*m, -m:m);
X = a * (I + J/2);
Y = a * (sqrt(3)/2) * J;
inside = X.^2 + Y.^2 <= R^2 * (1 + 1e-12);
id = zeros(size(I));
id(inside) = 1:nnz(inside);
pos = [X(inside), Y(inside), zeros(nnz(inside), 1)];
A = id(1:end-1, 1:end-1);
B = id(2:end, 1:end-1);
C = id(1:end-1, 2:end);
D = id(2:end, 2:end);
tri = [A(:), B(:), C(:); B(:), D(:), C(:)];
tri = tri(all(tri > 0, 2), :);
E = bp_edges(tri);
deg = accumarray(E(:), 1, [size(pos, 1), 1]);
rim = find(deg < 6);
end