""" v3 waveguide-facet <-> fiber-tip 2PP lens relay: mode solve, Gaussian/ABCD lens design, exact ray-trace verification + misalignment sensitivity, STL export. Pipeline (see README.md / LENS_DESIGN.md for the full writeup): 1. Solve the Si3N4/SiO2/air strip-waveguide fundamental mode (effective-index method: 1D asymmetric slab in z, then 1D symmetric slab in y using the vertical neff as the effective core index) to get the launch mode's 1/e^2 mode field diameter (MFD). 2. Reuse this repo's 630HP fiber numbers (dbt_tir_sil_fiber_coupler_v1 / v2_lens_codesign) as the target: MFD_fiber = 5.0 um, NA = 0.12. 3. Design a two-lens relay (IP-S resin, n=1.51 @ 780 nm): Lens 1 -- printed directly on the waveguide facet, single aspheric resin/air surface, collimates the (idealized point-source) waveguide mode. Lens 2 -- printed on/near the fiber tip, single aspheric air/resin surface, focuses the collimated relay beam down onto the fiber core plane. Vertex radii of curvature are solved from Gaussian-beam (ABCD / complex q-parameter) propagation so the mode WAIST and WAVEFRONT match the fiber's fundamental mode at the target plane; conic constants are fixed at the exact (aplanatic, zero on-axis spherical aberration) collimating value for each interface. 4. Exact sequential ray trace (Snell's law, numpy) of the same surfaces: - verifies on-axis collimation/refocusing for a fan of ray angles, - sweeps lateral misalignment, axial (gap) defocus, and tilt to produce a spot-size/efficiency-vs-perturbation figure. 5. STL export of both lens solids (+ a support-ring stub for the taller, steeper Lens 1 print) for 2PP printing. Run: python fiber_lens_raytrace.py """ from __future__ import annotations import datetime as dt import json from dataclasses import asdict, dataclass from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np from scipy.optimize import brentq, minimize try: from stl import mesh as stlmesh except ImportError as exc: # pragma: no cover raise RuntimeError("numpy-stl (`pip install numpy-stl`) is required for STL export") from exc HERE = Path(__file__).resolve().parent STL_DIR = HERE / "stl" FIG_DIR = HERE / "figures" DATA_DIR = HERE / "data" REPORT_DIR = HERE / "reports" for d in (STL_DIR, FIG_DIR, DATA_DIR, REPORT_DIR): d.mkdir(parents=True, exist_ok=True) MATERIAL_COLORS = { "air": "#ffffff", "sio2": "#e2e2e6", "sin": "#1f3b73", "resist": "#ee8c44", "core": "#cc2d2d", "clad": "#f4d6d6", "ray": "#222222", "beam": "#cc2d2d", "gridline": "#909090", } plt.rcParams.update({ "figure.dpi": 110, "savefig.dpi": 150, "figure.facecolor": "white", "savefig.facecolor": "white", "savefig.bbox": "tight", "font.size": 10, "axes.titlesize": 11, "axes.titleweight": "semibold", "axes.spines.top": False, "axes.spines.right": False, "axes.grid": True, "grid.alpha": 0.28, "grid.linewidth": 0.5, "legend.frameon": False, }) # ============================================================================ # 0. Design constants # ============================================================================ LAM0 = 0.780 # um, DBT:anthracene ZPL design wavelength (780 nm), matches # dbt_tir_sil_fiber_coupler_v1 / single_molecule_device_v1 / v2_lens_codesign # --- waveguide (Si3N4 core / SiO2 under-clad / air top+side-clad) --- N_SIN = 2.00 N_SIO2 = 1.45 N_AIR = 1.00 WG_THICKNESS = 0.300 # um, from device_architectures/v3 wafer stack (fixed by fab) WG_WIDTH = 0.400 # um, chosen -- see justification in main(): matches the width # already used in nanophotonic_devices/emitter_coupling/ # anthracene_slab_wg_coupling_v1 (400x300 nm Si3N4 strip), whose # FDTD mode solve reports fundamental TE n_eff = 1.6720 @ 780 nm, # which we use below as an independent cross-check of the EIM here. FDTD_NEFF_REFERENCE = 1.6720 # from anthracene_slab_wg_coupling_v1 metrics.json (TE0, mode_index=0) # --- fiber: 630HP, reused as-is from dbt_tir_sil_fiber_coupler_v1 / v2_lens_codesign --- FIBER_MFD = 5.0 # um @ 780 nm FIBER_W0 = FIBER_MFD / 2.0 FIBER_NA_DATASHEET = 0.12 # 1%-power datasheet NA (larger than the Gaussian-mode NA below) FIBER_N_CORE = 1.462 FIBER_N_CLAD = 1.457 # --- 2PP resin lens material --- N_LENS = 1.51 # Nanoscribe IP-S @ ~780 nm, same value used elsewhere in this repo # (v2_lens_codesign/shared/params.py: N_IPS) # --- relay geometry (fabrication/assembly design choices, NOT lens-shape free params) --- H1 = 6.0 # um, Lens-1 print height above the waveguide facet GAP = 15.0 # um, free-space gap between the two lenses across the cleave + # v-groove assembly tolerance (assumption -- flag for the GDS/mechanical # design to confirm against the actual cleave-to-v-groove distance) # Lens 2's standoff to the fiber facet is NOT a free/bounded design knob -- it is solved # (see design_lens2 / waist_after_lens2) as the physical distance to the beam's own # natural waist, which for this low-NA relay comes out to several tens of um. # ============================================================================ # 1. Waveguide mode: effective-index method (1D asymmetric slab x2) # ============================================================================ def slab_te0_neff(lam: float, d: float, n_core: float, n_top: float, n_bot: float) -> float: """Fundamental (m=0) TE effective index of a 3-layer asymmetric slab.""" k0 = 2 * np.pi / lam n_clad_max = max(n_top, n_bot) def eq(neff): kappa = k0 * np.sqrt(n_core**2 - neff**2) g_top = k0 * np.sqrt(max(neff**2 - n_top**2, 1e-12)) g_bot = k0 * np.sqrt(max(neff**2 - n_bot**2, 1e-12)) return kappa * d - (np.arctan(g_top / kappa) + np.arctan(g_bot / kappa)) lo, hi = n_clad_max + 1e-6, n_core - 1e-6 scan = np.linspace(hi, lo, 4000) vals = [eq(n) for n in scan] for i in range(len(scan) - 1): if np.sign(vals[i]) != np.sign(vals[i + 1]): return brentq(eq, scan[i + 1], scan[i]) raise RuntimeError("no guided TE0 mode found for this slab") def slab_mode_profile(lam, d, n_core, n_top, n_bot, neff, z): """TE field profile (Okamoto asymmetric-slab form), z=0 at bottom, z=d at top.""" k0 = 2 * np.pi / lam kappa = k0 * np.sqrt(n_core**2 - neff**2) g_top = k0 * np.sqrt(neff**2 - n_top**2) g_bot = k0 * np.sqrt(neff**2 - n_bot**2) phi = np.arctan(g_top / kappa) E = np.zeros_like(z) core, top, bot = (z >= 0) & (z <= d), z > d, z < 0 E[core] = np.cos(kappa * z[core] - kappa * d + phi) E[top] = np.cos(phi) * np.exp(-g_top * (z[top] - d)) E[bot] = np.cos(-kappa * d + phi) * np.exp(g_bot * z[bot]) return E def mfd_1e2(z: np.ndarray, E: np.ndarray) -> float: """1/e^2 intensity full width, by direct threshold crossing (no Gaussian assumption).""" I = E**2 I /= I.max() i0 = int(np.argmax(I)) def cross(i0, direction): i = i0 while 0 <= i + direction < len(I): j = i + direction if I[j] <= 1 / np.e**2: frac = (I[i] - 1 / np.e**2) / (I[i] - I[j]) return z[i] + frac * (z[j] - z[i]) i = j return z[-1] if direction > 0 else z[0] return cross(i0, +1) - cross(i0, -1) def solve_waveguide_mode() -> dict: z = np.linspace(-2, 2, 20000) neff_z = slab_te0_neff(LAM0, WG_THICKNESS, N_SIN, N_AIR, N_SIO2) Ez = slab_mode_profile(LAM0, WG_THICKNESS, N_SIN, N_AIR, N_SIO2, neff_z, z) mfd_z = mfd_1e2(z, Ez) y = np.linspace(-2, 2, 20000) neff_2d = slab_te0_neff(LAM0, WG_WIDTH, neff_z, N_AIR, N_AIR) Ey = slab_mode_profile(LAM0, WG_WIDTH, neff_z, N_AIR, N_AIR, neff_2d, y) mfd_y = mfd_1e2(y, Ey) mfd_circ = np.sqrt(mfd_y * mfd_z) return dict( neff_vertical_slab=neff_z, neff_2d_eim=neff_2d, fdtd_reference_neff=FDTD_NEFF_REFERENCE, eim_vs_fdtd_pct_diff=100 * (neff_2d - FDTD_NEFF_REFERENCE) / FDTD_NEFF_REFERENCE, mfd_y=mfd_y, mfd_z=mfd_z, mfd_circularized=mfd_circ, w0_wg=mfd_circ / 2, z_profile=(z, Ez), y_profile=(y, Ey), ) # ============================================================================ # 2. Gaussian-beam (complex q-parameter) ABCD propagation -- validated against # (a) flat-interface index scaling q2 = q1*n2/n1, (b) thin-lens EFL formula. # ============================================================================ def q_from_w(w, n, R=np.inf): if np.isinf(R): return 1j * np.pi * w**2 * n / LAM0 inv_q = 1.0 / R - 1j * LAM0 / (np.pi * n * w**2) return 1.0 / inv_q def w_and_R_from_q(q, n): inv_q = 1.0 / q R = np.inf if inv_q.real == 0 else 1.0 / inv_q.real w = np.sqrt(-LAM0 / (np.pi * n * inv_q.imag)) return w, R def propagate(q, d): return q + d def refract_q(q1, n1, n2, R): """Spherical/paraxial interface, radius R (>0: center of curvature on the n2 side).""" inv_q2 = (n1 / q1 + (n1 - n2) / R) / n2 return 1.0 / inv_q2 def gaussian_overlap(w1, R1, w2, R2, n_common=1.0): """Standard two-Gaussian-mode overlap (waist-mismatch + wavefront-curvature-mismatch), e.g. Neumann, 'Single-Mode Fibers', eq. for Gaussian-beam-to-fiber coupling. inv_R is 0 for a flat wavefront (R=inf) -- a mismatch between one flat and one curved wavefront IS penalized (only R1==R2==inf gives zero curvature penalty).""" k = 2 * np.pi * n_common / LAM0 inv_R1 = 0.0 if not np.isfinite(R1) else 1.0 / R1 inv_R2 = 0.0 if not np.isfinite(R2) else 1.0 / R2 term_curv = (k * w1 * w2 / 2.0 * (inv_R1 - inv_R2)) ** 2 denom = (w1 / w2 + w2 / w1) ** 2 + term_curv return 4.0 / denom def fresnel_T_normal(n1, n2): R = ((n1 - n2) / (n1 + n2)) ** 2 return 1.0 - R # ============================================================================ # 3. Lens design: solve R1 (paraxial collimation), then optimize R2 + standoff2 # ============================================================================ def solve_R1(w0_wg, h1): """R1_sag (exact-conic 'sag equation' convention, R>0, matching lens_asphere.py / the sag()/height() functions below) such that the WG-mode waist (in resin, at the facet) collimates after propagating h1 in resin and refracting resin->air at Lens-1's vertex. Note on sign conventions: refract_q() below uses the standard paraxial convention R_optical>0 <=> center of curvature on the transmission (n2) side. Lens 1's physical shape is a dome whose apex is the FARTHEST point from the source (apex_z=h1, height decreasing outward, i.e. height(r)=h1-sag(r)); for that orientation the center of curvature sits on the *source* (n1) side, so R_optical = -R_sag. We solve in R_optical space (bracketing negative R, verified numerically) and return R_sag.""" q_source = q_from_w(w0_wg, N_LENS, R=np.inf) def residual(R_optical): q_at_l1 = propagate(q_source, h1) q_after = refract_q(q_at_l1, N_LENS, N_AIR, R_optical) return (1.0 / q_after).real # zero when output wavefront is flat (collimated) lo, hi = -500.0, -0.02 R1_optical = brentq(residual, lo, hi) return -R1_optical # R_sag, positive def propagate_system(R1, h1, R2, gap, standoff2, w0_wg): """R1, R2 are given in the 'sag equation' convention (R>0, matching sag()/height() below). Lens 1's dome apex faces AWAY from its source (height=h1-sag(r)) so its paraxial-optical R has the opposite sign (R_optical=-R1); Lens 2's dome apex faces TOWARD its source/the gap (height=apex_z+sag(r)) so its paraxial-optical R has the same sign (R_optical=+R2). Both signs verified numerically against the standard single-surface focal-length formula -- see scratch validation in the design notes.""" q = q_from_w(w0_wg, N_LENS, R=np.inf) q = propagate(q, h1) q = refract_q(q, N_LENS, N_AIR, -R1) q = propagate(q, gap) q = refract_q(q, N_AIR, N_LENS, R2) q = propagate(q, standoff2) w, R = w_and_R_from_q(q, N_LENS) return w, R, q def waist_after_lens2(R1, h1, gap, R2, w0_wg): """For a candidate R2, propagate to Lens 2's vertex, refract, then find the distance (standoff2) to the beam's OWN natural waist beyond the vertex (where the wavefront is automatically flat, R=inf) and the waist size there. This decouples the two design targets (flat wavefront <-> waist size) instead of a joint 2D search that can stall at box bounds.""" q = q_from_w(w0_wg, N_LENS, R=np.inf) q = propagate(q, h1) q = refract_q(q, N_LENS, N_AIR, -R1) q = propagate(q, gap) q = refract_q(q, N_AIR, N_LENS, R2) standoff2 = -q.real # distance to the natural waist (Re(q)=0) from the L2 vertex q_waist = q + standoff2 w0, _ = w_and_R_from_q(q_waist, N_LENS) return standoff2, w0 def design_lens2(R1, h1, gap, w0_wg): """Solve for R2 such that Lens 2's own natural output waist (which by construction forms with a flat wavefront, R=inf -- exactly the fiber's own waist-plane condition) has size w0_fiber. This is the 'mode-overlap-integral' design step the ray-trace-only route can't rigorously provide, reduced to a 1D root-find (R2) instead of a fragile 2D joint optimization.""" def residual(R2): _, w0 = waist_after_lens2(R1, h1, gap, R2, w0_wg) return w0 - FIBER_W0 lo, hi = 0.2, 2000.0 r_lo, r_hi = residual(lo), residual(hi) if np.sign(r_lo) == np.sign(r_hi): # fall back to a coarse scan to find a bracket scan = np.geomspace(lo, hi, 400) vals = [residual(r) for r in scan] R2 = None for i in range(len(scan) - 1): if np.sign(vals[i]) != np.sign(vals[i + 1]): R2 = brentq(residual, scan[i], scan[i + 1]) break if R2 is None: raise RuntimeError("could not bracket a Lens-2 solution for the target fiber waist") else: R2 = brentq(residual, lo, hi) standoff2, w0 = waist_after_lens2(R1, h1, gap, R2, w0_wg) standoff2 = max(standoff2, 0.0) w, R, _ = propagate_system(R1, h1, R2, gap, standoff2, w0_wg) eta = gaussian_overlap(w, R, FIBER_W0, np.inf, n_common=N_LENS) return dict(R2=R2, standoff2=standoff2, w_at_fiber=w, R_at_fiber=R, eta_gaussian=eta) # ============================================================================ # 4. Exact aspheric surfaces + sequential ray tracer (Snell's law) # ============================================================================ def conic_const(n_in, n_out): """Exact aplanatic (zero spherical aberration) conic constant for a point source in n_in collimated into n_out through a single refracting surface.""" return -(n_out / n_in) ** 2 def sag(r, R, k, A4=0.0): c = 1.0 / R disc = 1 - (1 + k) * c**2 * r**2 disc = np.maximum(disc, 1e-9) return c * r**2 / (1 + np.sqrt(disc)) + A4 * r**4 def surface_r_max_valid(R, k): c = 1.0 / R if 1 + k <= 0: return np.inf return 0.999 / (c * np.sqrt(1 + k)) def refract_ray_vec(d_in, normal_out, n_in, n_out): """Vector Snell's law. d_in, normal_out unit vectors; normal_out points from n_in medium into n_out medium. Returns unit refracted direction or None (TIR).""" d_in = d_in / np.linalg.norm(d_in) normal_out = normal_out / np.linalg.norm(normal_out) cosi = -np.dot(normal_out, d_in) if cosi < 0: normal_out = -normal_out cosi = -np.dot(normal_out, d_in) eta = n_in / n_out k2 = 1 - eta**2 * (1 - cosi**2) if k2 < 0: return None cost = np.sqrt(k2) d_out = eta * d_in + (eta * cosi - cost) * normal_out return d_out / np.linalg.norm(d_out) def surface_height(r, R, k, A4, apex_z, orientation): """orientation=+1 (Lens-1-style): dome apex is the FARTHEST point from its own source, height decreases outward -> height(r) = apex_z - sag(r). orientation=-1 (Lens-2-style, mirror image): dome apex is the point CLOSEST to its incoming beam (faces the gap), height increases outward -> apex_z + sag(r).""" return apex_z - orientation * sag(r, R, k, A4) def surface_normal_outward(r, R, k, A4, orientation, dr=1e-5): """Outward normal (pointing from the lens solid into the medium beyond the surface), consistent with surface_height()'s orientation convention. Verified against the aplanatic-conic reference implementation in dbt_tir_sil_fiber_coupler_v1 / v2_lens_codesign/dual/lens_asphere.py (trace(): normal=(sag', 1) for the orientation=+1 / Lens-1 case).""" sp = (sag(r + dr, R, k, A4) - sag(max(r - dr, 0), R, k, A4)) / (dr if r < dr else 2 * dr) if orientation > 0: n = np.array([sp, 1.0]) else: n = np.array([sp, -1.0]) return n / np.linalg.norm(n) def intersect_surface(origin, direction, R, k, A4, apex_z, orientation, z_search_max): """origin, direction in the meridional (r,z) plane (2D), surface per surface_height(). Bisection along the ray parameter t for the ray-surface crossing.""" def f(t): p = origin + t * direction r = abs(p[0]) return p[1] - surface_height(r, R, k, A4, apex_z, orientation) t_lo, t_hi = 1e-6, z_search_max / max(abs(direction[1]), 1e-6) * 3 f_lo, f_hi = f(t_lo), f(t_hi) if np.sign(f_lo) == np.sign(f_hi): return None for _ in range(60): t_mid = 0.5 * (t_lo + t_hi) if np.sign(f(t_mid)) == np.sign(f_lo): t_lo = t_mid else: t_hi = t_mid t = 0.5 * (t_lo + t_hi) return origin + t * direction def trace_ray_points(theta_launch, field_offset, R1, k1, A4_1, h1, gap, R2, k2, A4_2, standoff2, n_lens=N_LENS): """2D meridional ray trace: point source at (field_offset, 0) inside the resin at the WG facet -> Lens1 (resin/air, apex at z=h1, orientation +1) -> gap -> Lens2 (air/resin, apex at z=h1+gap, orientation -1, mirror-facing L1) -> fiber plane at z = h1+gap+apex2_sag+standoff2. Returns the list of waypoints [origin, p1, p2, p_fiber], or None if TIR/miss.""" origin = np.array([field_offset, 0.0]) direction = np.array([np.sin(theta_launch), np.cos(theta_launch)]) p1 = intersect_surface(origin, direction, R1, k1, A4_1, apex_z=h1, orientation=+1, z_search_max=h1 * 1.5) if p1 is None or p1[1] > h1 * 1.3 or p1[1] < 0: return None r1 = abs(p1[0]) normal1 = surface_normal_outward(r1, R1, k1, A4_1, orientation=+1) if p1[0] < 0: normal1 = np.array([-normal1[0], normal1[1]]) d2 = refract_ray_vec(direction, normal1, n_lens, N_AIR) if d2 is None: return None z2_vertex = h1 + gap rho2_edge = surface_r_max_valid(R2, k2) apex2_sag = sag(min(rho2_edge, 50.0), R2, k2) p2 = intersect_surface(p1, d2, R2, k2, A4_2, apex_z=z2_vertex, orientation=-1, z_search_max=(z2_vertex - p1[1]) + apex2_sag + 5) if p2 is None: return None r2 = abs(p2[0]) normal2 = surface_normal_outward(r2, R2, k2, A4_2, orientation=-1) if p2[0] < 0: normal2 = np.array([-normal2[0], normal2[1]]) d3 = refract_ray_vec(d2, normal2, N_AIR, n_lens) if d3 is None: return None z_fiber = z2_vertex + apex2_sag + standoff2 if abs(d3[1]) < 1e-9: return None t = (z_fiber - p2[1]) / d3[1] p_fiber = p2 + t * d3 return [origin, p1, p2, p_fiber] def trace_ray_through_relay(theta_launch, field_offset, R1, k1, A4_1, h1, gap, R2, k2, A4_2, standoff2, n_lens=N_LENS): """As trace_ray_points(), but returns only the final (r, z) hit point at the fiber plane (or None).""" pts = trace_ray_points(theta_launch, field_offset, R1, k1, A4_1, h1, gap, R2, k2, A4_2, standoff2, n_lens=n_lens) return None if pts is None else pts[-1] def rms_spot_at_fiber(theta_max, n_rays, field_offset, design, n_field=1): thetas = np.linspace(-theta_max, theta_max, n_rays) hits = [] for th in thetas: p = trace_ray_through_relay(th, field_offset, **design) if p is not None: hits.append(p[0]) if len(hits) < 2: return np.nan, 0 hits = np.array(hits) centroid = hits.mean() rms = np.sqrt(np.mean((hits - centroid) ** 2)) return rms, len(hits) # ============================================================================ # 5. STL export -- solid of revolution, numpy-stl only (no trimesh dependency) # ============================================================================ def revolve_to_stl(r_max, R, k, A4, base_thickness, n_r, n_theta, out_path: Path, apex_up=True): """Build a watertight solid-of-revolution STL for a plano-convex asphere: flat base of given thickness, curved top surface z=sag(r) (clipped to r_max).""" r_vals = np.linspace(0, r_max, n_r) theta_vals = np.linspace(0, 2 * np.pi, n_theta, endpoint=False) z_top = sag(r_vals, R, k, A4) if not apex_up: z_top = z_top[::-1] * 0 + (z_top.max() - z_top) # unused path, kept simple # vertices: [base ring*n_r_outer... ] simpler approach: build as stacked rings # top surface vertices (including center apex point at r=0) top_verts = [] for ri, r in enumerate(r_vals): z = z_top[ri] + base_thickness if r < 1e-9: top_verts.append([(0.0, 0.0, z)]) else: ring = [(r * np.cos(t), r * np.sin(t), z) for t in theta_vals] top_verts.append(ring) bot_z = 0.0 bot_ring = [(r_max * np.cos(t), r_max * np.sin(t), bot_z) for t in theta_vals] triangles = [] # side wall of the flat base (cylinder from bot ring to the r_max top ring) top_rim = top_verts[-1] for j in range(n_theta): jn = (j + 1) % n_theta a, b = bot_ring[j], bot_ring[jn] c, d = top_rim[j], top_rim[jn] triangles.append([a, b, c]) triangles.append([b, d, c]) # bottom cap (fan from center) center_bot = (0.0, 0.0, bot_z) for j in range(n_theta): jn = (j + 1) % n_theta triangles.append([center_bot, bot_ring[jn], bot_ring[j]]) # curved top surface, ring-to-ring for ri in range(len(r_vals) - 1): ring_a = top_verts[ri] ring_b = top_verts[ri + 1] if len(ring_a) == 1: apex = ring_a[0] for j in range(n_theta): jn = (j + 1) % n_theta triangles.append([apex, ring_b[j], ring_b[jn]]) else: for j in range(n_theta): jn = (j + 1) % n_theta a, b = ring_a[j], ring_a[jn] c, d = ring_b[j], ring_b[jn] triangles.append([a, c, b]) triangles.append([b, c, d]) tri_arr = np.array(triangles, dtype=np.float64) solid = stlmesh.Mesh(np.zeros(tri_arr.shape[0], dtype=stlmesh.Mesh.dtype)) solid.vectors[:] = tri_arr solid.save(str(out_path)) return solid def build_support_ring_stl(r_inner, r_outer, height, n_theta, out_path: Path): """Simple annular support ring (Nanoscribe-style scaffold) for the steep, tall Lens-1 print, standing off the substrate to key the print in place / relieve write-order proximity effects at the steep flank.""" theta_vals = np.linspace(0, 2 * np.pi, n_theta, endpoint=False) inner = [(r_inner * np.cos(t), r_inner * np.sin(t)) for t in theta_vals] outer = [(r_outer * np.cos(t), r_outer * np.sin(t)) for t in theta_vals] triangles = [] for j in range(n_theta): jn = (j + 1) % n_theta i0, i1, o0, o1 = inner[j], inner[jn], outer[j], outer[jn] # bottom triangles.append([(*i0, 0.0), (*i1, 0.0), (*o0, 0.0)]) triangles.append([(*i1, 0.0), (*o1, 0.0), (*o0, 0.0)]) # top triangles.append([(*i0, height), (*o0, height), (*i1, height)]) triangles.append([(*i1, height), (*o0, height), (*o1, height)]) # outer wall triangles.append([(*o0, 0.0), (*o1, 0.0), (*o0, height)]) triangles.append([(*o1, 0.0), (*o1, height), (*o0, height)]) # inner wall triangles.append([(*i0, 0.0), (*i0, height), (*i1, 0.0)]) triangles.append([(*i1, 0.0), (*i0, height), (*i1, height)]) tri_arr = np.array(triangles, dtype=np.float64) solid = stlmesh.Mesh(np.zeros(tri_arr.shape[0], dtype=stlmesh.Mesh.dtype)) solid.vectors[:] = tri_arr solid.save(str(out_path)) return solid # ============================================================================ # 6. Main # ============================================================================ def main(): ts = dt.datetime.now().strftime("%Y%m%d_%H%M%S") print("=" * 78) print("v3 fiber-coupling lens relay -- waveguide facet -> 2PP asphere pair -> 630HP fiber") print("=" * 78) # --- 1. waveguide mode --- mode = solve_waveguide_mode() print(f"\n[1] Waveguide mode ({WG_WIDTH*1e3:.0f} x {WG_THICKNESS*1e3:.0f} nm Si3N4 strip, air-clad):") print(f" vertical-slab neff (infinite width) = {mode['neff_vertical_slab']:.4f}") print(f" 2D EIM neff (TE0-like) = {mode['neff_2d_eim']:.4f}") print(f" FDTD reference neff (400x300nm, anthracene_slab_wg_coupling_v1) = {FDTD_NEFF_REFERENCE:.4f}") print(f" EIM vs FDTD difference = {mode['eim_vs_fdtd_pct_diff']:+.2f}%") print(f" MFD_y (horizontal) = {mode['mfd_y']*1e3:.0f} nm, MFD_z (vertical) = {mode['mfd_z']*1e3:.0f} nm") print(f" circularized MFD = {mode['mfd_circularized']*1e3:.0f} nm -> w0_wg = {mode['w0_wg']*1e3:.0f} nm") w0_wg = mode["w0_wg"] theta0_wg = LAM0 / (np.pi * N_LENS * w0_wg) zR_wg = np.pi * w0_wg**2 * N_LENS / LAM0 print(f" Gaussian divergence half-angle (in resin) = {np.degrees(theta0_wg):.1f} deg," f" Rayleigh range zR = {zR_wg:.3f} um") # --- 2. fiber target --- theta0_fiber = LAM0 / (np.pi * N_AIR * FIBER_W0) print(f"\n[2] Target: 630HP fiber, MFD = {FIBER_MFD} um (w0 = {FIBER_W0} um)," f" datasheet NA = {FIBER_NA_DATASHEET}") print(f" Gaussian-mode-equivalent divergence half-angle = {np.degrees(theta0_fiber):.2f} deg" f" (vs. datasheet NA half-angle {np.degrees(np.arcsin(FIBER_NA_DATASHEET)):.2f} deg -- " f"expected: datasheet NA is the larger 1%-power definition)") # --- 3. lens design --- R1 = solve_R1(w0_wg, H1) k1 = conic_const(N_LENS, N_AIR) l2 = design_lens2(R1, H1, GAP, w0_wg) R2, standoff2, w_fiber_plane, R_fiber_plane, eta_gauss = ( l2["R2"], l2["standoff2"], l2["w_at_fiber"], l2["R_at_fiber"], l2["eta_gaussian"]) k2 = conic_const(N_AIR, N_LENS) rho1_natural = surface_r_max_valid(R1, k1) rho2_natural = surface_r_max_valid(R2, k2) # apertures actually used: L1 -> its own natural edge (largest useful NA); # L2 -> clipped to whichever is smaller of its natural edge and the L1 output beam # radius at the L2 vertex (avoid vignetting *and* avoid an oversized L2 aperture) w_at_l2, _, _ = propagate_system(R1, H1, 1e9, GAP, 0.0, w0_wg) # beam radius at L2 plane (~no lens) rho1_use = rho1_natural rho2_use = min(rho2_natural, 3.0 * w_at_l2) print(f"\n[3] Lens design (IP-S resin, n={N_LENS} @ {LAM0*1e3:.0f} nm; " f"exact aplanatic conic on each surface):") print(f" Lens 1 (on waveguide facet): R1={R1:.3f} um, k1={k1:.4f}, " f"h_apex={H1:.1f} um, aperture radius (natural)={rho1_natural:.3f} um") print(f" free-space relay gap = {GAP:.1f} um") print(f" Lens 2 (on fiber tip): R2={R2:.3f} um, k2={k2:.4f}, " f"standoff to fiber facet={standoff2:.3f} um, aperture used={rho2_use:.3f} um " f"(natural edge {rho2_natural:.3f} um)") print(f" beam radius at Lens-2 vertex plane (from Lens 1) = {w_at_l2:.3f} um") print(f" predicted spot at fiber plane: w = {w_fiber_plane:.3f} um " f"(target w0_fiber = {FIBER_W0:.3f} um), wavefront R = {R_fiber_plane}") print(f" Gaussian-mode overlap efficiency (paraxial ABCD, waist+curvature match) " f"= {eta_gauss*100:.1f}%") # Fresnel transmission budget (normal-incidence estimate at each interface) T_wg_resin = fresnel_T_normal(mode["neff_2d_eim"], N_LENS) T_l1 = fresnel_T_normal(N_LENS, N_AIR) T_l2 = fresnel_T_normal(N_AIR, N_LENS) T_resin_fiber = fresnel_T_normal(N_LENS, FIBER_N_CORE) T_total = T_wg_resin * T_l1 * T_l2 * T_resin_fiber eta_combined = eta_gauss * T_total print(f"\n Fresnel transmission (normal-incidence estimate, uncoated):") print(f" WG(neff={mode['neff_2d_eim']:.3f})->resin: {T_wg_resin*100:.1f}% " f"resin->air (L1): {T_l1*100:.1f}% air->resin (L2): {T_l2*100:.1f}% " f"resin->fiber-core: {T_resin_fiber*100:.1f}%") print(f" combined Fresnel transmission = {T_total*100:.1f}%") print(f" COMBINED estimate (Gaussian-mode overlap x normal-incidence Fresnel) " f"= {eta_combined*100:.1f}%") print(" *** This is a paraxial/normal-incidence approximation, NOT a rigorous ***") print(" *** coupling-efficiency number -- see README caveats. ***") # --- 4. exact ray trace: on-axis verification, A4 aberration correction, and # misalignment/defocus sweeps --- design = dict(R1=R1, k1=k1, A4_1=0.0, h1=H1, gap=GAP, R2=R2, k2=k2, A4_2=0.0, standoff2=standoff2) # angle from the source (WG facet) to Lens 1's aperture edge = acceptance half-angle z_edge = H1 - sag(rho1_natural, R1, k1) theta_accept = np.arctan2(rho1_natural, z_edge) print(f"\n[4] Ray trace: on-axis point-source acceptance half-angle captured by Lens 1's " f"full natural aperture = {np.degrees(theta_accept):.1f} deg") rms0_before, n_hit0 = rms_spot_at_fiber(theta_accept * 0.98, 41, 0.0, design) print(f" on-axis, full-aperture ray fan ({n_hit0}/41 rays), A4_2=0 (pure paraxial-" f"matched R2 + exact point-source conic k2): RMS spot at fiber plane = " f"{rms0_before*1e3:.2f} nm") print(f" Lens 1 alone is EXACTLY aplanatic for its own point source at any angle " f"(that's what k1 guarantees); this residual comes from Lens 2, whose R2 was " f"matched via PARAXIAL Gaussian propagation, not an exact wide-angle condition " f"-- i.e. real, ray-trace-revealed spherical aberration, not numerical noise.") # optimize a quartic (A4) correction on Lens 2 to reduce this real aberration -- # exactly the "optimize the aspheric surface profile via ray tracing" step. def rms_for_A4(A4_2): d = dict(design); d["A4_2"] = A4_2 rms, n = rms_spot_at_fiber(theta_accept * 0.98, 41, 0.0, d) return rms if n >= 30 else 1e3 A4_scale = 1.0 / rho2_use**3 # natural scale for the quartic term at this aperture res_A4 = minimize(lambda x: rms_for_A4(x[0] * A4_scale), x0=[0.0], method="Nelder-Mead", options=dict(xatol=1e-6, fatol=1e-9, maxiter=200)) A4_2_opt = res_A4.x[0] * A4_scale design["A4_2"] = A4_2_opt rms0, n_hit0 = rms_spot_at_fiber(theta_accept * 0.98, 41, 0.0, design) print(f" after optimizing Lens 2's quartic term (A4_2 = {A4_2_opt:.3e} um^-3): " f"RMS spot = {rms0*1e3:.2f} nm ({100*(1-rms0/max(rms0_before,1e-9)):.0f}% reduction)") # off-axis field point at +-w0_wg (finite source size / coma-like blur) rms_field, n_field = rms_spot_at_fiber(theta_accept * 0.9, 31, w0_wg, design) print(f" off-axis field point (+{w0_wg*1e3:.0f} nm, i.e. mode-edge launch): " f"RMS spot radius = {rms_field*1e3:.2f} nm (this is the real, finite-source-size " f"blur a point-source-designed conic cannot remove)") # misalignment sweeps # NOTE: modeled as a source-offset proxy (see comment below), which only stays # geometrically valid while the offset is small compared to Lens 1's own aperture # (rho1 ~ 2.7 um); keep the ray-traced sweep within that regime. lateral_offsets = np.linspace(0, 1.2, 13) # um, L1-L2 axis lateral misalignment proxy gap_errors = np.linspace(-5.0, 5.0, 13) # um, gap length error (defocus) rms_lateral, rms_gap = [], [] for dx in lateral_offsets: # model as a lateral shift of the source relative to lens 1's axis is degenerate # with shifting lens2's axis; use field_offset as the proxy for relative decenter rms, n = rms_spot_at_fiber(theta_accept * 0.8, 25, dx, design) rms_lateral.append(rms if n > 5 else np.nan) for dg in gap_errors: d2 = dict(design); d2["gap"] = GAP + dg rms, n = rms_spot_at_fiber(theta_accept * 0.8, 25, 0.0, d2) rms_gap.append(rms if n > 5 else np.nan) rms_lateral = np.array(rms_lateral) rms_gap = np.array(rms_gap) # approximate coupling-efficiency-vs-perturbation using the Gaussian overlap with a # lateral-shift penalty (standard shifted-Gaussian overlap) and gap-induced defocus def eta_vs_lateral(dx): w = w_fiber_plane return gaussian_overlap(w, R_fiber_plane, FIBER_W0, np.inf, n_common=N_LENS) * \ np.exp(-2 * dx**2 / (w**2 + FIBER_W0**2)) def eta_vs_gap(dg): w, R, _ = propagate_system(R1, H1, R2, GAP + dg, standoff2, w0_wg) return gaussian_overlap(w, R, FIBER_W0, np.inf, n_common=N_LENS) eta_lateral_curve = np.array([eta_vs_lateral(dx) for dx in lateral_offsets]) eta_gap_curve = np.array([eta_vs_gap(dg) for dg in gap_errors]) # --- 5. figures --- fig_geom_path = FIG_DIR / f"geometry_{ts}.png" fig_sens_path = FIG_DIR / f"sensitivity_{ts}.png" fig_mode_path = FIG_DIR / f"waveguide_mode_{ts}.png" plot_geometry(design, mode, w0_wg, theta_accept, fig_geom_path, rho2_use=rho2_use) plot_sensitivity(lateral_offsets, rms_lateral, eta_lateral_curve, gap_errors, rms_gap, eta_gap_curve, fig_sens_path) plot_waveguide_mode(mode, fig_mode_path) print(f"\n[5] Figures written: {fig_geom_path.name}, {fig_sens_path.name}, {fig_mode_path.name}") # --- 6. STL export --- stl1_path = STL_DIR / "lens1_waveguide_facet.stl" stl2_path = STL_DIR / "lens2_fiber_tip.stl" support_path = STL_DIR / "lens1_support_ring.stl" revolve_to_stl(rho1_use, R1, k1, 0.0, base_thickness=1.0, n_r=100, n_theta=96, out_path=stl1_path) revolve_to_stl(rho2_use, R2, k2, design["A4_2"], base_thickness=1.0, n_r=100, n_theta=96, out_path=stl2_path) # steep-lens support: L1's edge slope is checked below and a scaffold ring is always # written (cheap) sized to L1's footprint; whether it is *needed* is reported. max_slope_deg = max_surface_slope_deg(rho1_use, R1, k1) build_support_ring_stl(r_inner=rho1_use * 1.15, r_outer=rho1_use * 1.6, height=H1 * 0.4, n_theta=64, out_path=support_path) print(f"\n[6] STL exported: {stl1_path.name} (Lens 1, apex height {H1:.1f} um, " f"aperture r={rho1_use:.2f} um, max local flank angle from vertical " f"{max_slope_deg:.0f} deg), {stl2_path.name} (Lens 2, apex height " f"{sag(rho2_use, R2, k2):.2f} um, aperture r={rho2_use:.2f} um), " f"{support_path.name} (support ring stub)") if max_slope_deg > 50: print(f" Lens 1's flank exceeds a 50 deg overhang guideline " f"({max_slope_deg:.0f} deg from vertical) -- support ring RECOMMENDED for print.") else: print(f" Lens 1's flank ({max_slope_deg:.0f} deg from vertical) is within typical " f"2PP self-supporting overhang limits -- support ring included as a conservative " f"anchor/proximity-relief feature, not strictly required.") # --- 7. data + report --- results = dict( timestamp=ts, wavelength_um=LAM0, waveguide=dict(width_um=WG_WIDTH, thickness_um=WG_THICKNESS, n_core=N_SIN, n_top=N_AIR, n_bot=N_SIO2, neff_vertical_slab=mode["neff_vertical_slab"], neff_2d_eim=mode["neff_2d_eim"], fdtd_reference_neff=FDTD_NEFF_REFERENCE, eim_vs_fdtd_pct_diff=mode["eim_vs_fdtd_pct_diff"], mfd_y_um=mode["mfd_y"], mfd_z_um=mode["mfd_z"], mfd_circularized_um=mode["mfd_circularized"], w0_wg_um=w0_wg), fiber=dict(type="630HP (reused from dbt_tir_sil_fiber_coupler_v1 / v2_lens_codesign)", mfd_um=FIBER_MFD, w0_um=FIBER_W0, na_datasheet=FIBER_NA_DATASHEET, n_core=FIBER_N_CORE, n_clad=FIBER_N_CLAD), lens_material=dict(name="Nanoscribe IP-S (assumed)", n=N_LENS), lens1=dict(R_um=R1, k=k1, A4=0.0, h_apex_um=H1, aperture_radius_um=rho1_use, natural_aperture_radius_um=rho1_natural, max_flank_angle_deg=max_slope_deg), gap_um=GAP, lens2=dict(R_um=R2, k=k2, A4=design["A4_2"], standoff_to_fiber_um=standoff2, aperture_radius_um=rho2_use, natural_aperture_radius_um=rho2_natural), performance=dict( w_at_fiber_plane_um=w_fiber_plane, gaussian_mode_overlap_efficiency=eta_gauss, fresnel_transmission=dict(wg_to_resin=T_wg_resin, lens1_resin_to_air=T_l1, lens2_air_to_resin=T_l2, resin_to_fiber_core=T_resin_fiber, total=T_total), combined_estimate_efficiency=eta_combined, onaxis_rms_spot_um=rms0, offaxis_field_rms_spot_um=rms_field, lateral_offsets_um=lateral_offsets.tolist(), rms_spot_vs_lateral_um=[None if np.isnan(v) else v for v in rms_lateral], eta_vs_lateral=eta_lateral_curve.tolist(), gap_errors_um=gap_errors.tolist(), rms_spot_vs_gap_error_um=[None if np.isnan(v) else v for v in rms_gap], eta_vs_gap_error=eta_gap_curve.tolist(), ), caveats=[ "Waveguide mode solved by the effective-index method (1D slab x2), not a full " "2D vectorial mode solver; validated to within " f"{abs(mode['eim_vs_fdtd_pct_diff']):.1f}% of an existing FDTD mode solve in this " "repo (anthracene_slab_wg_coupling_v1) for the same waveguide cross-section.", "The waveguide mode's circularized waist (w0~0.20 um) is only ~0.4x the " "wavelength-in-resin (lambda/n=0.52um); its Rayleigh range (~0.25um) is far " "smaller than the print working distances used here. Paraxial Gaussian-beam " "propagation and geometric point-source ray tracing are both rough " "approximations in this near-field regime -- true performance requires a " "full-wave (FDTD/BPM/eigenmode-expansion) launch, not ray/Gaussian optics.", "Coupling efficiency reported here is a paraxial Gaussian-mode-overlap estimate " "(waist + wavefront-curvature matching) times a normal-incidence Fresnel " "transmission estimate. It ignores angle-dependent Fresnel loss (relevant " "given the >45 deg ray angles at Lens 1), the waveguide mode's non-Gaussian, " "elliptical shape (the design uses a circularized/rotationally-symmetric " "approximation of an intrinsically asymmetric mode), and any coating.", "Lens 1 is EXACTLY aplanatic for its own idealized point source at any ray angle " "(by construction: k1 is the exact collimating conic constant). Lens 2's vertex " "radius was matched via PARAXIAL Gaussian propagation, not an exact wide-angle " "condition, so the full-aperture ray trace reveals genuine (not numerical-noise) " "residual spherical aberration -- " f"{rms0*1e3:.0f} nm RMS on-axis after a quartic (A4) " "correction on Lens 2 (vs. a target fiber waist of " f"{FIBER_W0*1e3:.0f} nm, so a real but secondary effect). A fuller wide-angle " "optimization (more polynomial terms, or solving Lens 2 as an exact finite-" "conjugate conic) would reduce this further. Ray tracing by itself still does " "NOT produce a meaningful diffraction-limited efficiency number -- that is what " "the Gaussian-overlap estimate above is for, with its own caveats.", "The lateral-misalignment sweep models a Lens1-Lens2 decenter as an equivalent " "source-point offset (a standard small-perturbation proxy), which is only valid " "while the offset stays well inside Lens 1's own aperture (~2.7 um); it is not a " "full independent-decenter tolerance analysis of the two separately-printed/" "positioned optics.", ], ) data_path = DATA_DIR / f"results_{ts}.json" with open(data_path, "w") as f: json.dump(results, f, indent=2) print(f"\n[7] Results written: {data_path.relative_to(HERE)}") write_report(ts, results, fig_geom_path, fig_sens_path, fig_mode_path) print(f" Report written: reports/report_{ts}.md") print("\nDone.") return results def max_surface_slope_deg(r_max, R, k, n=400): r = np.linspace(1e-4, r_max * 0.999, n) dr = r[1] - r[0] z = sag(r, R, k) slope = np.gradient(z, dr) # dz/dr, local tangent slope # angle of the local surface tangent measured from VERTICAL (0 deg = vertical wall, # 90 deg = flat horizontal) -- overhang concern grows as this angle -> 90 deg for an # inward-curling profile; for a monotonic outward dome the flank steepens toward the # rim, so report the angle from vertical at the rim as the printability figure. angle_from_vertical = np.degrees(np.arctan2(1.0, np.abs(slope))) return float(np.min(angle_from_vertical)) # smallest angle-from-vertical = steepest point # ============================================================================ # 7. Plotting # ============================================================================ def plot_geometry(design, mode, w0_wg, theta_accept, out_path: Path, rho2_use=None): """Two-panel schematic: the full path (h1 + gap + L2-height + standoff2 ~ 90 um) is dominated by the long low-NA standoff between Lens 2 and the fiber, so a single equal-aspect plot would make the lenses illegible. Panel A zooms on Lens 1 at the waveguide facet (the wide-angle launch); panel B zooms on Lens 2 and the fiber (the gentle refocusing stage); the long uneventful collimated-beam run between them is called out with a break marker rather than drawn to scale.""" R1, k1, A4_1, h1 = design["R1"], design["k1"], design["A4_1"], design["h1"] gap, R2, k2, A4_2, standoff2 = design["gap"], design["R2"], design["k2"], design["A4_2"], design["standoff2"] rho1 = surface_r_max_valid(R1, k1) rho2 = rho2_use if rho2_use is not None else min(surface_r_max_valid(R2, k2), 3.0 * rho1) z2_vertex = h1 + gap apex2 = sag(rho2, R2, k2, A4_2) z_fiber = z2_vertex + apex2 + standoff2 fig, (axA, axB) = plt.subplots(1, 2, figsize=(12.5, 5.5)) # --- Panel A: Lens 1 on the waveguide facet --- axA.add_patch(plt.Rectangle((-1.5, -1.5), 1.5, 3.0, color=MATERIAL_COLORS["sin"], alpha=0.85, label="Si3N4 waveguide")) axA.text(-2.3, 0.0, "waveguide\nfacet", fontsize=8, ha="center", va="center") r1 = np.linspace(0, rho1, 200) z1 = h1 - sag(r1, R1, k1, A4_1) axA.fill_between(np.r_[-r1[::-1], r1], np.r_[z1[::-1], z1], 0, color=MATERIAL_COLORS["resist"], alpha=0.6, label="Lens 1 (IP-S, on facet)") thetas = np.linspace(-theta_accept * 0.95, theta_accept * 0.95, 17) for th in thetas: pts = trace_ray_points(th, 0.0, R1, k1, A4_1, h1, gap, R2, k2, A4_2, standoff2) if pts is None: continue pts = np.array(pts[:2]) # just source -> Lens-1 hit, for this zoomed panel axA.plot(pts[:, 0], pts[:, 1], color=MATERIAL_COLORS["ray"], lw=0.6, alpha=0.7) # a short collimated stub above Lens 1 showing the ray continues parallel to axis p1 = pts[-1] axA.plot([p1[0], p1[0]], [p1[1], h1 * 1.4], color=MATERIAL_COLORS["ray"], lw=0.6, alpha=0.4) axA.set_title(f"Lens 1: WG mode (w0={w0_wg*1e3:.0f} nm) -> collimated\n" f"(h_apex={h1:.1f} um, aperture r={rho1:.2f} um, capture " f"{np.degrees(theta_accept):.0f} deg half-angle)") axA.set_xlabel("r (um)"); axA.set_ylabel("z (um)") axA.set_xlim(-rho1 * 2.2, rho1 * 2.2) axA.set_ylim(-1.8, h1 * 1.4) axA.set_aspect("equal") axA.legend(loc="lower right", fontsize=7.5) # --- Panel B: Lens 2 + fiber --- r2 = np.linspace(0, rho2, 200) z2_curve_local = sag(r2, R2, k2, A4_2) # local coords: 0 at L2 vertex z_fiber_local = apex2 + standoff2 axB.fill_between(np.r_[-r2[::-1], r2], np.r_[z2_curve_local[::-1], z2_curve_local], z_fiber_local, color=MATERIAL_COLORS["resist"], alpha=0.6, label="Lens 2 (IP-S, on fiber tip)") clad_half = max(rho2 * 1.1, FIBER_W0 * 1.6) axB.add_patch(plt.Rectangle((-clad_half, z_fiber_local), 2 * clad_half, 1.5, color=MATERIAL_COLORS["clad"], alpha=0.8, label="fiber cladding (schematic width)")) axB.add_patch(plt.Rectangle((-FIBER_W0, z_fiber_local), 2 * FIBER_W0, 1.5, color=MATERIAL_COLORS["core"], alpha=0.85, label="630HP fiber core")) for th in thetas: pts = trace_ray_points(th, 0.0, R1, k1, A4_1, h1, gap, R2, k2, A4_2, standoff2) if pts is None: continue pts = np.array(pts[1:]) # Lens-1 hit -> Lens-2 hit -> fiber plane pts_local = pts.copy() pts_local[:, 1] -= z2_vertex axB.plot(pts_local[:, 0], pts_local[:, 1], color=MATERIAL_COLORS["ray"], lw=0.6, alpha=0.7) axB.set_title(f"Lens 2: collimated (r~{rho2:.1f} um) -> fiber waist\n" f"(standoff to fiber = {standoff2:.1f} um; total L1-vertex-to-fiber " f"path = {z_fiber:.0f} um)") axB.set_xlabel("r (um)"); axB.set_ylabel("z from Lens-2 vertex (um)") axB.set_xlim(-clad_half * 1.3, clad_half * 1.3) axB.set_ylim(-2, z_fiber_local * 1.08) # not equal-aspect: this panel's z-extent (dominated by the low-NA standoff run) is # much larger than its r-extent -- forcing equal aspect would make Lens 2 illegible. axB.legend(loc="upper right", fontsize=7.5) fig.suptitle("v3 fiber-coupling lens relay -- waveguide facet -> Lens 1 -> " f"[{gap:.0f} um gap, collimated] -> Lens 2 -> [{standoff2:.0f} um standoff] -> fiber " "(panels zoomed independently; long collimated run not to scale between them)") fig.tight_layout() fig.savefig(out_path) plt.close(fig) def plot_sensitivity(lateral, rms_lateral, eta_lateral, gap_err, rms_gap, eta_gap, out_path): fig, axes = plt.subplots(1, 2, figsize=(11, 4.3)) ax = axes[0] ax.plot(lateral, rms_lateral * 1e3, "o-", color=MATERIAL_COLORS["ray"], label="geometric RMS spot") ax.set_xlabel("Lens1-Lens2 lateral misalignment (um)") ax.set_ylabel("RMS spot radius (nm)", color=MATERIAL_COLORS["ray"]) ax2 = ax.twinx() ax2.plot(lateral, eta_lateral * 100, "s--", color=MATERIAL_COLORS["beam"], label="Gaussian-overlap efficiency estimate") ax2.set_ylabel("coupling efficiency estimate (%)", color=MATERIAL_COLORS["beam"]) ax.set_title("Sensitivity to lateral misalignment") ax.grid(True, alpha=0.28) ax = axes[1] ax.plot(gap_err, rms_gap * 1e3, "o-", color=MATERIAL_COLORS["ray"]) ax.set_xlabel("gap length error (um)") ax.set_ylabel("RMS spot radius (nm)", color=MATERIAL_COLORS["ray"]) ax2 = ax.twinx() ax2.plot(gap_err, eta_gap * 100, "s--", color=MATERIAL_COLORS["beam"]) ax2.set_ylabel("coupling efficiency estimate (%)", color=MATERIAL_COLORS["beam"]) ax.set_title("Sensitivity to gap-length (defocus) error") ax.grid(True, alpha=0.28) fig.suptitle("Spot size (ray trace, left axis) and coupling-efficiency estimate " "(paraxial Gaussian overlap, right axis) vs. assembly perturbation") fig.tight_layout() fig.savefig(out_path) plt.close(fig) def plot_waveguide_mode(mode, out_path): z, Ez = mode["z_profile"] y, Ey = mode["y_profile"] fig, axes = plt.subplots(1, 2, figsize=(10, 4)) axes[0].plot(z * 1e3, (Ez / np.abs(Ez).max()) ** 2, color=MATERIAL_COLORS["beam"]) axes[0].axvspan(0, WG_THICKNESS * 1e3, color=MATERIAL_COLORS["sin"], alpha=0.25, label="Si3N4 core") axes[0].set_xlabel("z (nm)"); axes[0].set_ylabel("normalized intensity") axes[0].set_title(f"Vertical (z) mode, MFD={mode['mfd_z']*1e3:.0f} nm") axes[0].legend(fontsize=8) axes[1].plot(y * 1e3, (Ey / np.abs(Ey).max()) ** 2, color=MATERIAL_COLORS["beam"]) axes[1].axvspan(-WG_WIDTH / 2 * 1e3, WG_WIDTH / 2 * 1e3, color=MATERIAL_COLORS["sin"], alpha=0.25, label="Si3N4 core (EIM)") axes[1].set_xlabel("y (nm)"); axes[1].set_ylabel("normalized intensity") axes[1].set_title(f"Horizontal (y, EIM) mode, MFD={mode['mfd_y']*1e3:.0f} nm") axes[1].legend(fontsize=8) fig.suptitle(f"{WG_WIDTH*1e3:.0f}x{WG_THICKNESS*1e3:.0f} nm Si3N4 strip waveguide TE0 mode @ " f"{LAM0*1e3:.0f} nm (effective-index method)") fig.tight_layout() fig.savefig(out_path) plt.close(fig) def write_report(ts, results, fig_geom, fig_sens, fig_mode): p = REPORT_DIR / f"report_{ts}.md" perf = results["performance"] lines = [ f"# v3 fiber-coupling lens relay -- run {ts}", "", f"Wavelength: {results['wavelength_um']*1e3:.0f} nm. Waveguide: " f"{results['waveguide']['width_um']*1e3:.0f}x{results['waveguide']['thickness_um']*1e3:.0f} nm " f"Si3N4 strip, air-clad. Fiber: {results['fiber']['type']}.", "", "## Waveguide mode", f"- EIM neff (2D) = {results['waveguide']['neff_2d_eim']:.4f}, " f"FDTD reference (anthracene_slab_wg_coupling_v1) = {results['waveguide']['fdtd_reference_neff']:.4f} " f"({results['waveguide']['eim_vs_fdtd_pct_diff']:+.2f}%)", f"- MFD_y = {results['waveguide']['mfd_y_um']*1e3:.0f} nm, " f"MFD_z = {results['waveguide']['mfd_z_um']*1e3:.0f} nm, " f"circularized w0 = {results['waveguide']['w0_wg_um']*1e3:.0f} nm", "", "## Lens prescription", f"- Lens 1 (waveguide facet): R = {results['lens1']['R_um']:.3f} um, " f"k = {results['lens1']['k']:.4f}, apex height = {results['lens1']['h_apex_um']:.1f} um, " f"aperture radius = {results['lens1']['aperture_radius_um']:.2f} um", f"- gap = {results['gap_um']:.1f} um", f"- Lens 2 (fiber tip): R = {results['lens2']['R_um']:.3f} um, " f"k = {results['lens2']['k']:.4f}, standoff to fiber facet = " f"{results['lens2']['standoff_to_fiber_um']:.2f} um, " f"aperture radius = {results['lens2']['aperture_radius_um']:.2f} um", "", "## Predicted performance", f"- Gaussian-mode overlap efficiency (paraxial, waist+curvature match): " f"{perf['gaussian_mode_overlap_efficiency']*100:.1f}%", f"- Combined with normal-incidence Fresnel transmission " f"({perf['fresnel_transmission']['total']*100:.1f}%): " f"{perf['combined_estimate_efficiency']*100:.1f}%", f"- On-axis ray-trace RMS spot (point-source limit): {perf['onaxis_rms_spot_um']*1e3:.2f} nm", f"- Off-axis (mode-edge) ray-trace RMS spot: {perf['offaxis_field_rms_spot_um']*1e3:.2f} nm", "", "## Figures", f"![geometry]({fig_geom.relative_to(REPORT_DIR.parent)})", f"![sensitivity]({fig_sens.relative_to(REPORT_DIR.parent)})", f"![mode]({fig_mode.relative_to(REPORT_DIR.parent)})", "", "## Caveats", ] + [f"- {c}" for c in results["caveats"]] p.write_text("\n".join(lines)) if __name__ == "__main__": main()