#!/usr/bin/env python3 """ Atomizer Zernike HTML Generator - ANNULAR APERTURE VERSION ============================================================ This version properly handles mirrors with a central hole by: 1. Masking out the central obscuration from Zernike fitting 2. Using inner/outer radius ratio (obscuration ratio) in calculations 3. Properly visualizing the annular aperture without filling the hole For M1 Mirror: Inner radius = 271.5mm / 2 = 135.75mm (diameter = 271.5mm) Usage: conda activate atomizer python zernike_html_generator_annular.py "path/to/solution.op2" # Or specify custom inner radius (diameter/2) python zernike_html_generator_annular.py "path/to/solution.op2" --inner-radius 135.75 Author: Atomizer Created: 2025-01-06 """ import sys import os from pathlib import Path from math import factorial from datetime import datetime import numpy as np from numpy.linalg import LinAlgError import argparse # Add Atomizer root to path ATOMIZER_ROOT = Path(__file__).parent.parent if str(ATOMIZER_ROOT) not in sys.path: sys.path.insert(0, str(ATOMIZER_ROOT)) try: import plotly.graph_objects as go from plotly.subplots import make_subplots from matplotlib.tri import Triangulation from pyNastran.op2.op2 import OP2 from pyNastran.bdf.bdf import BDF except ImportError as e: print(f"ERROR: Missing dependency: {e}") print("Run: conda activate atomizer") sys.exit(1) # Import the rigorous OPD extractor try: from optimization_engine.extractors.extract_zernike_figure import ZernikeOPDExtractor USE_OPD_METHOD = True print("[INFO] Using rigorous OPD method (accounts for lateral displacement)") except ImportError: USE_OPD_METHOD = False print("[WARN] OPD extractor not available, falling back to simple Z-only method") # ============================================================================ # Configuration # ============================================================================ N_MODES = 50 AMP = 0.5 # visual scale for residual plot (0.5x = reduced deformation) PANCAKE = 3.0 # Z-axis range multiplier for camera view PLOT_DOWNSAMPLE = 10000 FILTER_LOW_ORDERS = 4 # piston, tip, tilt, defocus # Default inner radius for M1 Mirror central hole (271.5mm diameter -> 135.75mm radius) DEFAULT_INNER_RADIUS_MM = 135.75 # Surface plot style COLORSCALE = 'Turbo' SURFACE_LIGHTING = True SHOW_ZERNIKE_BAR = True REQUIRED_SUBCASES = [90, 20, 40, 60] # Displacement unit in OP2 -> nm scale for WFE = 2*Disp_Z DISP_SRC_UNIT = "mm" NM_PER_UNIT = 1e6 if DISP_SRC_UNIT == "mm" else 1e9 # ============================================================================ # Zernike Math (with annular mask support) # ============================================================================ def noll_indices(j: int): if j < 1: raise ValueError("Noll index j must be >= 1") count = 0 n = 0 while True: if n == 0: ms = [0] elif n % 2 == 0: ms = [0] + [m for k in range(1, n//2 + 1) for m in (-2*k, 2*k)] else: ms = [m for k in range(0, (n+1)//2) for m in (-(2*k+1), (2*k+1))] for m in ms: count += 1 if count == j: return n, m n += 1 def zernike_noll(j, r, th): n, m = noll_indices(j) R = np.zeros_like(r) for s in range((n-abs(m))//2 + 1): c = ((-1)**s * factorial(n-s) / (factorial(s) * factorial((n+abs(m))//2 - s) * factorial((n-abs(m))//2 - s))) R += c * r**(n-2*s) if m == 0: return R return R * (np.cos(m*th) if m > 0 else np.sin(-m*th)) def compute_zernike_coeffs_annular(X, Y, vals, n_modes, inner_radius_mm, chunk_size=100000): """ Fit Zernike coefficients to surface data with ANNULAR APERTURE masking. Points inside the central obscuration (r < inner_radius) are EXCLUDED from fitting. Args: X, Y: Node coordinates (mm) vals: Surface values at each node n_modes: Number of Zernike modes inner_radius_mm: Inner radius of annular aperture (mm) chunk_size: For memory efficiency Returns: (coefficients, R_outer, R_inner, obscuration_ratio) """ Xc, Yc = X - np.mean(X), Y - np.mean(Y) R_outer = float(np.max(np.hypot(Xc, Yc))) # Compute radial distance from center (in mm, before normalization) r_mm = np.hypot(Xc, Yc) # Normalize to unit disk r = (r_mm / R_outer).astype(np.float32) th = np.arctan2(Yc, Xc).astype(np.float32) # Compute normalized inner radius r_inner_normalized = inner_radius_mm / R_outer # ANNULAR MASK: r must be between inner and outer radius # Exclude points inside central hole AND outside outer radius mask = (r >= r_inner_normalized) & (r <= 1.0) & ~np.isnan(vals) n_excluded = np.sum(r < r_inner_normalized) n_total = len(r) print(f" [ANNULAR] Inner radius: {inner_radius_mm:.2f} mm (normalized: {r_inner_normalized:.4f})") print(f" [ANNULAR] Excluded {n_excluded} nodes inside central hole ({100*n_excluded/n_total:.1f}%)") print(f" [ANNULAR] Using {np.sum(mask)} nodes for fitting") if not np.any(mask): raise RuntimeError("No valid points in annular region for Zernike fitting.") idx = np.nonzero(mask)[0] m = int(n_modes) G = np.zeros((m, m), dtype=np.float64) h = np.zeros((m,), dtype=np.float64) v = vals.astype(np.float64) for start in range(0, len(idx), chunk_size): sl = idx[start:start+chunk_size] r_b, th_b, v_b = r[sl], th[sl], v[sl] Zb = np.column_stack([zernike_noll(j, r_b, th_b).astype(np.float32) for j in range(1, m+1)]) G += (Zb.T @ Zb).astype(np.float64) h += (Zb.T @ v_b).astype(np.float64) try: coeffs = np.linalg.solve(G, h) except LinAlgError: coeffs = np.linalg.lstsq(G, h, rcond=None)[0] return coeffs, R_outer, inner_radius_mm, r_inner_normalized def zernike_common_name(n: int, m: int) -> str: names = { (0, 0): "Piston", (1, -1): "Tilt X", (1, 1): "Tilt Y", (2, 0): "Defocus", (2, -2): "Astig 45 deg", (2, 2): "Astig 0 deg", (3, -1): "Coma X", (3, 1): "Coma Y", (3, -3): "Trefoil X", (3, 3): "Trefoil Y", (4, 0): "Primary Spherical", (4, -2): "Sec Astig X", (4, 2): "Sec Astig Y", (4, -4): "Quadrafoil X", (4, 4): "Quadrafoil Y", (5, -1): "Sec Coma X", (5, 1): "Sec Coma Y", (5, -3): "Sec Trefoil X", (5, 3): "Sec Trefoil Y", (5, -5): "Pentafoil X", (5, 5): "Pentafoil Y", (6, 0): "Sec Spherical", } return names.get((n, m), f"Z(n={n}, m={m})") def zernike_label_for_j(j: int) -> str: n, m = noll_indices(j) return f"J{j:02d} - {zernike_common_name(n, m)} (n={n}, m={m})" # ============================================================================ # File I/O # ============================================================================ def find_geometry_file(op2_path: Path) -> Path: """Find matching BDF/DAT file for an OP2.""" folder = op2_path.parent base = op2_path.stem for ext in ['.dat', '.bdf']: cand = folder / (base + ext) if cand.exists(): return cand for f in folder.iterdir(): if f.suffix.lower() in ['.dat', '.bdf']: return f raise FileNotFoundError(f"No .dat or .bdf geometry file found for {op2_path}") def read_geometry(dat_path: Path) -> dict: bdf = BDF() bdf.read_bdf(str(dat_path)) return {int(nid): node.get_position() for nid, node in bdf.nodes.items()} def read_displacements(op2_path: Path) -> dict: """Read displacement data organized by subcase.""" op2 = OP2() op2.read_op2(str(op2_path)) if not op2.displacements: raise RuntimeError("No displacement data found in OP2 file") result = {} for key, darr in op2.displacements.items(): data = darr.data dmat = data[0] if data.ndim == 3 else (data if data.ndim == 2 else None) if dmat is None: continue ngt = darr.node_gridtype.astype(int) node_ids = ngt if ngt.ndim == 1 else ngt[:, 0] isubcase = getattr(darr, 'isubcase', None) label = str(isubcase) if isubcase else str(key) result[label] = { 'node_ids': node_ids.astype(int), 'disp': dmat.copy() } return result # ============================================================================ # Data Processing (with annular support) # ============================================================================ def build_wfe_arrays(label: str, node_ids, dmat, node_geo): """Build X, Y, WFE arrays for a subcase.""" X, Y, WFE = [], [], [] for nid, vec in zip(node_ids, dmat): geo = node_geo.get(int(nid)) if geo is None: continue X.append(geo[0]) Y.append(geo[1]) wfe = vec[2] * 2.0 * NM_PER_UNIT # Z-disp to WFE (nm) WFE.append(wfe) return np.array(X), np.array(Y), np.array(WFE) def compute_relative_wfe(X1, Y1, WFE1, node_ids1, X2, Y2, WFE2, node_ids2): """Compute relative WFE: WFE1 - WFE2 for common nodes.""" ref_map = {int(nid): (x, y, w) for nid, x, y, w in zip(node_ids2, X2, Y2, WFE2)} X_rel, Y_rel, WFE_rel = [], [], [] for nid, x, y, w in zip(node_ids1, X1, Y1, WFE1): nid = int(nid) if nid in ref_map: _, _, w_ref = ref_map[nid] X_rel.append(x) Y_rel.append(y) WFE_rel.append(w - w_ref) return np.array(X_rel), np.array(Y_rel), np.array(WFE_rel) def compute_rms_metrics_annular(X, Y, W_nm, inner_radius_mm): """Compute RMS metrics with ANNULAR APERTURE masking.""" coeffs, R_outer, R_inner, r_inner_norm = compute_zernike_coeffs_annular( X, Y, W_nm, N_MODES, inner_radius_mm ) Xc = X - np.mean(X) Yc = Y - np.mean(Y) r_mm = np.hypot(Xc, Yc) r = r_mm / R_outer th = np.arctan2(Yc, Xc) # Create annular mask for RMS calculation annular_mask = r >= r_inner_norm Z = np.column_stack([zernike_noll(j, r, th) for j in range(1, N_MODES+1)]) # Apply coefficients to get low-order contribution W_low = Z[:, :FILTER_LOW_ORDERS].dot(coeffs[:FILTER_LOW_ORDERS]) W_res_filt = W_nm - W_low # J1-J3 filtered W_j1to3 = Z[:, :3].dot(coeffs[:3]) W_res_filt_j1to3 = W_nm - W_j1to3 # Compute RMS ONLY over annular region (excluding central hole) global_rms = float(np.sqrt(np.mean(W_nm[annular_mask]**2))) filtered_rms = float(np.sqrt(np.mean(W_res_filt[annular_mask]**2))) rms_filter_j1to3 = float(np.sqrt(np.mean(W_res_filt_j1to3[annular_mask]**2))) return { 'coefficients': coeffs, 'R_outer': R_outer, 'R_inner': R_inner, 'obscuration_ratio': r_inner_norm, 'global_rms': global_rms, 'filtered_rms': filtered_rms, 'rms_filter_j1to3': rms_filter_j1to3, 'W_res_filt': W_res_filt, 'annular_mask': annular_mask, 'n_annular_nodes': int(np.sum(annular_mask)), 'n_total_nodes': len(W_nm), } def compute_mfg_metrics(coeffs): """Manufacturing aberration magnitudes from Zernike coefficients.""" defocus = float(abs(coeffs[3])) astigmatism = float(np.sqrt(coeffs[4]**2 + coeffs[5]**2)) coma = float(np.sqrt(coeffs[6]**2 + coeffs[7]**2)) trefoil = float(np.sqrt(coeffs[8]**2 + coeffs[9]**2)) spherical = float(abs(coeffs[10])) if len(coeffs) > 10 else 0.0 higher_order_rms = float(np.sqrt(np.sum(coeffs[3:]**2))) return { 'defocus_nm': defocus, 'astigmatism_rms': astigmatism, 'coma_rms': coma, 'trefoil_rms': trefoil, 'spherical_nm': spherical, 'higher_order_rms': higher_order_rms, } # ============================================================================ # HTML Generation (with annular visualization) # ============================================================================ def generate_html_annular( title: str, X: np.ndarray, Y: np.ndarray, W_nm: np.ndarray, rms_data: dict, inner_radius_mm: float, is_relative: bool = False, ref_title: str = "20 deg", abs_pair: tuple = None, is_manufacturing: bool = False, mfg_metrics: dict = None, correction_metrics: dict = None, ) -> str: """Generate HTML with proper annular aperture visualization.""" coeffs = rms_data['coefficients'] global_rms = rms_data['global_rms'] filtered_rms = rms_data['filtered_rms'] W_res_filt = rms_data['W_res_filt'] annular_mask = rms_data['annular_mask'] labels = [zernike_label_for_j(j) for j in range(1, N_MODES+1)] coeff_abs = np.abs(coeffs) # Apply annular mask BEFORE downsampling X_ann = X[annular_mask] Y_ann = Y[annular_mask] W_ann = W_res_filt[annular_mask] # Downsample for display n = len(X_ann) if n > PLOT_DOWNSAMPLE: rng = np.random.default_rng(42) sel = rng.choice(n, size=PLOT_DOWNSAMPLE, replace=False) Xp, Yp, Wp = X_ann[sel], Y_ann[sel], W_ann[sel] else: Xp, Yp, Wp = X_ann, Y_ann, W_ann res_amp = AMP * Wp max_amp = float(np.max(np.abs(res_amp))) if res_amp.size else 1.0 # Triangulate with constrained triangulation to respect the hole mesh_traces = [] try: tri = Triangulation(Xp, Yp) # Filter out triangles that span across the central hole # A triangle spans the hole if any edge crosses the center if tri.triangles is not None and len(tri.triangles) > 0: # Get triangle centroids tri_x = Xp[tri.triangles].mean(axis=1) tri_y = Yp[tri.triangles].mean(axis=1) # Compute centroid distance from center Xc_mean = np.mean(X) Yc_mean = np.mean(Y) tri_r = np.hypot(tri_x - Xc_mean, tri_y - Yc_mean) # Filter triangles: keep only those with centroid outside inner radius # Also filter triangles that are too large (span across hole) max_edge_length = 2 * inner_radius_mm # Triangles spanning hole would be large valid_triangles = [] for i, t in enumerate(tri.triangles): # Check centroid is outside hole if tri_r[i] < inner_radius_mm * 0.8: # Some margin continue # Check all vertices are outside hole vx = Xp[t] - Xc_mean vy = Yp[t] - Yc_mean vr = np.hypot(vx, vy) if np.any(vr < inner_radius_mm * 0.9): continue # Check triangle isn't spanning the hole (edge length check) p0, p1, p2 = Xp[t] + 1j*Yp[t] edges = [abs(p1-p0), abs(p2-p1), abs(p0-p2)] if max(edges) > max_edge_length: continue valid_triangles.append(t) if valid_triangles: valid_triangles = np.array(valid_triangles) i, j, k = valid_triangles.T mesh_traces.append(go.Mesh3d( x=Xp, y=Yp, z=res_amp, i=i, j=j, k=k, intensity=res_amp, colorscale=COLORSCALE, opacity=1.0, flatshading=False, lighting=dict( ambient=0.4, diffuse=0.8, specular=0.3, roughness=0.5, fresnel=0.2 ), lightposition=dict(x=100, y=200, z=300), showscale=True, colorbar=dict( title=dict(text="Residual (nm)", side="right"), thickness=15, len=0.6, tickformat=".1f" ), hovertemplate="X: %{x:.1f}
Y: %{y:.1f}
Residual: %{z:.2f} nm" )) # Add a circle to show the inner hole boundary theta_circle = np.linspace(0, 2*np.pi, 100) hole_x = Xc_mean + inner_radius_mm * np.cos(theta_circle) hole_y = Yc_mean + inner_radius_mm * np.sin(theta_circle) hole_z = np.zeros_like(hole_x) mesh_traces.append(go.Scatter3d( x=hole_x, y=hole_y, z=hole_z, mode='lines', line=dict(color='white', width=3), name='Central Hole', showlegend=True, hoverinfo='name' )) except Exception as e: print(f"Triangulation warning: {e}") # Fallback scatter if mesh failed if not mesh_traces: mesh_traces.append(go.Scatter3d( x=Xp, y=Yp, z=res_amp, mode='markers', marker=dict(size=2, color=res_amp, colorscale=COLORSCALE, showscale=True), showlegend=False )) title_suffix = f" (relative to {ref_title})" if is_relative else " (absolute)" # Annular info for title obscuration = rms_data.get('obscuration_ratio', 0) * 100 annular_info = f" | Annular: {inner_radius_mm:.1f}mm inner ({obscuration:.1f}% obscuration)" # Build figure layout if is_manufacturing and mfg_metrics and correction_metrics: fig = make_subplots( rows=5, cols=1, specs=[[{"type":"scene"}], [{"type":"table"}], [{"type":"table"}], [{"type":"table"}], [{"type":"xy"}]], row_heights=[0.38, 0.12, 0.12, 0.18, 0.20], vertical_spacing=0.025, subplot_titles=[ f"Surface Residual - {title}{title_suffix}{annular_info}", "RMS Metrics (Annular Aperture)", "Mode Magnitudes at 90 deg", "Pre-Correction (90 deg - 20 deg)", "|Zernike Coefficients| (nm)" ] ) elif SHOW_ZERNIKE_BAR: fig = make_subplots( rows=4, cols=1, specs=[[{"type":"scene"}], [{"type":"table"}], [{"type":"table"}], [{"type":"xy"}]], row_heights=[0.45, 0.12, 0.25, 0.18], vertical_spacing=0.03, subplot_titles=[ f"Surface Residual - {title}{title_suffix}{annular_info}", "RMS Metrics (Annular Aperture)", f"Zernike Coefficients ({N_MODES} modes)", "|Zernike Coefficients| (nm)" ] ) else: fig = make_subplots( rows=3, cols=1, specs=[[{"type":"scene"}], [{"type":"table"}], [{"type":"table"}]], row_heights=[0.55, 0.15, 0.30], vertical_spacing=0.03, subplot_titles=[ f"Surface Residual - {title}{title_suffix}{annular_info}", "RMS Metrics (Annular Aperture)", f"Zernike Coefficients ({N_MODES} modes)" ] ) for tr in mesh_traces: fig.add_trace(tr, row=1, col=1) fig.update_scenes( camera=dict( eye=dict(x=1.2, y=1.2, z=0.8), up=dict(x=0, y=0, z=1) ), xaxis=dict( title="X (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)', showbackground=True, backgroundcolor='rgba(240,240,240,0.9)' ), yaxis=dict( title="Y (mm)", showgrid=True, gridcolor='rgba(128,128,128,0.3)', showbackground=True, backgroundcolor='rgba(240,240,240,0.9)' ), zaxis=dict( title="Residual (nm)", range=[-max_amp * PANCAKE, max_amp * PANCAKE], showgrid=True, gridcolor='rgba(128,128,128,0.3)', showbackground=True, backgroundcolor='rgba(230,230,250,0.9)' ), aspectmode='manual', aspectratio=dict(x=1, y=1, z=0.4) ) # RMS table with annular info n_ann = rms_data.get('n_annular_nodes', 0) n_tot = rms_data.get('n_total_nodes', 0) if is_relative and abs_pair: abs_global, abs_filtered = abs_pair fig.add_trace(go.Table( header=dict(values=["Metric", "Relative (nm)", "Absolute (nm)", "Notes"], align="left", fill_color='#1f2937', font=dict(color='white')), cells=dict(values=[ ["Global RMS", "Filtered RMS (J1-J4 removed)", "Annular Nodes"], [f"{global_rms:.2f}", f"{filtered_rms:.2f}", f"{n_ann}"], [f"{abs_global:.2f}", f"{abs_filtered:.2f}", f"{n_tot} total"], ["Annular mask applied", "Central hole excluded", f"{inner_radius_mm:.1f}mm inner radius"], ], align="left", fill_color='#374151', font=dict(color='white')) ), row=2, col=1) elif is_manufacturing and mfg_metrics and correction_metrics: fig.add_trace(go.Table( header=dict(values=["Metric", "Value (nm)", "Notes"], align="left", fill_color='#1f2937', font=dict(color='white')), cells=dict(values=[ ["Global RMS", "Filtered RMS (J1-J4)", "Annular Nodes"], [f"{global_rms:.2f}", f"{filtered_rms:.2f}", f"{n_ann} / {n_tot}"], ["Annular mask applied", "Central hole excluded", f"Inner R = {inner_radius_mm:.1f}mm"] ], align="left", fill_color='#374151', font=dict(color='white')) ), row=2, col=1) fig.add_trace(go.Table( header=dict(values=["Mode", "Value (nm)"], align="left", fill_color='#1f2937', font=dict(color='white')), cells=dict(values=[ ["MFG_90 Objective (90-20, J1-J3 filtered)", "Astigmatism (J5+J6)", "Coma (J7+J8)", "Trefoil (J9+J10)", "Spherical (J11)"], [f"{rms_data['rms_filter_j1to3']:.2f}", f"{mfg_metrics['astigmatism_rms']:.2f}", f"{mfg_metrics['coma_rms']:.2f}", f"{mfg_metrics['trefoil_rms']:.2f}", f"{mfg_metrics['spherical_nm']:.2f}"] ], align="left", fill_color='#374151', font=dict(color='white')) ), row=3, col=1) fig.add_trace(go.Table( header=dict(values=["Aberration", "Magnitude (nm)"], align="left", fill_color='#1f2937', font=dict(color='white')), cells=dict(values=[ ["Defocus (J4)", "Astigmatism (J5+J6)", "Coma (J7+J8)", "Trefoil (J9+J10)", "Spherical (J11)"], [f"{correction_metrics['defocus_nm']:.2f}", f"{correction_metrics['astigmatism_rms']:.2f}", f"{correction_metrics['coma_rms']:.2f}", f"{correction_metrics['trefoil_rms']:.2f}", f"{correction_metrics['spherical_nm']:.2f}"] ], align="left", fill_color='#374151', font=dict(color='white')) ), row=4, col=1) else: fig.add_trace(go.Table( header=dict(values=["Metric", "Value (nm)", "Notes"], align="left", fill_color='#1f2937', font=dict(color='white')), cells=dict(values=[ ["Global RMS", "Filtered RMS (J1-J4 removed)", "Annular Nodes"], [f"{global_rms:.2f}", f"{filtered_rms:.2f}", f"{n_ann} / {n_tot}"], ["Annular mask applied", "Central hole excluded", f"Inner R = {inner_radius_mm:.1f}mm"] ], align="left", fill_color='#374151', font=dict(color='white')) ), row=2, col=1) # Zernike coefficients table if not (is_manufacturing and mfg_metrics and correction_metrics): fig.add_trace(go.Table( header=dict(values=["Noll j", "Label", "|Coeff| (nm)"], align="left", fill_color='#1f2937', font=dict(color='white')), cells=dict(values=[ list(range(1, N_MODES+1)), labels, [f"{c:.3f}" for c in coeff_abs] ], align="left", fill_color='#374151', font=dict(color='white')) ), row=3, col=1) # Bar chart if SHOW_ZERNIKE_BAR: bar_row = 5 if (is_manufacturing and mfg_metrics and correction_metrics) else 4 fig.add_trace( go.Bar( x=coeff_abs.tolist(), y=labels, orientation='h', marker_color='#6366f1', hovertemplate="%{y}
|Coeff| = %{x:.3f} nm", showlegend=False ), row=bar_row, col=1 ) height = 1500 if (is_manufacturing and mfg_metrics and correction_metrics) else 1300 fig.update_layout( width=1400, height=height, margin=dict(t=60, b=20, l=20, r=20), paper_bgcolor='#111827', plot_bgcolor='#1f2937', font=dict(color='white'), title=dict( text=f"Atomizer Zernike Analysis (ANNULAR) - {title}", x=0.5, font=dict(size=18) ) ) return fig.to_html(include_plotlyjs='cdn', full_html=True) # ============================================================================ # Main # ============================================================================ def find_op2_file(working_dir=None): """Find the most recent OP2 file in the working directory.""" if working_dir is None: working_dir = Path.cwd() else: working_dir = Path(working_dir) op2_files = list(working_dir.glob("*solution*.op2")) + list(working_dir.glob("*.op2")) if not op2_files: op2_files = list(working_dir.glob("**/*solution*.op2")) if not op2_files: return None return max(op2_files, key=lambda p: p.stat().st_mtime) def main_annular(op2_path: Path, inner_radius_mm: float): """Generate all 3 HTML files with ANNULAR aperture handling.""" print("=" * 70) print(" ATOMIZER ZERNIKE HTML GENERATOR (ANNULAR APERTURE)") print("=" * 70) print(f"\nOP2 File: {op2_path.name}") print(f"Directory: {op2_path.parent}") print(f"\n[ANNULAR] Central hole inner radius: {inner_radius_mm:.2f} mm") print(f"[ANNULAR] Central hole diameter: {2*inner_radius_mm:.2f} mm") # Find geometry print("\nFinding geometry file...") geo_path = find_geometry_file(op2_path) print(f"Found: {geo_path.name}") # Read data print("\nReading geometry...") node_geo = read_geometry(geo_path) print(f"Loaded {len(node_geo)} nodes") print("\nReading displacements...") displacements = read_displacements(op2_path) print(f"Found subcases: {list(displacements.keys())}") # Map subcases - find required angles by name (robust to extra subcases) required_angles = ['90', '20', '40', '60'] subcase_map = {} # First, try direct angle matching (e.g., "20", "40", "60", "90") for angle in required_angles: if angle in displacements: subcase_map[angle] = angle # If we didn't find all angles, try numeric IDs (e.g., "1"=90, "2"=20, "3"=40, "4"=60) if len(subcase_map) < 4: if all(str(i) in displacements for i in range(1, 5)): subcase_map = {'90': '1', '20': '2', '40': '3', '60': '4'} print(f"[INFO] Using numeric subcases: 1=90°, 2=20°, 3=40°, 4=60°") # Check if we found all required angles if len(subcase_map) < 4: missing = [a for a in required_angles if a not in subcase_map] print(f"[ERROR] Missing required subcases: {missing}") print(f"[ERROR] Available subcases: {list(displacements.keys())}") print(f"[ERROR] Required subcases must be named: {required_angles} (or numeric 1,2,3,4)") return print(f"[INFO] Using subcase mapping: {subcase_map}") output_dir = op2_path.parent base = op2_path.stem timestamp = datetime.now().strftime("%Y%m%d_%H%M%S") html_files = [] # ======================================================================== # Process subcases with ANNULAR masking # ======================================================================== print("\nProcessing subcases with ANNULAR aperture masking...") # Reference: 20 deg ref_label = subcase_map['20'] ref_data = displacements[ref_label] X_ref, Y_ref, WFE_ref = build_wfe_arrays( '20', ref_data['node_ids'], ref_data['disp'], node_geo ) print(f"\n20 deg (Reference):") rms_ref = compute_rms_metrics_annular(X_ref, Y_ref, WFE_ref, inner_radius_mm) print(f" Global RMS = {rms_ref['global_rms']:.2f} nm, Filtered = {rms_ref['filtered_rms']:.2f} nm") # Manufacturing: 90 deg mfg_label = subcase_map['90'] mfg_data = displacements[mfg_label] X_90, Y_90, WFE_90 = build_wfe_arrays( '90', mfg_data['node_ids'], mfg_data['disp'], node_geo ) print(f"\n90 deg (Manufacturing):") rms_90 = compute_rms_metrics_annular(X_90, Y_90, WFE_90, inner_radius_mm) mfg_metrics = compute_mfg_metrics(rms_90['coefficients']) print(f" Global RMS = {rms_90['global_rms']:.2f} nm, Filtered = {rms_90['filtered_rms']:.2f} nm") # ======================================================================== # 1. Generate 40 deg vs 20 deg (relative) # ======================================================================== print("\n" + "-" * 70) print("Generating 40 deg vs 20 deg...") sc_40_label = subcase_map['40'] sc_40_data = displacements[sc_40_label] X_40, Y_40, WFE_40 = build_wfe_arrays( '40', sc_40_data['node_ids'], sc_40_data['disp'], node_geo ) X_40_rel, Y_40_rel, WFE_40_rel = compute_relative_wfe( X_40, Y_40, WFE_40, sc_40_data['node_ids'], X_ref, Y_ref, WFE_ref, ref_data['node_ids'] ) print(f"40 deg (Relative to 20 deg):") rms_40_abs = compute_rms_metrics_annular(X_40, Y_40, WFE_40, inner_radius_mm) rms_40_rel = compute_rms_metrics_annular(X_40_rel, Y_40_rel, WFE_40_rel, inner_radius_mm) html_40 = generate_html_annular( title="40 deg (Annular)", X=X_40_rel, Y=Y_40_rel, W_nm=WFE_40_rel, rms_data=rms_40_rel, inner_radius_mm=inner_radius_mm, is_relative=True, ref_title="20 deg", abs_pair=(rms_40_abs['global_rms'], rms_40_abs['filtered_rms']) ) path_40 = output_dir / f"{base}_{timestamp}_40_vs_20_ANNULAR.html" path_40.write_text(html_40, encoding='utf-8') html_files.append(path_40) print(f" Created: {path_40.name}") print(f" Relative: Global={rms_40_rel['global_rms']:.2f}, Filtered={rms_40_rel['filtered_rms']:.2f}") # ======================================================================== # 2. Generate 60 deg vs 20 deg (relative) # ======================================================================== print("\n" + "-" * 70) print("Generating 60 deg vs 20 deg...") sc_60_label = subcase_map['60'] sc_60_data = displacements[sc_60_label] X_60, Y_60, WFE_60 = build_wfe_arrays( '60', sc_60_data['node_ids'], sc_60_data['disp'], node_geo ) X_60_rel, Y_60_rel, WFE_60_rel = compute_relative_wfe( X_60, Y_60, WFE_60, sc_60_data['node_ids'], X_ref, Y_ref, WFE_ref, ref_data['node_ids'] ) print(f"60 deg (Relative to 20 deg):") rms_60_abs = compute_rms_metrics_annular(X_60, Y_60, WFE_60, inner_radius_mm) rms_60_rel = compute_rms_metrics_annular(X_60_rel, Y_60_rel, WFE_60_rel, inner_radius_mm) html_60 = generate_html_annular( title="60 deg (Annular)", X=X_60_rel, Y=Y_60_rel, W_nm=WFE_60_rel, rms_data=rms_60_rel, inner_radius_mm=inner_radius_mm, is_relative=True, ref_title="20 deg", abs_pair=(rms_60_abs['global_rms'], rms_60_abs['filtered_rms']) ) path_60 = output_dir / f"{base}_{timestamp}_60_vs_20_ANNULAR.html" path_60.write_text(html_60, encoding='utf-8') html_files.append(path_60) print(f" Created: {path_60.name}") print(f" Relative: Global={rms_60_rel['global_rms']:.2f}, Filtered={rms_60_rel['filtered_rms']:.2f}") # ======================================================================== # 3. Generate 90 deg Manufacturing # ======================================================================== print("\n" + "-" * 70) print("Generating 90 deg Manufacturing...") X_90_rel, Y_90_rel, WFE_90_rel = compute_relative_wfe( X_90, Y_90, WFE_90, mfg_data['node_ids'], X_ref, Y_ref, WFE_ref, ref_data['node_ids'] ) print(f"90 deg (Relative to 20 deg for correction):") rms_90_rel = compute_rms_metrics_annular(X_90_rel, Y_90_rel, WFE_90_rel, inner_radius_mm) correction_metrics = compute_mfg_metrics(rms_90_rel['coefficients']) html_90 = generate_html_annular( title="90 deg Manufacturing (Annular)", X=X_90, Y=Y_90, W_nm=WFE_90, rms_data=rms_90_rel, inner_radius_mm=inner_radius_mm, is_relative=False, is_manufacturing=True, mfg_metrics=mfg_metrics, correction_metrics=correction_metrics ) path_90 = output_dir / f"{base}_{timestamp}_90_mfg_ANNULAR.html" path_90.write_text(html_90, encoding='utf-8') html_files.append(path_90) print(f" Created: {path_90.name}") print(f" Absolute: Global={rms_90['global_rms']:.2f}, Filtered={rms_90['filtered_rms']:.2f}") print(f" Optician Workload (J1-J3): {rms_90['rms_filter_j1to3']:.2f} nm") # ======================================================================== # Summary # ======================================================================== print("\n" + "=" * 70) print("SUMMARY (ANNULAR APERTURE)") print("=" * 70) print(f"\nCentral hole: {2*inner_radius_mm:.1f} mm diameter ({inner_radius_mm:.2f} mm radius)") print(f"Obscuration ratio: {rms_40_rel['obscuration_ratio']*100:.1f}%") print(f"\nGenerated {len(html_files)} HTML files:") for f in html_files: print(f" - {f.name}") print("\n" + "-" * 70) print("OPTIMIZATION OBJECTIVES (ANNULAR)") print("-" * 70) print(f" 40-20 Filtered RMS: {rms_40_rel['filtered_rms']:.2f} nm") print(f" 60-20 Filtered RMS: {rms_60_rel['filtered_rms']:.2f} nm") print(f" MFG 90 (J1-J3): {rms_90_rel['rms_filter_j1to3']:.2f} nm") # Weighted sums ws_v4 = 5*rms_40_rel['filtered_rms'] + 5*rms_60_rel['filtered_rms'] + 2*rms_90_rel['rms_filter_j1to3'] ws_v5 = 5*rms_40_rel['filtered_rms'] + 5*rms_60_rel['filtered_rms'] + 3*rms_90_rel['rms_filter_j1to3'] print(f"\n V4 Weighted Sum (5/5/2): {ws_v4:.2f}") print(f" V5 Weighted Sum (5/5/3): {ws_v5:.2f}") print("\n" + "=" * 70) print("DONE") print("=" * 70) return html_files if __name__ == '__main__': parser = argparse.ArgumentParser( description='Atomizer Zernike HTML Generator - ANNULAR APERTURE VERSION', epilog='For M1 Mirror with 271.5mm central hole diameter, use --inner-radius 135.75' ) parser.add_argument('op2_file', nargs='?', help='Path to OP2 results file') parser.add_argument('--inner-radius', '-r', type=float, default=DEFAULT_INNER_RADIUS_MM, help=f'Inner radius of central hole in mm (default: {DEFAULT_INNER_RADIUS_MM}mm for 271.5mm diameter)') parser.add_argument('--inner-diameter', '-d', type=float, default=None, help='Inner diameter of central hole in mm (alternative to --inner-radius)') args = parser.parse_args() # Handle diameter vs radius inner_radius = args.inner_radius if args.inner_diameter is not None: inner_radius = args.inner_diameter / 2.0 print(f"[INFO] Using inner diameter {args.inner_diameter}mm -> radius {inner_radius}mm") # Find OP2 file if args.op2_file: op2_path = Path(args.op2_file) if not op2_path.exists(): print(f"ERROR: File not found: {op2_path}") sys.exit(1) else: print("No OP2 file specified, searching...") op2_path = find_op2_file() if op2_path is None: print("ERROR: No OP2 file found in current directory.") print("Usage: python zernike_html_generator_annular.py ") sys.exit(1) print(f"Found: {op2_path}") try: main_annular(op2_path, inner_radius) except Exception as e: print(f"\nERROR: {e}") import traceback traceback.print_exc() sys.exit(1)