#!/usr/bin/env python3 """ HPC-027: Bicone Apex Angle Sweep — EM FDTD ============================================= Date: 2026-03-28 Depends on: HPC-026 (bicone 30° gave 132x EM concentration) Sweeps bicone half-angles from 10° to 60° to find the optimal apex geometry for photovoltaic concentration. Uses the SAME EM engine (Yee FDTD, Maxwell's equations) validated in HPC-025/026. Critical for: Geometric Photovoltaic Generator provisional patent. The question: is 30° optimal, or is there a sharper/wider angle that concentrates EM even more? Usage: python3 HPC-027_bicone_apex_sweep.py python3 HPC-027_bicone_apex_sweep.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_hpc027" RESULTS_DIR.mkdir(exist_ok=True) # ============================================================================== # EM FDTD ENGINE — Yee Algorithm (Maxwell's Equations) # Identical to HPC-025/026 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] Fields: Ex, Ey, Ez, Hx, Hy, Hz (6 vector components on staggered grid) Boundaries: PEC via high-conductivity regions (sigma = 1e7 S/m) """ 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) # Faraday: dH/dt = -(1/mu_0) curl(E) 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, :])) # Ampere: eps * dE/dt + sigma * E = curl(H) 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) # ============================================================================== # BICONE GEOMETRY # ============================================================================== def make_bicone(half_angle_deg, n_eq=12): """Create bicone with given half-angle. Returns (vertices, faces).""" 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 def angular_deficit_per_vertex(half_angle_deg): """Compute angular deficit at apex of bicone.""" ha = math.radians(half_angle_deg) # Cone apex deficit: 360 - (n_eq * face_angle_at_apex) # For a bicone with n_eq equatorial vertices, each face subtends # an angle at the apex. Total solid angle deficit increases with sharpness. # Approximation: deficit ≈ 360 - 360*sin(ha) return 360 * (1 - math.sin(ha)) 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) # ============================================================================== # 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 ANGLE # ============================================================================== def run_angle(half_angle, N=64, steps=1500, is_sphere=False): """Run EM concentration test for a bicone at given half-angle.""" deficit = angular_deficit_per_vertex(half_angle) if not is_sphere else 0 name = f"Bicone_{half_angle}deg" if not is_sphere else "Sphere_control" print(f"\n {name} (half-angle={half_angle}°, deficit≈{deficit:.0f}°)") if is_sphere: verts = get_sphere_points(42) faces = None else: verts, faces = make_bicone(half_angle, 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) sim = EMFDTD3D(N, dx) n_interior = build_em_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 (cavity too sharp for grid)") return None center = N // 2 omega = 2 * np.pi * freq t_width = 5.0 / freq # Vertex positions (apex and equatorial 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) # Equatorial midpoints (baseline) 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 // 5) == 0: print(f" {100*step//steps}% E={sim.total_energy():.3e}") elapsed = time_mod.monotonic() - t_start # Analyze 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 best_vtx = max(vtx_rms, key=vtx_rms.get) if vtx_rms else "none" 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) # Identify apex vs equatorial vertex concentration apex_rms = [vtx_rms.get("v00", 0), vtx_rms.get("v01", 0)] eq_vtx_rms = [v for k, v in vtx_rms.items() if k not in ("v00", "v01")] apex_avg = np.mean(apex_rms) if apex_rms else 0 eq_vtx_avg = np.mean(eq_vtx_rms) if eq_vtx_rms else 0 apex_to_eq_ratio = apex_avg / eq_vtx_avg if eq_vtx_avg > 0 else 0 result = { "half_angle_deg": half_angle, "deficit_approx_deg": round(deficit, 1), "n_vertices": len(verts), "n_interior": n_interior, "max_vertex_rms": max_vtx_rms, "best_vertex": best_vtx, "mean_equatorial_rms": mean_eq_rms, "concentration_ratio": concentration, "pi_per_cell": pi_per_cell, "apex_avg_rms": apex_avg, "equatorial_vertex_avg_rms": eq_vtx_avg, "apex_to_equatorial_ratio": apex_to_eq_ratio, "elapsed_s": round(elapsed, 1), } print(f" Concentration: {concentration:.1f}x, " f"Apex/Eq: {apex_to_eq_ratio:.1f}x, " f"PI/cell: {pi_per_cell:.3e}, {elapsed:.1f}s") return result # ============================================================================== # MAIN # ============================================================================== def main(): parser = argparse.ArgumentParser(description="HPC-027: Bicone Apex Angle Sweep (EM)") parser.add_argument("--grid", type=int, default=64, help="Grid size N (default 64 = 64^3)") parser.add_argument("--steps", type=int, default=1500, help="Simulation steps (default 1500)") args = parser.parse_args() HALF_ANGLES = [10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60] print("=" * 70) print("HPC-027: BICONE APEX ANGLE SWEEP — EM FDTD") print(f" Grid: {args.grid}³, Steps: {args.steps}") print(f" Angles: {HALF_ANGLES}") print(f" Date: {datetime.now(timezone.utc).isoformat()}") print(f" Engine: EMFDTD3D (Yee algorithm, Maxwell's equations)") print(f" Fields: Ex, Ey, Ez, Hx, Hy, Hz (6-component vector)") print(f" Purpose: Find optimal bicone concentrator geometry for PV patent") print("=" * 70) results = [] # Sphere control print("\n--- SPHERE CONTROL ---") sphere_result = run_angle(0, N=args.grid, steps=args.steps, is_sphere=True) if sphere_result: sphere_result["half_angle_deg"] = 0 sphere_result["deficit_approx_deg"] = 0 results.append(sphere_result) # Bicone sweep print("\n--- BICONE ANGLE SWEEP ---") for angle in HALF_ANGLES: result = run_angle(angle, N=args.grid, steps=args.steps) if result: results.append(result) # Summary print(f"\n{'='*70}") print("BICONE APEX ANGLE SWEEP RESULTS") print(f"{'='*70}") print(f"{'Half-Angle':>12} {'Deficit':>8} {'Interior':>9} {'Conc.':>8} " f"{'Apex/Eq':>8} {'PI/cell':>12}") print("-" * 62) sphere_pi = results[0]["pi_per_cell"] if results else 1 for r in results: angle_str = f"{r['half_angle_deg']}°" if r['half_angle_deg'] > 0 else "Sphere" ratio_to_sphere = r["pi_per_cell"] / sphere_pi if sphere_pi > 0 else 0 print(f"{angle_str:>12} {r['deficit_approx_deg']:>7.0f}° {r['n_interior']:>9} " f"{r['concentration_ratio']:>7.1f}x " f"{r['apex_to_equatorial_ratio']:>7.1f}x " f"{r['pi_per_cell']:>12.3e} ({ratio_to_sphere:.0f}x sphere)") # Find optimal if len(results) > 1: best = max(results[1:], key=lambda r: r["concentration_ratio"]) print(f"\n OPTIMAL: {best['half_angle_deg']}° half-angle → " f"{best['concentration_ratio']:.1f}x concentration, " f"{best['pi_per_cell'] / sphere_pi:.0f}x sphere") # Save timestamp = datetime.now(timezone.utc).strftime("%Y%m%d_%H%M%S") out_path = RESULTS_DIR / f"hpc027_apex_sweep_{timestamp}.json" with open(out_path, "w") as f: json.dump(results, f, indent=2) print(f"\nResults: {out_path}") if __name__ == "__main__": main()