lamp: paint ξ in the u-chart (1920×1080, 2× supersampled) and trace the |ξ| = 1 cone edge

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FAqjv3n4yABxvJPJHFRhVA

src/lamp/xi_image.py

@@ -0,0 +1,51 @@
+# Paint xi in the u-chart, u = 1/(4(s - 1/2)), as a 1920x1080 PNG (supersampled 2x), and trace the |xi| = 1 contour (the edge of the dark cone).
+import numpy as np, json
+from scipy.special import loggamma
+from PIL import Image
+W, H, SS = 1920, 1080, 2
+HW_Y = 0.016                     # half-height of the window in u
+PPU = H*SS/(2*HW_Y)              # pixels per unit of u (supersampled)
+BC = [1/12,-1/720,1/30240,-1/1209600,1/47900160,-691/1307674368000,1/74724249600,-3617/10670622842880000,43867/5109094217170944000,-174611/802857662698291200000]
+N = 160
+def log_zeta(s):                 # Euler-Maclaurin, valid here because every s in view has Re s >= 1/2 after reflection
+    n = np.arange(1, N, dtype=float)
+    tot = np.zeros_like(s)
+    for k in n: tot += k**(-s)
+    tot += N**(1-s)/(s-1) + 0.5*N**(-s)
+    p = s.copy(); sc = N**(-s)/N
+    for k, c in enumerate(BC, 1):
+        tot += c*p*sc; p = p*(s+2*k-1)*(s+2*k); sc = sc/(N*N)
+    return np.log(tot)
+def log_xi(s):
+    s = np.where(s.real < 0.5, 1-s, s)
+    return np.log(0.5*s*(s-1)) - s/2*np.log(np.pi) + loggamma(s/2) + log_zeta(s)
+ys, xs = np.mgrid[0:H*SS, 0:W*SS]
+ux = (xs - W*SS/2 + 0.5)/PPU; uy = -(ys - H*SS/2 + 0.5)/PPU
+u = ux + 1j*uy
+r2 = ux*ux + uy*uy
+s = 0.5 + 0.25*np.conj(u)/np.maximum(r2, 1e-12)      # 1/(4u) = conj(u)/(4|u|^2)
+L = np.full(u.shape, np.nan+0j)
+ok = r2 > (1/(4*900))**2                              # |s - 1/2| < 900: within the Euler-Maclaurin budget
+L[ok] = log_xi(s[ok])
+mag, ph = L.real, L.imag
+h = (ph/(2*np.pi)) % 1.0
+v = 0.14 + 0.74*(0.5 + np.arctan(0.7*mag/2.5)/np.pi)
+v = np.where(ok, v, 0.16); h = np.where(ok, h, 0.0); sat = np.where(ok, 0.62, 0.0)
+# hsv -> rgb
+i = np.floor(h*6).astype(int) % 6; f = h*6 - np.floor(h*6)
+p = v*(1-sat); q = v*(1-f*sat); t = v*(1-(1-f)*sat)
+R = np.choose(i, [v, q, p, p, t, v]); G = np.choose(i, [t, v, v, q, p, p]); B = np.choose(i, [p, p, t, v, v, q])
+img = np.stack([R, G, B], -1)
+img = img.reshape(H, SS, W, SS, 3).mean(axis=(1, 3))          # downsample
+Image.fromarray((img*255).clip(0, 255).astype(np.uint8)).save('xi.png')
+# cone edge: for each row, the first column right of the strut where |xi| crosses 1, in u-units
+m = mag.reshape(H, SS, W, SS).mean(axis=(1, 3))
+pts = []
+cx = W/2
+for row in range(0, H, 6):
+    line = m[row, int(cx):]
+    idx = np.argmax(line > 0) if np.any(line > 0) else -1
+    if idx > 0:
+        pts.append([round(idx/(PPU/SS), 6), round(-(row - H/2 + 0.5)/(PPU/SS), 6)])
+json.dump({'hw_y': HW_Y, 'hw_x': HW_Y*W/H, 'edge': pts}, open('cone.json', 'w'))
+print('done', img.shape, len(pts))