1function [pos, tri, rim, trip] = 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.
9% trip bending triplets [i j k]: three beads in a row along one of the
10% three lattice directions, j in the middle (for sys.bends)
11m = ceil(R / a) + 2;
12[I, J] = ndgrid(-2*m:2*m, -m:m);
13X = a * (I + J/2);
14Y = a * (sqrt(3)/2) * J;
15inside = X.^2 + Y.^2 <= R^2 * (1 + 1e-12);
16id = zeros(size(I));
17id(inside) = 1:nnz(inside);
18pos = [X(inside), Y(inside), zeros(nnz(inside), 1)];
19A = id(1:end-1, 1:end-1);
20B = id(2:end, 1:end-1);
21C = id(1:end-1, 2:end);
22D = id(2:end, 2:end);
23tri = [A(:), B(:), C(:); B(:), D(:), C(:)];
24tri = tri(all(tri > 0, 2), :);
25E = bp_edges(tri);
26deg = accumarray(E(:), 1, [size(pos, 1), 1]);
27rim = find(deg < 6);
28% triplets along the three lattice directions: (1,0), (0,1), (1,-1) in (I,J)
29t1 = [reshape(id(1:end-2, :), [], 1), reshape(id(2:end-1, :), [], 1), reshape(id(3:end, :), [], 1)];
30t2 = [reshape(id(:, 1:end-2), [], 1), reshape(id(:, 2:end-1), [], 1), reshape(id(:, 3:end), [], 1)];
31t3 = [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)];
32trip = [t1; t2; t3];
33trip = trip(all(trip > 0, 2), :);
34end