diff --git a/democapsid/src/democapsid.py b/democapsid/src/democapsid.py index dc8d318..02866e2 100755 --- a/democapsid/src/democapsid.py +++ b/democapsid/src/democapsid.py @@ -16,21 +16,37 @@ SQRT5 = np.sqrt(5) PHI = (1 + SQRT5) / 2 -# ⬢ lattice unit -HEXAGON = np.array([ +R6 = 1 +r6 = R6 * (SQRT3 / 2) + +HEXAGON_1 = np.array([ [0, 1], - [SQRT3 / 2, 0.5], - [SQRT3 / 2, -0.5], + [r6, 0.5], + [r6, -0.5], [0, -1], - [-(SQRT3 / 2), -0.5], - [-(SQRT3 / 2), 0.5] + [-r6, -0.5], + [-r6, 0.5] ] ) -# tilers -TILERS = [ - lambda coor: HEXAGON + coor +HEXAGON_2 = [ + [R6 / 2, r6], + [R6, 0], + [R6 / 2, -r6], + [-R6 / 2, -r6], + [-R6, 0], + [-R6 / 2, r6] ] -# icosahedron symmetry config +LATTICE = { + "hex": ( + lambda R: (R * (SQRT3 / 2)) * np.array([[2, 0], [1, SQRT3]]), + [lambda coor: HEXAGON + coor] + ), + "trihex": ( + lambda R: np.array([[2 * R, 0], [R, R * SQRT3]]), + [lambda coor: HEXAGON_2 + coor] + ) +} + ICO_CONFIG = ( (), (), @@ -55,7 +71,6 @@ ) ) - def iter_ring(elements): yield from ((elements[i], elements[(i + 1) % len(elements)]) for i in range(len(elements))) if len(elements) > 1 else () @@ -99,13 +114,6 @@ def in_triangle(p, q1, q2, q3): ) -def update_coor_to_id(coor_to_id, coors): - for coor in map(tuple, coors): - identifier = coor_to_id.get(coor, len(coor_to_id)) - coor_to_id[coor] = identifier - yield identifier - - def intersection(p1, q1, p2, q2): # http://paulbourke.net/geometry/pointlineplane/edge_intersection.py @@ -338,15 +346,14 @@ def main(argv): args = parse_args(argv[1:]) h, k, H, K, s, c = args.h, args.k, args.H, args.K, args.symmetry, args.sphericity - # hexagonal tile unit inradius - r = SQRT3 / 2 - # hexagonal tile lattice basis - basis = np.array([[2 * r, 0], [r, SQRT3 * r]]) + # tile + tile = "trihex" + + # lattice basis + basis = LATTICE[tile][0](1) v1, v2, v3 = (*basis, basis[1] @ rmat(np.pi / 3)) + # Caspar-Klug vectors - # self.q1 = np.array([0, 0]), self.C1, self.C4, self.C1 + self.C4 - # self.q2 = np.array([0, 0]), self.C2, self.C1, self.C2 + self.C1 - # self.q3 = np.array([0, 0]), self.C3, self.C2, self.C3 + self.C2 ckv = [ [h, k] @ basis, [H, K] @ np.stack([v2, v3]), @@ -354,9 +361,6 @@ def main(argv): [k, -h] @ np.stack([v1, v3]) ] # Caspar-Klug triangles - # self.C1 = h * self.a1 + k * self.a2 # \vec{C}_{T} - # self.C2 = H * self.A1 + K * self.A2 # \vec{C}_{Q} - # self.C3 = (-h - k) * self.a1 + h * self.a2 # \vec{C}^{120^{\circ}}_{T} ckt = [ [], [np.array([0, 0]), ckv[0], ckv[3]], @@ -379,8 +383,8 @@ def main(argv): mesh = [] for coor in lattice_coordinates: # process tile subunits - for tiler in TILERS: - path = list(iter_ring(tiler(coor @ basis))) + for calc_tile in LATTICE[tile][1]: + path = list(iter_ring(calc_tile(coor @ basis))) vertices = [] # iterate polygon edges for src, tar in path: @@ -391,15 +395,15 @@ def main(argv): # add point that at the intersetion of the polygon and triangle edges if (x := intersection(src, tar, *edge)).any(): vertices.append((np.append(x, 1), 1)) - # keep edges if they occuron the tile polygon path - edges = [ - (s1, t1) - for s1, t1 in iter_ring(list(range(len(vertices)))) - if any(on_same_line(vertices[s1][0], vertices[t1][0], np.append(s2, 1), np.append(t2, 1)) for s2, t2 in path) - ] - # if there are only two edges and they point to each other, then only keep one - edges = [edges[0]] if len(edges) == 2 and edges[0] == edges[1][::-1] else edges - edges and mesh.append(([ele[0] for ele in vertices], edges)) + # keep edges if they occur on the tile polygon path + edges = [ + (s1, t1) + for s1, t1 in iter_ring(list(range(len(vertices)))) + if any(on_same_line(vertices[s1][0], vertices[t1][0], np.append(s2, 1), np.append(t2, 1)) for s2, t2 in path) + ] + # if there are only two edges and they point to each other, then only keep one + edges = [edges[0]] if len(edges) == 2 and edges[0] == edges[1][::-1] else edges + edges and mesh.append(([ele[0] for ele in vertices], edges)) mesh and meshes.append(mesh) meshes3d = [[], [], [], []]