#!/usr/bin/env python3 """ HPC-029: Real Materials Validation — EM FDTD =============================================== Date: 2026-03-28 Depends on: HPC-026 (PEC baseline), HPC-027 (optimal bicone angle) Validates geometric concentration and compute chain using REAL material properties instead of idealized PEC. Two material scenarios: Scenario A — Microwave PoC (what we're physically building): Walls: aluminum foil (sigma=3.77e7 S/m) Interior: air (eps_r=1.0) Shell: PLA plastic (eps_r=2.7, structural only) Frequency: 1.5 GHz Scenario B — Silicon Photonic Chip (target product): Body: silicon (eps_r=11.68) Cladding: SiO2 (eps_r=3.9) Cavity: air-etched void (eps_r=1.0) Reflectors: gold (sigma=4.1e7 S/m) Frequency: 193 THz (1550nm telecom) Tests: 1. Bicone concentration (35° optimal) — PEC vs aluminum vs copper vs silicon 2. Icosahedron vertex differentiation — PEC vs real materials 3. Compute chain (3-cavity) — PEC vs real materials Critical for: patent credibility. Shows results hold with real physics. Usage: python3 HPC-029_real_materials_validation.py python3 HPC-029_real_materials_validation.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 RESULTS_DIR = Path(__file__).parent / "results_hpc029" RESULTS_DIR.mkdir(exist_ok=True) # ============================================================================== # REAL MATERIAL PROPERTIES # ============================================================================== MATERIALS = { # --- Conductors --- "PEC": {"sigma": 1e7, "eps_r": 1.0, "desc": "Perfect Electric Conductor (idealized)"}, "aluminum": {"sigma": 3.77e7, "eps_r": 1.0, "desc": "Aluminum foil (PoC lining)"}, "copper": {"sigma": 5.96e7, "eps_r": 1.0, "desc": "Copper tape (PoC alternative)"}, "gold": {"sigma": 4.1e7, "eps_r": 1.0, "desc": "Gold contacts (chip-scale)"}, # --- Dielectrics --- "air": {"sigma": 0.0, "eps_r": 1.0, "desc": "Air (cavity interior)"}, "PLA": {"sigma": 1e-14, "eps_r": 2.7, "desc": "PLA plastic (3D print shell)"}, "silicon": {"sigma": 11.7, "eps_r": 11.68, "desc": "Intrinsic silicon (chip body)"}, "SiO2": {"sigma": 1e-14, "eps_r": 3.9, "desc": "Silicon dioxide (cladding)"}, "doped_Si": {"sigma": 1e4, "eps_r": 11.68, "desc": "Heavily doped silicon (reflector)"}, } # Material configurations for testing CONFIGS = { "PEC_baseline": { "wall": "PEC", "interior": "air", "shell": None, "desc": "Idealized PEC walls + vacuum (original tests)", }, "aluminum_PoC": { "wall": "aluminum", "interior": "air", "shell": "PLA", "desc": "Aluminum foil + air + PLA shell (physical PoC)", }, "copper_PoC": { "wall": "copper", "interior": "air", "shell": "PLA", "desc": "Copper tape + air + PLA shell (PoC alternative)", }, "silicon_photonic": { "wall": "doped_Si", "interior": "air", "shell": "SiO2", "desc": "Doped silicon reflectors + air cavity + SiO2 cladding (chip)", }, "gold_photonic": { "wall": "gold", "interior": "air", "shell": "SiO2", "desc": "Gold reflectors + air cavity + SiO2 cladding (chip alt)", }, } # ============================================================================== # EM FDTD ENGINE — Yee Algorithm (Maxwell's Equations) # Same validated engine as HPC-025/026/027/028 # ============================================================================== class EMFDTD3D: """3D EM FDTD — Yee algorithm solving Maxwell's curl equations. dH/dt = -(1/mu_0) curl(E) [Faraday's law] eps * dE/dt + sigma * E = curl(H) [Ampere's law with conduction loss] """ 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_material(self, mask, sigma, eps_r): """Set material properties for masked region.""" self.sigma[mask] = sigma self.eps_r[mask] = eps_r self._compute_coefficients() 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) # ============================================================================== # GEOMETRIES # ============================================================================== def get_sphere_points(n=42): 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) def get_icosahedron(): v = [] v.append([0, 0, 1]); v.append([0, 0, -1]) theta = math.atan(0.5); z_u = math.sin(theta); r_u = math.cos(theta) for k in range(5): a = 2 * math.pi * k / 5 v.append([r_u * math.cos(a), r_u * math.sin(a), z_u]) for k in range(5): a = 2 * math.pi * k / 5 + math.pi / 5 v.append([r_u * math.cos(a), r_u * math.sin(a), -z_u]) v = np.array(v) f = [(0,2,3),(0,3,4),(0,4,5),(0,5,6),(0,6,2), (1,8,7),(1,9,8),(1,10,9),(1,11,10),(1,7,11), (2,7,3),(3,7,8),(3,8,4),(4,8,9),(4,9,5),(5,9,10),(5,10,6),(6,10,11),(6,11,2),(2,11,7)] return v, f def make_bicone(half_angle_deg=35, n_eq=12): ha = math.radians(half_angle_deg) v = [[0, 0, 1.0 / math.cos(ha)], [0, 0, -1.0 / math.cos(ha)]] r_eq = math.tan(ha) for k in range(n_eq): a = 2 * math.pi * k / n_eq v.append([r_eq * math.cos(a), r_eq * math.sin(a), 0]) v = np.array(v, dtype=float) v /= np.max(np.linalg.norm(v, axis=1)) from scipy.spatial import ConvexHull hull = ConvexHull(v) f = hull.simplices.tolist() return v, f # ============================================================================== # BUILD CAVITY WITH REAL MATERIALS # ============================================================================== def build_material_cavity(sim, verts, faces, R_cells, config_name, wall_cells=2, shell_cells=2, is_sphere=False): """Build cavity with real material properties.""" config = CONFIGS[config_name] wall_mat = MATERIALS[config["wall"]] int_mat = MATERIALS[config["interior"]] shell_mat = MATERIALS[config["shell"]] if config["shell"] else None N = sim.N center = N // 2 R_inner = R_cells R_wall = R_cells + wall_cells R_shell = R_wall + shell_cells if shell_mat else R_wall 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_wall = (r2 > R_inner**2) & (r2 <= R_wall**2) inside_shell = (r2 > R_wall**2) & (r2 <= R_shell**2) if shell_mat else np.zeros_like(inside_inner) # Set interior sim.set_material(inside_inner, int_mat["sigma"], int_mat["eps_r"]) # Set wall sim.set_material(inside_wall, wall_mat["sigma"], wall_mat["eps_r"]) # Set shell if shell_mat and np.any(inside_shell): sim.set_material(inside_shell, shell_mat["sigma"], shell_mat["eps_r"]) return int(inside_inner.sum()) # Polyhedral cavity face_normals = [] face_dist = { "inner": [], "wall": [], "shell": [], } 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["wall"].append(d * R_wall) face_dist["shell"].append(d * R_shell) 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 = {} for region in ["inner", "wall", "shell"]: mask = 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] mask &= (sd <= face_dist[region][i]) inside[region] = mask interior_mask = inside["inner"] wall_mask = inside["wall"] & ~inside["inner"] shell_mask = inside["shell"] & ~inside["wall"] if shell_mat else np.zeros_like(interior_mask) # Set materials sim.set_material(interior_mask, int_mat["sigma"], int_mat["eps_r"]) sim.set_material(wall_mask, wall_mat["sigma"], wall_mat["eps_r"]) if shell_mat and np.any(shell_mask): sim.set_material(shell_mask, shell_mat["sigma"], shell_mat["eps_r"]) return int(interior_mask.sum()) # ============================================================================== # TEST 1: BICONE CONCENTRATION WITH REAL MATERIALS # ============================================================================== def test_bicone_concentration(N=64, steps=1500): """Test bicone (35°) concentration with different materials.""" print(f"\n{'='*70}") print("TEST 1: BICONE (35°) CONCENTRATION — REAL MATERIALS") print(f"{'='*70}") verts, faces = make_bicone(35, n_eq=12) 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) results = [] for config_name, config in CONFIGS.items(): print(f"\n --- {config_name}: {config['desc']} ---") sim = EMFDTD3D(N, dx) n_interior = build_material_cavity(sim, verts, faces, R_cells, config_name) print(f" Interior: {n_interior} cells") if n_interior < 10: print(f" SKIP: too few interior cells") continue center = N // 2 omega = 2 * np.pi * freq t_width = 5.0 / freq # Probes 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) eq_pos = {} for k in range(3): angle = 2 * math.pi * k / 3 + math.pi / 6 r = R_cells * 0.5 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) for k in range(3): angle = 2 * math.pi * k / 3 sx = max(1, min(N-2, int(center + R_cells * 0.5 * math.cos(angle)))) sy = max(1, min(N-2, int(center + R_cells * 0.5 * 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 // 4) == 0: print(f" {100*step//steps}% E={sim.total_energy():.3e}") elapsed = time_mod.monotonic() - t_start 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 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 pi_per_cell = max_vtx_rms**2 / max(n_interior, 1) result = { "test": "bicone_35deg", "config": config_name, "materials": config["desc"], "wall_sigma": MATERIALS[config["wall"]]["sigma"], "wall_eps_r": MATERIALS[config["wall"]]["eps_r"], "interior_eps_r": MATERIALS[config["interior"]]["eps_r"], "n_interior": n_interior, "max_vertex_rms": max_vtx_rms, "mean_equatorial_rms": mean_eq_rms, "concentration_ratio": concentration, "pi_per_cell": pi_per_cell, "total_energy_final": sim.total_energy(), "elapsed_s": round(elapsed, 1), } results.append(result) print(f" Concentration: {concentration:.1f}x, " f"PI/cell: {pi_per_cell:.3e}, E_total={sim.total_energy():.3e}, {elapsed:.1f}s") # Summary print(f"\n BICONE CONCENTRATION — MATERIAL COMPARISON:") print(f" {'Config':>20} {'Wall σ':>12} {'Conc.':>8} {'PI/cell':>12} {'E_total':>10}") print(f" {'-'*66}") baseline_conc = results[0]["concentration_ratio"] if results else 1 for r in results: pct = r["concentration_ratio"] / baseline_conc * 100 if baseline_conc > 0 else 0 print(f" {r['config']:>20} {r['wall_sigma']:>11.2e} " f"{r['concentration_ratio']:>7.1f}x {r['pi_per_cell']:>12.3e} " f"{r['total_energy_final']:>9.2e} ({pct:.0f}% of PEC)") return results # ============================================================================== # TEST 2: ICOSAHEDRON VERTEX DIFFERENTIATION WITH REAL MATERIALS # ============================================================================== def test_icosahedron_differentiation(N=64, steps=1500): """Test icosahedron vertex selection with real materials.""" print(f"\n{'='*70}") print("TEST 2: ICOSAHEDRON VERTEX DIFFERENTIATION — REAL MATERIALS") print(f"{'='*70}") verts, faces = get_icosahedron() 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) results = [] for config_name in ["PEC_baseline", "aluminum_PoC", "copper_PoC"]: config = CONFIGS[config_name] print(f"\n --- {config_name}: {config['desc']} ---") sim = EMFDTD3D(N, dx) n_interior = build_material_cavity(sim, verts, faces, R_cells, config_name) print(f" Interior: {n_interior} cells") center = N // 2 omega = 2 * np.pi * freq t_width = 5.0 / freq 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) vtx_series = {name: [] for name in vtx_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) sx = max(1, min(N-2, int(center + R_cells * 0.5))) sim.Ez[sx, center, 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)) if step > 0 and step % (steps // 4) == 0: print(f" {100*step//steps}% E={sim.total_energy():.3e}") elapsed = time_mod.monotonic() - t_start 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 sorted_vtx = sorted(vtx_rms.values(), reverse=True) top_6 = np.mean(sorted_vtx[:6]) if len(sorted_vtx) >= 6 else 0 bottom_6 = np.mean(sorted_vtx[6:]) if len(sorted_vtx) >= 12 else 0 diff_ratio = top_6 / bottom_6 if bottom_6 > 0 else float('inf') n_active = sum(1 for v in sorted_vtx if v > 0.1 * sorted_vtx[0]) result = { "test": "icosahedron_vertex_diff", "config": config_name, "materials": config["desc"], "n_interior": n_interior, "top_6_avg_rms": top_6, "bottom_6_avg_rms": bottom_6, "differentiation_ratio": diff_ratio, "active_vertices": n_active, "total_vertices": 12, "elapsed_s": round(elapsed, 1), } results.append(result) print(f" Top 6 avg: {top_6:.3e}, Bottom 6 avg: {bottom_6:.3e}") print(f" Differentiation: {diff_ratio:.1f}x, Active: {n_active}/12, {elapsed:.1f}s") return results # ============================================================================== # TEST 3: SPHERE CONTROL WITH REAL MATERIALS # ============================================================================== def test_sphere_control(N=64, steps=1500): """Sphere control — should show NO concentration regardless of material.""" print(f"\n{'='*70}") print("TEST 3: SPHERE CONTROL — REAL MATERIALS") print(f"{'='*70}") verts = get_sphere_points(42) 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) results = [] for config_name in ["PEC_baseline", "aluminum_PoC", "copper_PoC"]: config = CONFIGS[config_name] print(f"\n --- {config_name} ---") sim = EMFDTD3D(N, dx) n_interior = build_material_cavity(sim, verts, None, R_cells, config_name, is_sphere=True) center = N // 2 omega = 2 * np.pi * freq t_width = 5.0 / freq vtx_pos = {} for i, v in enumerate(verts[:12]): vn = v * 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) eq_pos = {} for k in range(3): angle = 2 * math.pi * k / 3 + math.pi / 6 r = R_cells * 0.5 eq_pos[f"eq{k}"] = ( max(1, min(N-2, int(center + r * math.cos(angle)))), max(1, min(N-2, int(center + r * math.sin(angle)))), center) vtx_series = {n: [] for n in vtx_pos} eq_series = {n: [] for n in eq_pos} 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) for k in range(3): angle = 2 * math.pi * k / 3 sx = max(1, min(N-2, int(center + R_cells * 0.5 * math.cos(angle)))) sy = max(1, min(N-2, int(center + R_cells * 0.5 * math.sin(angle)))) sim.Ez[sx, sy, center] += source sim.step() if step % 10 == 0: for vn, (vx, vy, vz) in vtx_pos.items(): vtx_series[vn].append(sim.field_at(vx, vy, vz)) for en, (ex, ey, ez) in eq_pos.items(): eq_series[en].append(sim.field_at(ex, ey, ez)) vtx_rms = [float(np.sqrt(np.mean(np.array(s[len(s)//5:])**2))) for s in vtx_series.values()] eq_rms = [float(np.sqrt(np.mean(np.array(s[len(s)//5:])**2))) for s in eq_series.values()] conc = np.mean(vtx_rms) / np.mean(eq_rms) if np.mean(eq_rms) > 0 else 0 results.append({ "test": "sphere_control", "config": config_name, "concentration_ratio": conc, }) print(f" Sphere concentration: {conc:.2f}x (should be ~1.0)") return results # ============================================================================== # MAIN # ============================================================================== def main(): parser = argparse.ArgumentParser( description="HPC-029: Real Materials Validation (EM)") parser.add_argument("--grid", type=int, default=64) parser.add_argument("--steps", type=int, default=1500) args = parser.parse_args() print("=" * 70) print("HPC-029: REAL MATERIALS VALIDATION — 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: Validate concentration with real material properties") print(f" Materials tested: {list(CONFIGS.keys())}") print("=" * 70) all_results = {} # Test 1: Bicone concentration all_results["bicone_concentration"] = test_bicone_concentration( N=args.grid, steps=args.steps) # Test 2: Icosahedron differentiation all_results["icosahedron_differentiation"] = test_icosahedron_differentiation( N=args.grid, steps=args.steps) # Test 3: Sphere control all_results["sphere_control"] = test_sphere_control( N=args.grid, steps=args.steps) # Final comparison print(f"\n{'='*70}") print("FINAL COMPARISON — REAL MATERIALS vs PEC") print(f"{'='*70}") bicone = all_results["bicone_concentration"] ico = all_results["icosahedron_differentiation"] sphere = all_results["sphere_control"] print("\n BICONE (35°) CONCENTRATION:") for r in bicone: print(f" {r['config']:>20}: {r['concentration_ratio']:.1f}x") print("\n ICOSAHEDRON VERTEX DIFFERENTIATION:") for r in ico: print(f" {r['config']:>20}: {r['differentiation_ratio']:.1f}x " f"({r['active_vertices']}/{r['total_vertices']} active)") print("\n SPHERE CONTROL:") for r in sphere: print(f" {r['config']:>20}: {r['concentration_ratio']:.2f}x " f"(should be ~1.0)") # Save timestamp = datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S") out_path = RESULTS_DIR / f"hpc029_real_materials_{timestamp}.json" save_data = { "metadata": { "test": "HPC-029", "engine": "EMFDTD3D (Yee FDTD, Maxwell's equations)", "grid": args.grid, "steps": args.steps, "date": datetime.now(timezone.utc).isoformat(), "materials": {k: v for k, v in MATERIALS.items()}, "configs": {k: v for k, v in CONFIGS.items()}, }, "results": all_results, } with open(out_path, "w") as f: json.dump(save_data, f, indent=2, default=float) print(f"\nResults: {out_path}") if __name__ == "__main__": main()