function [pos, tri, rim, trip] = 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.
%   trip bending triplets [i j k]: three beads in a row along one of the
%        three lattice directions, j in the middle (for sys.bends)
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);
% triplets along the three lattice directions: (1,0), (0,1), (1,-1) in (I,J)
t1 = [reshape(id(1:end-2, :), [], 1), reshape(id(2:end-1, :), [], 1), reshape(id(3:end, :), [], 1)];
t2 = [reshape(id(:, 1:end-2), [], 1), reshape(id(:, 2:end-1), [], 1), reshape(id(:, 3:end), [], 1)];
t3 = [reshape(id(1:end-2, 3:end), [], 1), reshape(id(2:end-1, 2:end-1), [], 1), reshape(id(3:end, 1:end-2), [], 1)];
trip = [t1; t2; t3];
trip = trip(all(trip > 0, 2), :);
end