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