#!/usr/bin/env python3 """ HPC-030: Wigner-Seitz Cell Internal Resonance — EM FDTD ========================================================= Date: 2026-03-28 Depends on: HPC-026/027 (EM concentration by geometry), cipher archetypes Tests the Wigner-Seitz cells (Voronoi cells) of BCC, FCC, HCP lattices as EM cavities. Measures whether internal resonance patterns correspond to known material properties of each archetype. Key insight: atoms are the walls, voids are the cavities. The cipher classifies the external lattice; this test measures the internal resonance. Engine: EMFDTD3D (Yee FDTD, Maxwell's curl equations) Fields: Ex, Ey, Ez, Hx, Hy, Hz (6-component vector) Usage: python3 HPC-030_wigner_seitz_resonance.py python3 HPC-030_wigner_seitz_resonance.py --grid 96 --steps 2000 """ import argparse import json import math import time as time_mod from datetime import datetime, timezone from pathlib import Path import numpy as np C_LIGHT = 3.0e8 MU_0 = 4 * np.pi * 1e-7 EPS_0 = 8.854e-12 PHI = (1 + math.sqrt(5)) / 2 RESULTS_DIR = Path(__file__).parent / "results_hpc030" RESULTS_DIR.mkdir(exist_ok=True) # ============================================================================== # EM FDTD ENGINE — Yee Algorithm (Maxwell's Equations) # Same validated engine as HPC-025/026/027/028/029 # ============================================================================== class EMFDTD3D: """3D EM FDTD — Yee algorithm solving Maxwell's curl equations.""" def __init__(self, N, dx, dt=None): self.N = N self.dx = dx self.dt = dt or (0.99 * dx / (C_LIGHT * np.sqrt(3))) self.eps_r = np.ones((N, N, N), dtype=np.float32) self.sigma = np.zeros((N, N, N), dtype=np.float32) self.Ex = np.zeros((N, N, N), dtype=np.float32) self.Ey = np.zeros((N, N, N), dtype=np.float32) self.Ez = np.zeros((N, N, N), dtype=np.float32) self.Hx = np.zeros((N, N, N), dtype=np.float32) self.Hy = np.zeros((N, N, N), dtype=np.float32) self.Hz = np.zeros((N, N, N), dtype=np.float32) self._compute_coefficients() def _compute_coefficients(self): dt = self.dt eps = self.eps_r * EPS_0 self.Ca = (1 - self.sigma * dt / (2 * eps)) / (1 + self.sigma * dt / (2 * eps)) self.Cb = (dt / eps) / (1 + self.sigma * dt / (2 * eps)) / self.dx def set_conductor(self, mask): self.sigma[mask] = 1e7 self._compute_coefficients() def step(self): dt, dx = self.dt, self.dx coeff_h = dt / (MU_0 * dx) self.Hx[:, :-1, :-1] -= coeff_h * ( (self.Ez[:, 1:, :-1] - self.Ez[:, :-1, :-1]) - (self.Ey[:, :-1, 1:] - self.Ey[:, :-1, :-1])) self.Hy[:-1, :, :-1] -= coeff_h * ( (self.Ex[:-1, :, 1:] - self.Ex[:-1, :, :-1]) - (self.Ez[1:, :, :-1] - self.Ez[:-1, :, :-1])) self.Hz[:-1, :-1, :] -= coeff_h * ( (self.Ey[1:, :-1, :] - self.Ey[:-1, :-1, :]) - (self.Ex[:-1, 1:, :] - self.Ex[:-1, :-1, :])) Ca, Cb = self.Ca, self.Cb self.Ex[:, 1:, 1:] = (Ca[:, 1:, 1:] * self.Ex[:, 1:, 1:] + Cb[:, 1:, 1:] * ((self.Hz[:, 1:, 1:] - self.Hz[:, :-1, 1:]) - (self.Hy[:, 1:, 1:] - self.Hy[:, 1:, :-1]))) self.Ey[1:, :, 1:] = (Ca[1:, :, 1:] * self.Ey[1:, :, 1:] + Cb[1:, :, 1:] * ((self.Hx[1:, :, 1:] - self.Hx[1:, :, :-1]) - (self.Hz[1:, :, 1:] - self.Hz[:-1, :, 1:]))) self.Ez[1:, 1:, :] = (Ca[1:, 1:, :] * self.Ez[1:, 1:, :] + Cb[1:, 1:, :] * ((self.Hy[1:, 1:, :] - self.Hy[:-1, 1:, :]) - (self.Hx[1:, 1:, :] - self.Hx[1:, :-1, :]))) def field_at(self, ix, iy, iz): return float(np.sqrt( self.Ex[ix, iy, iz]**2 + self.Ey[ix, iy, iz]**2 + self.Ez[ix, iy, iz]**2)) def total_energy(self): E_e = 0.5 * EPS_0 * np.sum(self.eps_r * (self.Ex**2 + self.Ey**2 + self.Ez**2)) E_h = 0.5 * MU_0 * np.sum(self.Hx**2 + self.Hy**2 + self.Hz**2) return float(E_e + E_h) # ============================================================================== # WIGNER-SEITZ CELL GEOMETRY DEFINITIONS # ============================================================================== def get_truncated_octahedron(): """BCC Wigner-Seitz cell: truncated octahedron. 14 faces (8 hexagons + 6 squares), 24 vertices, 36 edges. Angular deficit: 30° per vertex (uniform). Vertices at permutations of (0, ±1, ±2) normalized to unit sphere. """ verts = [] # All permutations of (0, ±1, ±2) coords = [0, 1, -1, 2, -2] for perm in [(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)]: for s1 in [1, -1]: for s2 in [1, -1]: v = [0.0, 0.0, 0.0] v[perm[0]] = 0.0 v[perm[1]] = s1 * 1.0 v[perm[2]] = s2 * 2.0 verts.append(v) # Remove duplicates unique = [] for v in verts: is_dup = False for u in unique: if all(abs(v[i] - u[i]) < 0.01 for i in range(3)): is_dup = True break if not is_dup: unique.append(v) verts = np.array(unique, dtype=float) verts /= np.max(np.linalg.norm(verts, axis=1)) # normalize to unit sphere from scipy.spatial import ConvexHull hull = ConvexHull(verts) # Classify vertices by type (all equivalent for truncated octahedron) vertex_types = ["uniform"] * len(verts) return verts, hull.simplices.tolist(), vertex_types def get_rhombic_dodecahedron(): """FCC Wigner-Seitz cell: rhombic dodecahedron. 12 faces (congruent rhombuses), 14 vertices, 24 edges. Two vertex types: - 8 three-valent (cubic positions): deficit = +148.4° - 6 four-valent (octahedral positions): deficit = -77.9° """ # 8 three-valent vertices (cube vertices) verts_3v = [] for sx in [-1, 1]: for sy in [-1, 1]: for sz in [-1, 1]: verts_3v.append([sx, sy, sz]) # 6 four-valent vertices (octahedron vertices) verts_4v = [ [2, 0, 0], [-2, 0, 0], [0, 2, 0], [0, -2, 0], [0, 0, 2], [0, 0, -2], ] verts = np.array(verts_3v + verts_4v, dtype=float) verts /= np.max(np.linalg.norm(verts, axis=1)) # normalize from scipy.spatial import ConvexHull hull = ConvexHull(verts) # Classify: first 8 are 3-valent, last 6 are 4-valent vertex_types = (["3-valent_148deg"] * 8) + (["4-valent_neg78deg"] * 6) return verts, hull.simplices.tolist(), vertex_types def get_hcp_wigner_seitz(): """HCP Wigner-Seitz cell: trapezo-rhombic dodecahedron. 12 faces (6 rectangles + 6 trapezoids), 18 vertices, 28 edges. Built from the dual of the HCP lattice points. Anisotropic — different along c-axis vs basal plane. """ # Ideal c/a = sqrt(8/3) c_over_a = math.sqrt(8.0 / 3.0) h = c_over_a / 2 # half-height # 6 vertices in the top hexagonal ring (z = +h/2) # 6 vertices in the bottom hexagonal ring (z = -h/2) # 6 vertices at the equator (z = 0) but offset verts = [] # Top ring (z = +h) for k in range(6): angle = 2 * math.pi * k / 6 r = 1.0 / math.sqrt(3) verts.append([r * math.cos(angle), r * math.sin(angle), h]) # Bottom ring (z = -h), rotated by 30° for k in range(6): angle = 2 * math.pi * k / 6 + math.pi / 6 r = 1.0 / math.sqrt(3) verts.append([r * math.cos(angle), r * math.sin(angle), -h]) # Equatorial ring (z = 0), at larger radius for k in range(6): angle = 2 * math.pi * k / 6 + math.pi / 6 r = 2.0 / math.sqrt(3) verts.append([r * math.cos(angle), r * math.sin(angle), 0]) verts = np.array(verts, dtype=float) verts /= np.max(np.linalg.norm(verts, axis=1)) from scipy.spatial import ConvexHull hull = ConvexHull(verts) # Classify: top ring, bottom ring, equatorial vertex_types = (["top_axial"] * 6) + (["bottom_axial"] * 6) + (["equatorial"] * 6) return verts, hull.simplices.tolist(), vertex_types def get_sphere_points(n=42): """Sphere control.""" pts = [] for i in range(n): theta = math.acos(1 - 2 * (i + 0.5) / n) phi_a = math.pi * (1 + math.sqrt(5)) * i pts.append([math.sin(theta) * math.cos(phi_a), math.sin(theta) * math.sin(phi_a), math.cos(theta)]) return np.array(pts), None, ["uniform"] * n # Also test the interstitial voids directly def get_regular_tetrahedron(): """FCC tetrahedral void: regular tetrahedron. Deficit = 180°/vertex.""" v = np.array([[1,1,1],[1,-1,-1],[-1,1,-1],[-1,-1,1]], dtype=float) v /= np.linalg.norm(v[0]) from scipy.spatial import ConvexHull hull = ConvexHull(v) return v, hull.simplices.tolist(), ["tet_vertex"] * 4 def get_regular_octahedron(): """FCC octahedral void: regular octahedron. Deficit = 120°/vertex.""" v = np.array([[1,0,0],[-1,0,0],[0,1,0],[0,-1,0],[0,0,1],[0,0,-1]], dtype=float) from scipy.spatial import ConvexHull hull = ConvexHull(v) return v, hull.simplices.tolist(), ["oct_vertex"] * 6 # ============================================================================== # BUILD CAVITY FROM CONVEX HULL # ============================================================================== def build_cavity(sim, verts, faces, R_cells, wall_cells=2, is_sphere=False): """Build PEC cavity from vertex/face geometry.""" N = sim.N center = N // 2 R_inner = R_cells R_outer = R_cells + wall_cells if is_sphere: coords = np.mgrid[0:N, 0:N, 0:N].astype(np.float32) coords[0] -= center; coords[1] -= center; coords[2] -= center r2 = coords[0]**2 + coords[1]**2 + coords[2]**2 inside_inner = r2 <= R_inner**2 inside_outer = r2 <= R_outer**2 wall = inside_outer & ~inside_inner sim.set_conductor(wall) return int(inside_inner.sum()) face_normals = [] face_dist_inner = [] face_dist_outer = [] for f in faces: v0, v1, v2 = verts[f[0]], verts[f[1]], verts[f[2]] normal = np.cross(v1 - v0, v2 - v0) norm = np.linalg.norm(normal) if norm < 1e-10: continue normal /= norm d = np.dot(v0, normal) if d < 0: normal = -normal d = -d face_normals.append(normal) face_dist_inner.append(d * R_inner) face_dist_outer.append(d * R_outer) face_normals = np.array(face_normals) coords = np.mgrid[0:N, 0:N, 0:N].astype(np.float32) coords[0] -= center; coords[1] -= center; coords[2] -= center inside_outer = np.ones((N, N, N), dtype=bool) inside_inner = np.ones((N, N, N), dtype=bool) for i in range(len(face_normals)): n = face_normals[i] sd = n[0] * coords[0] + n[1] * coords[1] + n[2] * coords[2] inside_outer &= (sd <= face_dist_outer[i]) inside_inner &= (sd <= face_dist_inner[i]) wall = inside_outer & ~inside_inner sim.set_conductor(wall) return int(inside_inner.sum()) # ============================================================================== # RUN ONE GEOMETRY # ============================================================================== def run_geometry(name, verts, faces, vertex_types, N=96, steps=2000, is_sphere=False): """Run EM cavity resonance test for one Wigner-Seitz cell.""" print(f"\n {name} ({len(verts)} vertices)") cavity_radius_m = 0.05 dx = 2.2 * cavity_radius_m / N * 2 R_cells = int(cavity_radius_m / dx) freq = C_LIGHT / (4 * cavity_radius_m) sim = EMFDTD3D(N, dx) n_interior = build_cavity(sim, verts, faces, R_cells, is_sphere=is_sphere) print(f" Interior: {n_interior} cells") if n_interior < 10: print(f" SKIP: too few interior cells") return None center = N // 2 omega = 2 * np.pi * freq t_width = 5.0 / freq # Place probes at each vertex (85% radius to stay inside wall) vtx_pos = {} max_r = np.max(np.linalg.norm(verts, axis=1)) for i, v in enumerate(verts): vn = v / max_r * R_cells * 0.85 ix = max(1, min(N-2, int(center + vn[0]))) iy = max(1, min(N-2, int(center + vn[1]))) iz = max(1, min(N-2, int(center + vn[2]))) vtx_pos[f"v{i:02d}"] = (ix, iy, iz) # Equatorial midpoints for baseline eq_pos = {} for k in range(6): angle = 2 * math.pi * k / 6 r = R_cells * 0.4 ix = max(1, min(N-2, int(center + r * math.cos(angle)))) iy = max(1, min(N-2, int(center + r * math.sin(angle)))) eq_pos[f"eq{k}"] = (ix, iy, center) vtx_series = {name: [] for name in vtx_pos} eq_series = {name: [] for name in eq_pos} t_start = time_mod.monotonic() for step in range(steps): t = step * sim.dt envelope = np.exp(-0.5 * ((t - 3 * t_width) / t_width) ** 2) source = envelope * np.sin(omega * t) # 3 sources at 120° in equatorial plane for k in range(3): angle = 2 * math.pi * k / 3 sx = max(1, min(N-2, int(center + R_cells * 0.3 * math.cos(angle)))) sy = max(1, min(N-2, int(center + R_cells * 0.3 * math.sin(angle)))) sim.Ez[sx, sy, center] += source sim.step() if step % 10 == 0: for vname, (vx, vy, vz) in vtx_pos.items(): vtx_series[vname].append(sim.field_at(vx, vy, vz)) for ename, (ex, ey, ez) in eq_pos.items(): eq_series[ename].append(sim.field_at(ex, ey, ez)) if step > 0 and step % (steps // 5) == 0: print(f" {100*step//steps}% E={sim.total_energy():.3e}") elapsed = time_mod.monotonic() - t_start # Analyze vertex RMS (skip first 20% for transient) vtx_rms = {} for vname, series in vtx_series.items(): arr = np.array(series[len(series)//5:]) vtx_rms[vname] = float(np.sqrt(np.mean(arr**2))) if len(arr) > 0 else 0 eq_rms_vals = [] for ename, series in eq_series.items(): arr = np.array(series[len(series)//5:]) eq_rms_vals.append(float(np.sqrt(np.mean(arr**2))) if len(arr) > 0 else 0) max_vtx_rms = max(vtx_rms.values()) if vtx_rms else 0 min_vtx_rms = min(vtx_rms.values()) if vtx_rms else 0 mean_vtx_rms = np.mean(list(vtx_rms.values())) if vtx_rms else 0 std_vtx_rms = np.std(list(vtx_rms.values())) if vtx_rms else 0 mean_eq_rms = np.mean(eq_rms_vals) if eq_rms_vals else 0 concentration = max_vtx_rms / mean_eq_rms if mean_eq_rms > 0 else 0 uniformity = 1.0 - (std_vtx_rms / mean_vtx_rms) if mean_vtx_rms > 0 else 0 differentiation = max_vtx_rms / min_vtx_rms if min_vtx_rms > 0 else float('inf') # Group by vertex type type_groups = {} for i, vtype in enumerate(vertex_types[:len(verts)]): vname = f"v{i:02d}" if vtype not in type_groups: type_groups[vtype] = [] type_groups[vtype].append(vtx_rms.get(vname, 0)) type_summary = {} for vtype, rms_list in type_groups.items(): type_summary[vtype] = { "count": len(rms_list), "mean_rms": float(np.mean(rms_list)), "std_rms": float(np.std(rms_list)), "max_rms": float(np.max(rms_list)), "min_rms": float(np.min(rms_list)), } # Anisotropy check (for HCP): compare axial vs equatorial axial_rms = [vtx_rms[f"v{i:02d}"] for i in range(len(verts)) if vertex_types[i] in ("top_axial", "bottom_axial")] equatorial_vtx_rms = [vtx_rms[f"v{i:02d}"] for i in range(len(verts)) if vertex_types[i] == "equatorial"] anisotropy = (np.mean(axial_rms) / np.mean(equatorial_vtx_rms) if axial_rms and equatorial_vtx_rms and np.mean(equatorial_vtx_rms) > 0 else None) result = { "geometry": name, "n_vertices": len(verts), "n_interior": n_interior, "max_vertex_rms": max_vtx_rms, "min_vertex_rms": min_vtx_rms, "mean_vertex_rms": mean_vtx_rms, "std_vertex_rms": std_vtx_rms, "mean_equatorial_rms": mean_eq_rms, "concentration_ratio": concentration, "uniformity": uniformity, "differentiation_ratio": differentiation, "vertex_type_summary": type_summary, "anisotropy_ratio": anisotropy, "total_energy_final": sim.total_energy(), "elapsed_s": round(elapsed, 1), "per_vertex_rms": {k: round(v, 8) for k, v in vtx_rms.items()}, } print(f" Concentration: {concentration:.1f}x") print(f" Uniformity: {uniformity:.3f} (1.0=perfectly uniform)") print(f" Differentiation: {differentiation:.1f}x (max/min vertex)") if anisotropy is not None: print(f" Anisotropy (axial/equatorial): {anisotropy:.3f}") for vtype, summary in type_summary.items(): print(f" {vtype}: mean={summary['mean_rms']:.3e} " f"±{summary['std_rms']:.3e} (n={summary['count']})") print(f" {elapsed:.1f}s") return result # ============================================================================== # MAIN # ============================================================================== def main(): parser = argparse.ArgumentParser( description="HPC-030: Wigner-Seitz Cell Internal Resonance (EM)") parser.add_argument("--grid", type=int, default=96) parser.add_argument("--steps", type=int, default=2000) args = parser.parse_args() print("=" * 70) print("HPC-030: WIGNER-SEITZ CELL INTERNAL RESONANCE — EM FDTD") print(f" Grid: {args.grid}³, Steps: {args.steps}") print(f" Date: {datetime.now(timezone.utc).isoformat()}") print(f" Engine: EMFDTD3D (Yee algorithm, Maxwell's equations)") print(f" Purpose: Bridge cipher archetypes to internal cavity physics") print("=" * 70) GEOMETRIES = [ ("Sphere_control", *get_sphere_points(42), True), ("BCC_truncated_octahedron", *get_truncated_octahedron(), False), ("FCC_rhombic_dodecahedron", *get_rhombic_dodecahedron(), False), ("HCP_trapezo_rhombic_dodec", *get_hcp_wigner_seitz(), False), ("FCC_tet_void", *get_regular_tetrahedron(), False), ("FCC_oct_void", *get_regular_octahedron(), False), ] results = [] for name, verts, faces, vtx_types, is_sphere in GEOMETRIES: result = run_geometry(name, verts, faces, vtx_types, N=args.grid, steps=args.steps, is_sphere=is_sphere) if result: results.append(result) # Summary print(f"\n{'='*70}") print("WIGNER-SEITZ CELL RESONANCE COMPARISON") print(f"{'='*70}") print(f"{'Geometry':>30} {'Verts':>6} {'Conc.':>7} {'Uniform':>8} " f"{'Diff.':>7} {'Aniso':>7}") print("-" * 70) for r in results: aniso = f"{r['anisotropy_ratio']:.3f}" if r['anisotropy_ratio'] else " N/A" print(f"{r['geometry']:>30} {r['n_vertices']:>6} " f"{r['concentration_ratio']:>6.1f}x {r['uniformity']:>7.3f} " f"{r['differentiation_ratio']:>6.1f}x {aniso:>7}") # Specific cipher predictions print(f"\n{'='*70}") print("CIPHER PREDICTION CHECKS") print(f"{'='*70}") bcc = next((r for r in results if "BCC" in r["geometry"]), None) fcc = next((r for r in results if "FCC_rhombic" in r["geometry"]), None) hcp = next((r for r in results if "HCP" in r["geometry"]), None) if bcc: print(f"\n P1: BCC UNIFORM DISTRIBUTION") print(f" Uniformity: {bcc['uniformity']:.3f} " f"(predicted: high, near 1.0)") print(f" Differentiation: {bcc['differentiation_ratio']:.1f}x " f"(predicted: low, near 1.0)") check = "PASS" if bcc['uniformity'] > 0.5 else "FAIL" print(f" Result: {check}") if fcc: print(f"\n P2: FCC SELECTIVE CONCENTRATION") ts = fcc['vertex_type_summary'] if '3-valent_148deg' in ts and '4-valent_neg78deg' in ts: ratio_3v_4v = (ts['3-valent_148deg']['mean_rms'] / ts['4-valent_neg78deg']['mean_rms'] if ts['4-valent_neg78deg']['mean_rms'] > 0 else float('inf')) print(f" 3-valent (148° deficit) avg: {ts['3-valent_148deg']['mean_rms']:.3e}") print(f" 4-valent (-78° deficit) avg: {ts['4-valent_neg78deg']['mean_rms']:.3e}") print(f" 3-valent / 4-valent ratio: {ratio_3v_4v:.1f}x " f"(predicted: > 3x)") check = "PASS" if ratio_3v_4v > 1.5 else "FAIL" print(f" Result: {check}") if hcp: print(f"\n P3: HCP ANISOTROPY") if hcp['anisotropy_ratio'] is not None: print(f" Axial/equatorial ratio: {hcp['anisotropy_ratio']:.3f} " f"(predicted: ≠ 1.0)") check = "PASS" if abs(hcp['anisotropy_ratio'] - 1.0) > 0.1 else "FAIL" print(f" Result: {check}") if bcc and fcc: print(f"\n P5: FCC vs BCC CONCENTRATION COMPARISON") print(f" FCC concentration: {fcc['concentration_ratio']:.1f}x") print(f" BCC concentration: {bcc['concentration_ratio']:.1f}x") if bcc['concentration_ratio'] > 0: ratio = fcc['concentration_ratio'] / bcc['concentration_ratio'] print(f" FCC/BCC ratio: {ratio:.1f}x " f"(should be > 1 if FCC concentrates more)") # Save timestamp = datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S") out_path = RESULTS_DIR / f"hpc030_wigner_seitz_{timestamp}.json" with open(out_path, "w") as f: json.dump(results, f, indent=2, default=float) print(f"\nResults: {out_path}") if __name__ == "__main__": main()