After fiddling around with different ideas, I'm back to voronoi diagrams, this time using the same code with a few enhancements. Using a spherical kernel and max instead of plus, I get the upper envelope of spheres centered at vertices, showing the Penrose tile as conjoined bubbles, and automatically revealing the Voronoi diagram at boundaries (all too familiar to me after working on many variations of these seven cell shapes).
Run this code with options 11 2000 --positive --max --sphere 0.5 --heatmap --bw --no-dots --pts pen5.txt
Code: Select all
import math
import numpy as np
from PIL import Image, ImageDraw
from collections import deque
# Parameters
# Usage: script.py [k [size [--pts file.txt]]]
import argparse, sys
parser = argparse.ArgumentParser()
parser.add_argument("k", nargs="?", type=int, default=5, help="number of equipotential lines")
parser.add_argument("size", nargs="?", type=int, default=800, help="image size in pixels")
parser.add_argument("--pts", metavar="FILE", default=None, help="text file with x y coords, one per line")
parser.add_argument("--2d", dest="two_d", action="store_true", help="use 2D log potential (ln r) instead of 3D (1/r)")
parser.add_argument("--gaussian", metavar="SIGMA", type=float, default=None,
help="use Gaussian potential exp(-r²/2σ²) with given sigma instead of electrostatic")
parser.add_argument("--yukawa", metavar="LAMBDA", type=float, default=None,
help="use exponential potential exp(-r/λ) with given decay length")
parser.add_argument("--sphere", metavar="R", type=float, default=None,
help="use spherical cap kernel sqrt(1-(r/R)²) clamped to 0 beyond radius R")
parser.add_argument("--positive", action="store_true",
help="set all charges to +1 (skips bipartite colouring; best with --gaussian or --yukawa)")
parser.add_argument("--max", dest="use_max", action="store_true",
help="use max instead of sum over charges (Voronoi-like ridges)")
parser.add_argument("--edges", action="store_true", help="draw adjacency edges of the unit-distance graph")
parser.add_argument("--heatmap", action="store_true", help="render potential as a colour heatmap instead of contour lines")
parser.add_argument("--no-dots", dest="no_dots", action="store_true", help="do not draw charge position dots")
parser.add_argument("--bw", action="store_true", help="greyscale heightmap (black=min, white=max) for 3D rendering")
args = parser.parse_args()
k = args.k
size = args.size
r_min = 0.01 # clamp radius to avoid singularity
margin = 0.5 # margin around point cloud bounding box
# ------------------------------------------------------------
# All vertices of the Penrose rhomb tiling (unit edge length).
# We build the unit-distance graph and 2-colour it by BFS,
# then assign negative charge to one colour class and positive
# to the other so the total charge remains zero.
# ------------------------------------------------------------
all_pts = [
# ---- originally labelled "black" (51 pts) ----
(0.000000, 0.000000), (1.000000, 0.000000), (0.500000, -0.363271),
(0.309017, -0.951057), (-0.500000, -0.363271), (-0.690983, -0.951057),
(-0.500000, -1.538842), (0.000000, -1.902113), (0.500000, -1.538842),
(1.000000, -1.902113), (0.190983, -2.489898), (0.690983, -2.853170),
(1.309017, -2.853170), (1.809017, -2.489898), (1.618034, -1.902113),
(2.118034, -1.538842), (2.427051, -2.489898), (2.927051, -2.126627),
(3.118034, -1.538842), (2.927051, -0.951057), (2.309017, -0.951057),
(2.118034, -0.363271), (3.118034, -0.363271), (2.927051, 0.224514),
(2.427051, 0.587785), (1.809017, 0.587785), (1.618034, -0.000000),
(1.309017, 0.951057), (0.690983, 0.951057), (0.190983, 0.587785),
(-0.809017, 0.587785), (-0.000000, 1.175571), (0.500000, 1.538842),
(1.118034, 1.538842), (2.118034, 1.538842), (2.927051, 0.951057),
(3.427051, 0.587785), (3.618034, -0.000000), (3.927051, -0.951057),
(3.618034, -1.902113), (3.427051, -2.489898), (2.927051, -2.853170),
(2.118034, -3.440955), (1.118034, -3.440955), (0.500000, -3.440955),
(0.000000, -3.077684), (-0.809017, -2.489898), (-1.118034, -1.538842),
(-1.309017, -0.951057), (-1.118034, -0.363271), (1.309017, -0.951057),
# ---- originally labelled "red" (25 pts) ----
(-0.190983, -0.587785), (-0.190983, -1.314328), (0.500000, -2.265384),
(1.190983, -2.489898), (2.309017, -2.126627), (2.736068, -1.538842),
(2.736068, -0.363271), (2.309017, 0.224514), (1.190983, 0.587785),
(0.500000, 0.363271), (-0.190983, 0.587785), (1.618034, 1.175571),
(2.309017, 0.951057), (3.427051, -0.587785), (3.427051, -1.314328),
(2.309017, -2.853170), (1.618034, -3.077684), (-0.190983, -2.489898),
(-0.618034, -1.902113), (-0.618034, -0.000000), (1.118034, -0.363271),
(0.690983, -0.951057), (1.118034, -1.538842), (1.809017, -1.314328),
(1.809017, -0.587785),
]
if args.pts:
with open(args.pts) as f:
all_pts = [tuple(float(v) for v in line.split()) for line in f if line.strip()]
print(f"Loaded {len(all_pts)} points from {args.pts}")
pts_arr = np.array(all_pts)
n = len(pts_arr)
# Build unit-distance adjacency (edge length = 1/phi ≈ 0.618, tolerance 0.05)
EDGE_LEN = 2 / (1 + np.sqrt(5)) # 1/phi = phi - 1 ≈ 0.61803
adj = [[] for _ in range(n)]
for i in range(n):
for j in range(i + 1, n):
if abs(np.linalg.norm(pts_arr[i] - pts_arr[j]) - EDGE_LEN) < 0.05:
adj[i].append(j)
adj[j].append(i)
if args.positive:
# All charges +1; skip bipartite colouring entirely
positions = pts_arr
charges = np.ones(n)
adj_final = adj # adjacency already in pts_arr order
n_A = 0 # no group A; all dots drawn as blue
else:
# BFS 2-colouring (graph is bipartite — verified)
color = [-1] * n
color[0] = 0
q = deque([0])
while q:
node = q.popleft()
for nb in adj[node]:
if color[nb] == -1:
color[nb] = 1 - color[node]
q.append(nb)
group_A_idx = [i for i in range(n) if color[i] == 0]
group_B_idx = [i for i in range(n) if color[i] == 1]
group_A = [all_pts[i] for i in group_A_idx] # negative (red dots)
group_B = [all_pts[i] for i in group_B_idx] # positive (blue dots)
# Remap adjacency to positions-array ordering (group_A first, then group_B)
old_to_new = {old_i: new_i for new_i, old_i in enumerate(group_A_idx + group_B_idx)}
adj_final = [[] for _ in range(n)]
for old_i, neighbors in enumerate(adj):
new_i = old_to_new[old_i]
adj_final[new_i] = [old_to_new[nb] for nb in neighbors]
n_A = len(group_A)
n_B = len(group_B)
q_A = -float(n_B) / n_A
q_B = 1.0
positions = np.array(group_A + group_B)
charges = np.array([q_A] * n_A + [q_B] * n_B)
if args.positive:
print(f"All charges +1: {len(positions)} points")
else:
print(f"Group A (negative): {n_A} x {q_A:.4f} = {n_A * q_A:.4f}")
print(f"Group B (positive): {n_B} x {q_B:.4f} = {n_B * q_B:.4f}")
print(f"Total charge: {charges.sum():.6f}")
# Bounding box + margin
x_lo = positions[:, 0].min() - margin
x_hi = positions[:, 0].max() + margin
y_lo = positions[:, 1].min() - margin
y_hi = positions[:, 1].max() + margin
xs = np.linspace(x_lo, x_hi, size)
ys = np.linspace(y_lo, y_hi, size)
X, Y = np.meshgrid(xs, ys)
# Compute potential
if args.use_max:
V = np.full((size, size), -np.inf)
else:
V = np.zeros((size, size))
for i in range(len(positions)):
dx = X - positions[i, 0]
dy = Y - positions[i, 1]
r2 = dx**2 + dy**2
if args.gaussian is not None:
s2 = args.gaussian ** 2
kernel = np.exp(-r2 / (2 * s2))
elif args.yukawa is not None:
kernel = np.exp(-np.sqrt(r2) / args.yukawa)
elif args.sphere is not None:
t = np.clip(1.0 - r2 / args.sphere ** 2, 0.0, None)
kernel = np.sqrt(t)
elif args.two_d:
kernel = np.log(np.maximum(np.sqrt(r2), r_min))
else:
kernel = 1.0 / np.maximum(np.sqrt(r2), r_min)
contribution = charges[i] * kernel
if args.use_max:
V = np.maximum(V, contribution)
else:
V += contribution
mode = (f"Gaussian σ={args.gaussian}" if args.gaussian is not None
else f"Yukawa λ={args.yukawa}" if args.yukawa is not None
else f"Spherical cap R={args.sphere}" if args.sphere is not None
else "2D logarithmic" if args.two_d else "3D Coulomb (1/r)")
print(f"Potential: {mode}")
# Potential range: full extent for Gaussian, percentile clip otherwise
if args.gaussian is not None or args.yukawa is not None or args.sphere is not None:
v_lo = V.min()
v_hi = V.max()
else:
v_lo = np.percentile(V, 2)
v_hi = np.percentile(V, 98)
# Use symmetric levels if potential straddles zero, direct range otherwise
if v_lo * v_hi < 0:
v_rng = max(abs(v_lo), abs(v_hi))
levels = np.linspace(-v_rng, v_rng, k + 2)[1:-1]
else:
levels = np.linspace(v_lo, v_hi, k + 2)[1:-1]
print(f"Potential range: {v_lo:.3f} to {v_hi:.3f}")
print(f"Equipotential levels: {[f'{l:.3f}' for l in levels]}")
# Render
if args.bw:
# Greyscale heightmap: black=min, white=max
V_clip = np.clip(V, v_lo, v_hi)
t = (V_clip - v_lo) / (v_hi - v_lo) # 0..1
grey = (t * 255).astype(np.uint8)
img = Image.fromarray(grey, "L").convert("RGB")
elif args.heatmap:
# Diverging heatmap: blue (negative) → white (zero) → red (positive)
V_clip = np.clip(V, v_lo, v_hi)
t = (V_clip - v_lo) / (v_hi - v_lo) # 0..1
# blue→white for t<0.5, white→red for t>0.5
R = np.where(t < 0.5, t * 2, 1.0)
G = np.where(t < 0.5, t * 2, 2.0 - t * 2)
B = np.where(t < 0.5, 1.0, 2.0 - t * 2)
rgb = np.stack([R, G, B], axis=-1)
img = Image.fromarray((rgb * 255).astype(np.uint8), "RGB")
else:
img = Image.new("RGB", (size, size), "white")
pixels = img.load()
line_color = (20, 20, 20)
for level in levels:
shifted = V - level
iy, ix = np.where(shifted[:-1, :] * shifted[1:, :] < 0)
for y, x in zip(iy, ix):
pixels[x, y] = line_color
iy, ix = np.where(shifted[:, :-1] * shifted[:, 1:] < 0)
for y, x in zip(iy, ix):
pixels[x, y] = line_color
# Draw charge positions
draw = ImageDraw.Draw(img)
r_dot = size / 300
def to_pixel(px, py):
ix = int((px - x_lo) / (x_hi - x_lo) * size)
iy = int((py - y_lo) / (y_hi - y_lo) * size)
return ix, iy
# Draw adjacency edges (drawn first so dots render on top)
if args.edges:
edge_color = (180, 180, 180)
for i in range(len(positions)):
for j in adj_final[i]:
if j > i:
x1, y1 = to_pixel(*positions[i])
x2, y2 = to_pixel(*positions[j])
draw.line([x1, y1, x2, y2], fill=edge_color, width=1)
if not args.no_dots:
for i, (px, py) in enumerate(positions):
cx, cy = to_pixel(px, py)
# Colour by graph 2-colouring class, not by charge sign
dot_color = (200, 30, 30) if i < n_A else (30, 30, 200) # red = group A, blue = group B
draw.ellipse([cx - r_dot, cy - r_dot, cx + r_dot, cy + r_dot],
fill=dot_color, outline=(0, 0, 0))
out_path = 'equipotential_penrose_3d.png'
img.save(out_path)
img.show()
print(f"Saved to {out_path}")
Here's the B/W intensity map.
The cartwheels jump out visually, and you can find the rhomb tiling by connecting centers of adjacent bubbles, omitting the internal diagonals.
I don't plan to spend 4 hours and 10m of filament on this, but here it is in the slicer, 12cm diameter, at a high enough resolution to capture detail. I'll start with some smaller ones. I scaled it vertically and these are ellipsoids, not spheres.
I could restore the smooth surface by coating the mold and get a lens-like effect in clear resin. I might try this weekend, probably with a smaller layout.