#!/usr/bin/env python3 """ HPC-028: Frequency Selectivity & Q-Factor — EM FDTD ===================================================== Date: 2026-03-28 Depends on: HPC-026 (bicone 132x at resonance), HPC-027 (optimal angle) Drives a bicone cavity at multiple frequencies (0.5x to 2.0x the resonant frequency) to measure how concentration varies with detuning. This gives: - Q-factor of the geometric concentrator - Bandwidth (how narrow is the concentration peak?) - Off-resonance rejection ratio Critical for: Geometric Photovoltaic Generator patent — determines whether the cavity acts as a narrowband concentrator (enabling frequency-tiling) or broadband (single concentrator sufficient). Also runs the same sweep for Sphere (control) and Tetrahedron (comparison). Usage: python3 HPC-028_frequency_selectivity.py python3 HPC-028_frequency_selectivity.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_hpc028" RESULTS_DIR.mkdir(exist_ok=True) # ============================================================================== # EM FDTD ENGINE — Yee Algorithm (Maxwell's Equations) # Identical to HPC-025/026/027 validated engine # ============================================================================== class EMFDTD3D: """3D EM FDTD solver using the standard Yee algorithm. Solves 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 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_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_tetrahedron(): v = np.array([[1,1,1],[1,-1,-1],[-1,1,-1],[-1,-1,1]], dtype=float) v /= np.linalg.norm(v[0]) f = [(0,1,2),(0,2,3),(0,3,1),(1,3,2)] return v, f def make_bicone(half_angle_deg=30, n_eq=12): ha = math.radians(half_angle_deg) apex_z = 1.0 / math.cos(ha) if ha < math.pi / 2 else 1.0 v = [[0, 0, apex_z], [0, 0, -apex_z]] r_eq = math.tan(ha) if ha < math.pi / 2 else 10.0 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 EM CAVITY # ============================================================================== def build_em_cavity(sim, verts, faces, R_cells, wall_cells=2, is_sphere=False): 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 FREQUENCY # ============================================================================== def run_single(geometry_name, verts, faces, freq_multiplier, N=64, steps=1500, is_sphere=False): """Run EM sim at a specific frequency and measure vertex concentration.""" cavity_radius_m = 0.05 dx = 2.2 * cavity_radius_m / N * 2 R_cells = int(cavity_radius_m / dx) base_freq = C_LIGHT / (4 * cavity_radius_m) freq = base_freq * freq_multiplier sim = EMFDTD3D(N, dx) n_interior = build_em_cavity(sim, verts, faces, R_cells, is_sphere=is_sphere) if n_interior < 10: return None center = N // 2 omega = 2 * np.pi * freq t_width = 5.0 / base_freq # envelope width based on base freq # Vertex 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)) 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) return { "geometry": geometry_name, "freq_multiplier": freq_multiplier, "freq_hz": freq, "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, "elapsed_s": round(elapsed, 1), } # ============================================================================== # MAIN # ============================================================================== def main(): parser = argparse.ArgumentParser( description="HPC-028: Frequency Selectivity & Q-Factor (EM)") parser.add_argument("--grid", type=int, default=64) parser.add_argument("--steps", type=int, default=1500) args = parser.parse_args() # Frequency multipliers: 0.3x to 2.5x in fine steps around 1.0 FREQ_MULTS = [0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0, 1.05, 1.1, 1.15, 1.2, 1.3, 1.4, 1.5, 1.7, 2.0, 2.5] GEOMETRIES = { "Bicone_30deg": {"fn": lambda: make_bicone(30), "is_sphere": False}, "Tetrahedron": {"fn": lambda: get_tetrahedron(), "is_sphere": False}, "Sphere": {"fn": lambda: (get_sphere_points(42), None), "is_sphere": True}, } print("=" * 70) print("HPC-028: FREQUENCY SELECTIVITY & Q-FACTOR — EM FDTD") print(f" Grid: {args.grid}³, Steps: {args.steps}") print(f" Geometries: {list(GEOMETRIES.keys())}") print(f" Freq multipliers: {FREQ_MULTS}") print(f" Date: {datetime.now(timezone.utc).isoformat()}") print(f" Engine: EMFDTD3D (Yee algorithm, Maxwell's equations)") print(f" Purpose: Measure Q-factor and bandwidth for PV patent") print("=" * 70) all_results = {} for geo_name, geo_spec in GEOMETRIES.items(): print(f"\n{'='*50}") print(f" GEOMETRY: {geo_name}") print(f"{'='*50}") verts, faces = geo_spec["fn"]() is_sphere = geo_spec["is_sphere"] geo_results = [] for mult in FREQ_MULTS: print(f"\n freq={mult:.2f}x", end="", flush=True) result = run_single(geo_name, verts, faces, mult, N=args.grid, steps=args.steps, is_sphere=is_sphere) if result: geo_results.append(result) print(f" → conc={result['concentration_ratio']:.1f}x " f"PI={result['pi_per_cell']:.3e}", end="") else: print(f" → SKIP", end="") all_results[geo_name] = geo_results print() # Summary for this geometry if geo_results: print(f"\n {geo_name} FREQUENCY RESPONSE:") print(f" {'Freq×':>8} {'Conc.':>8} {'PI/cell':>12}") print(f" {'-'*32}") peak_result = max(geo_results, key=lambda r: r["concentration_ratio"]) for r in geo_results: marker = " ← PEAK" if r == peak_result else "" print(f" {r['freq_multiplier']:>7.2f}x " f"{r['concentration_ratio']:>7.1f}x " f"{r['pi_per_cell']:>12.3e}{marker}") # Q-factor estimate from -3dB bandwidth peak_conc = peak_result["concentration_ratio"] half_power = peak_conc / math.sqrt(2) in_band = [r for r in geo_results if r["concentration_ratio"] >= half_power] if len(in_band) >= 2: bw = max(r["freq_multiplier"] for r in in_band) - \ min(r["freq_multiplier"] for r in in_band) q_est = peak_result["freq_multiplier"] / bw if bw > 0 else float('inf') print(f"\n -3dB bandwidth: {bw:.2f}x ({bw*100:.0f}% of fundamental)") print(f" Estimated Q-factor: {q_est:.1f}") print(f" Peak at: {peak_result['freq_multiplier']:.2f}x fundamental") else: print(f"\n Q-factor: very high (single-point peak)") # Cross-geometry comparison at resonance print(f"\n{'='*70}") print("CROSS-GEOMETRY COMPARISON AT EACH FREQUENCY") print(f"{'='*70}") print(f"{'Freq×':>8}", end="") for geo_name in GEOMETRIES: print(f" {geo_name:>14}", end="") print() print("-" * 50) for mult in FREQ_MULTS: print(f"{mult:>7.2f}x", end="") for geo_name in GEOMETRIES: gr = all_results.get(geo_name, []) match = [r for r in gr if abs(r["freq_multiplier"] - mult) < 0.001] if match: print(f" {match[0]['concentration_ratio']:>13.1f}x", end="") else: print(f" {'---':>14}", end="") print() # Save all results timestamp = datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S") out_path = RESULTS_DIR / f"hpc028_freq_selectivity_{timestamp}.json" save_data = { "metadata": { "test": "HPC-028", "engine": "EMFDTD3D (Yee FDTD, Maxwell's equations)", "grid": args.grid, "steps": args.steps, "date": datetime.now(timezone.utc).isoformat(), }, "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()