A few updates. After trying to make a 3d-printable surface, I realized the electrostatic potential was the wrong kind of terrain and switched to Gaussian drop-off and then exponential, which Claude helpfully suggested as Yukawa potential.
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("--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)
# 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) # 36
n_B = len(group_B) # 40
q_A = -float(n_B) / n_A # ≈ -1.1111, keeps total charge = 0
q_B = 1.0
positions = np.array(group_A + group_B)
charges = np.array([q_A] * n_A + [q_B] * n_B)
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
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.two_d:
kernel = np.log(np.maximum(np.sqrt(r2), r_min))
else:
kernel = 1.0 / np.maximum(np.sqrt(r2), r_min)
V += charges[i] * kernel
mode = (f"Gaussian σ={args.gaussian}" if args.gaussian is not None
else f"Yukawa λ={args.yukawa}" if args.yukawa 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:
v_lo = V.min()
v_hi = V.max()
else:
v_lo = np.percentile(V, 2)
v_hi = np.percentile(V, 98)
v_rng = max(abs(v_lo), abs(v_hi))
levels = np.linspace(-v_rng, v_rng, 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}")
You might wonder why I'm doing the 2-coloring. There's not a clear justification except I initially imagined this as electrostatic charge and wanted it to zero-out at a distance. But with faster attenuation, there's not as much need. I'll try that later. Another idea is to have positive potential and take max, not sum. In that case, I believe the voronoi boundaries will show up naturally. I'll try these out later today.
Here are vertices expanded out 5 levels.
