diff --git a/tools/generate_optical_report.py b/tools/generate_optical_report.py
new file mode 100644
index 00000000..d0bffd4f
--- /dev/null
+++ b/tools/generate_optical_report.py
@@ -0,0 +1,1223 @@
+#!/usr/bin/env python3
+"""
+Atomizer Optical Performance Report Generator
+===============================================
+
+Generates a comprehensive, CDR-ready HTML report for the optical
+performance of an M1 mirror design from FEA results (OP2 file).
+
+The report combines:
+1. Executive Summary with pass/fail vs design targets
+2. Per-Angle Wavefront Error Analysis (3D surface plots)
+3. Zernike Trajectory Analysis (mode-specific metrics across elevation)
+4. Sensitivity Matrix (axial vs lateral load response)
+5. Manufacturing Analysis (90° correction metrics)
+6. Full Zernike coefficient tables
+
+Usage:
+    conda activate atomizer
+    python generate_optical_report.py "path/to/solution.op2"
+
+    # With annular aperture
+    python generate_optical_report.py "path/to/solution.op2" --inner-radius 135.75
+
+    # Custom targets
+    python generate_optical_report.py "path/to/solution.op2" --target-40 4.0 --target-60 10.0 --target-mfg 20.0
+
+    # Include design parameters from study database
+    python generate_optical_report.py "path/to/solution.op2" --study-db "path/to/study.db" --trial 725
+
+Output:
+    Creates a single comprehensive HTML file:
+    {basename}_OPTICAL_REPORT_{timestamp}.html
+
+Author: Atomizer / Atomaste
+Created: 2026-01-29
+"""
+
+import sys
+import os
+import argparse
+import json
+from pathlib import Path
+from math import factorial
+from datetime import datetime
+import numpy as np
+from numpy.linalg import LinAlgError
+
+# 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 Atomizer extractors
+from optimization_engine.extractors.extract_zernike import (
+    compute_zernike_coefficients,
+    compute_rms_metrics,
+    compute_aberration_magnitudes,
+    compute_rms_with_custom_filter,
+    zernike_noll,
+    zernike_label,
+    zernike_name,
+    noll_indices,
+    read_node_geometry,
+    find_geometry_file,
+    extract_displacements_by_subcase,
+    UNIT_TO_NM,
+    DEFAULT_N_MODES,
+    DEFAULT_FILTER_ORDERS,
+)
+
+from optimization_engine.extractors.extract_zernike_trajectory import (
+    ZernikeTrajectoryExtractor,
+    MODE_GROUPS,
+    MODE_NAMES,
+    compute_trajectory_params,
+)
+
+
+# ============================================================================
+# Configuration
+# ============================================================================
+
+N_MODES = 50
+FILTER_LOW_ORDERS = 4
+PLOT_DOWNSAMPLE = 12000
+COLORSCALE = 'Turbo'
+
+# Default design targets (nm)
+DEFAULT_TARGETS = {
+    'wfe_40_20': 4.0,
+    'wfe_60_20': 10.0,
+    'mfg_90': 20.0,
+}
+
+# Default annular aperture for M1 (271.5mm central hole diameter)
+DEFAULT_INNER_RADIUS = 135.75  # mm
+
+DISP_UNIT = 'mm'
+NM_SCALE = UNIT_TO_NM[DISP_UNIT]
+
+# Subcase mapping: subcase_id -> angle
+SUBCASE_ANGLE_MAP = {
+    '1': 90, '2': 20, '3': 40, '4': 60,
+    '90': 90, '20': 20, '40': 40, '60': 60,
+}
+
+
+# ============================================================================
+# Data Extraction Helpers
+# ============================================================================
+
+def load_study_params(db_path: str, trial_id: int = None) -> dict:
+    """Load design parameters from study database."""
+    import sqlite3
+    conn = sqlite3.connect(db_path)
+    c = conn.cursor()
+
+    if trial_id is None:
+        # Find best trial by weighted sum
+        c.execute('''
+            SELECT t.trial_id, tua.key, tua.value_json
+            FROM trials t
+            JOIN trial_user_attributes tua ON t.trial_id = tua.trial_id
+            WHERE t.state = 'COMPLETE' AND tua.key = 'weighted_sum'
+            ORDER BY CAST(tua.value_json AS REAL) ASC
+            LIMIT 1
+        ''')
+        row = c.fetchone()
+        if row:
+            trial_id = row[0]
+        else:
+            conn.close()
+            return {}
+
+    # Get all attributes for the trial
+    c.execute('''
+        SELECT key, value_json
+        FROM trial_user_attributes
+        WHERE trial_id = ?
+    ''', (trial_id,))
+    attrs = {row[0]: json.loads(row[1]) for row in c.fetchall()}
+
+    # Get parameters
+    c.execute('''
+        SELECT tp.key, tp.value_json
+        FROM trial_params tp
+        WHERE tp.trial_id = ?
+    ''', (trial_id,))
+    params = {row[0]: json.loads(row[1]) for row in c.fetchall()}
+
+    conn.close()
+
+    return {
+        'trial_id': trial_id,
+        'attributes': attrs,
+        'parameters': params,
+    }
+
+
+def build_wfe_arrays(node_ids, disp, node_geo):
+    """Build X, Y, WFE arrays from displacement data."""
+    X, Y, WFE = [], [], []
+    for nid, vec in zip(node_ids, disp):
+        geo = node_geo.get(int(nid))
+        if geo is None:
+            continue
+        X.append(geo[0])
+        Y.append(geo[1])
+        WFE.append(vec[2] * 2.0 * NM_SCALE)
+    return np.array(X), np.array(Y), np.array(WFE)
+
+
+def compute_relative_wfe(X1, Y1, WFE1, nids1, X2, Y2, WFE2, nids2):
+    """Compute WFE1 - WFE2 for common nodes."""
+    ref_map = {int(n): w for n, w in zip(nids2, WFE2)}
+    Xr, Yr, Wr = [], [], []
+    for nid, x, y, w in zip(nids1, X1, Y1, WFE1):
+        nid = int(nid)
+        if nid in ref_map:
+            Xr.append(x)
+            Yr.append(y)
+            Wr.append(w - ref_map[nid])
+    return np.array(Xr), np.array(Yr), np.array(Wr)
+
+
+def zernike_fit(X, Y, W, n_modes=N_MODES, inner_radius=None):
+    """Compute Zernike fit with optional annular masking."""
+    Xc = X - np.mean(X)
+    Yc = Y - np.mean(Y)
+    R_outer = float(np.max(np.hypot(Xc, Yc)))
+    r = np.hypot(Xc, Yc) / R_outer
+    th = np.arctan2(Yc, Xc)
+
+    # Annular mask
+    if inner_radius is not None:
+        r_inner_norm = inner_radius / R_outer
+        mask = (r >= r_inner_norm) & (r <= 1.0) & ~np.isnan(W)
+    else:
+        mask = (r <= 1.0) & ~np.isnan(W)
+
+    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 = W.astype(np.float64)
+
+    for start in range(0, len(idx), 100000):
+        sl = idx[start:start+100000]
+        Zb = np.column_stack([zernike_noll(j, r[sl].astype(np.float32), th[sl].astype(np.float32)).astype(np.float32)
+                              for j in range(1, m+1)])
+        G += (Zb.T @ Zb).astype(np.float64)
+        h += (Zb.T @ v[sl]).astype(np.float64)
+
+    try:
+        coeffs = np.linalg.solve(G, h)
+    except LinAlgError:
+        coeffs = np.linalg.lstsq(G, h, rcond=None)[0]
+
+    # Compute residuals
+    Z_all = np.column_stack([zernike_noll(j, r.astype(np.float32), th.astype(np.float32))
+                              for j in range(1, m+1)])
+    W_low4 = Z_all[:, :FILTER_LOW_ORDERS].dot(coeffs[:FILTER_LOW_ORDERS])
+    W_low3 = Z_all[:, :3].dot(coeffs[:3])
+    W_res_j4 = W - W_low4  # J1-J4 removed
+    W_res_j3 = W - W_low3  # J1-J3 removed
+
+    global_rms = float(np.sqrt(np.mean(W[mask]**2)))
+    filtered_rms = float(np.sqrt(np.mean(W_res_j4[mask]**2)))
+    rms_j1to3 = float(np.sqrt(np.mean(W_res_j3[mask]**2)))
+
+    return {
+        'coefficients': coeffs,
+        'R_outer': R_outer,
+        'global_rms': global_rms,
+        'filtered_rms': filtered_rms,
+        'rms_j1to3': rms_j1to3,
+        'W_res_filt': W_res_j4,
+        'mask': mask,
+        'n_masked': int(np.sum(mask)),
+        'n_total': len(W),
+    }
+
+
+def aberration_magnitudes(coeffs):
+    """Get individual aberration magnitudes from Zernike coefficients."""
+    defocus = float(abs(coeffs[3]))
+    astig = 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
+    return {
+        'defocus': defocus, 'astigmatism': astig, 'coma': coma,
+        'trefoil': trefoil, 'spherical': spherical,
+    }
+
+
+# ============================================================================
+# HTML Report Generation
+# ============================================================================
+
+def status_badge(value, target, unit='nm'):
+    """Return pass/fail badge HTML."""
+    if value <= target:
+        return f'<span class="badge pass">✅ {value:.2f} {unit} ≤ {target:.1f}</span>'
+    ratio = value / target
+    if ratio < 1.5:
+        return f'<span class="badge warn">⚠️ {value:.2f} {unit} ({ratio:.1f}× target)</span>'
+    return f'<span class="badge fail">❌ {value:.2f} {unit} ({ratio:.1f}× target)</span>'
+
+
+def make_surface_plot(X, Y, W_res, mask, inner_radius=None, title="", amp=0.5, downsample=PLOT_DOWNSAMPLE):
+    """Create a 3D surface plot of residual WFE."""
+    Xm, Ym, Wm = X[mask], Y[mask], W_res[mask]
+
+    n = len(Xm)
+    if n > downsample:
+        rng = np.random.default_rng(42)
+        sel = rng.choice(n, size=downsample, replace=False)
+        Xp, Yp, Wp = Xm[sel], Ym[sel], Wm[sel]
+    else:
+        Xp, Yp, Wp = Xm, Ym, Wm
+
+    res_amp = amp * Wp
+    max_amp = float(np.max(np.abs(res_amp))) if res_amp.size else 1.0
+
+    traces = []
+    try:
+        tri = Triangulation(Xp, Yp)
+        if tri.triangles is not None and len(tri.triangles) > 0:
+            # Filter triangles spanning central hole
+            if inner_radius is not None:
+                cx, cy = np.mean(X), np.mean(Y)
+                valid = []
+                for t in tri.triangles:
+                    vx = Xp[t] - cx
+                    vy = Yp[t] - cy
+                    vr = np.hypot(vx, vy)
+                    if np.any(vr < inner_radius * 0.9):
+                        continue
+                    p0, p1, p2 = Xp[t] + 1j*Yp[t]
+                    if max(abs(p1-p0), abs(p2-p1), abs(p0-p2)) > 2*inner_radius:
+                        continue
+                    valid.append(t)
+                if valid:
+                    tri_arr = np.array(valid)
+                else:
+                    tri_arr = tri.triangles
+            else:
+                tri_arr = tri.triangles
+
+            i, j, k = tri_arr.T
+            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="nm", side="right"), thickness=12, len=0.5, tickformat=".1f"),
+                hovertemplate="X: %{x:.1f}<br>Y: %{y:.1f}<br>Residual: %{z:.2f} nm<extra></extra>"
+            ))
+
+            # Inner hole circle
+            if inner_radius:
+                theta_c = np.linspace(0, 2*np.pi, 80)
+                traces.append(go.Scatter3d(
+                    x=cx + inner_radius*np.cos(theta_c),
+                    y=cy + inner_radius*np.sin(theta_c),
+                    z=np.zeros(80),
+                    mode='lines', line=dict(color='white', width=2),
+                    name='Central Hole', showlegend=False, hoverinfo='name'
+                ))
+    except Exception:
+        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
+        ))
+
+    fig = go.Figure(data=traces)
+    fig.update_layout(
+        scene=dict(
+            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*3.0, max_amp*3.0],
+                       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),
+        ),
+        margin=dict(t=30, b=10, l=10, r=10),
+        paper_bgcolor='rgba(0,0,0,0)',
+        plot_bgcolor='rgba(0,0,0,0)',
+        font=dict(color='#e0e0e0'),
+        height=500,
+        width=700,
+    )
+    return fig.to_html(include_plotlyjs=False, full_html=False, div_id=f"surface_{title.replace(' ','_')}")
+
+
+def make_bar_chart(coeffs, title="Zernike Coefficients", max_modes=30):
+    """Create horizontal bar chart of Zernike coefficient magnitudes."""
+    n = min(len(coeffs), max_modes)
+    labels = [zernike_label(j) for j in range(1, n+1)]
+    vals = np.abs(coeffs[:n])
+
+    fig = go.Figure(go.Bar(
+        x=vals, y=labels, orientation='h',
+        marker_color='#6366f1',
+        hovertemplate="%{y}<br>|c| = %{x:.3f} nm<extra></extra>",
+    ))
+    fig.update_layout(
+        height=max(400, n*22),
+        margin=dict(t=30, b=10, l=200, r=20),
+        paper_bgcolor='rgba(0,0,0,0)',
+        plot_bgcolor='rgba(17,24,39,0.8)',
+        font=dict(color='#e0e0e0', size=10),
+        xaxis=dict(title="|Coefficient| (nm)", gridcolor='rgba(128,128,128,0.2)'),
+        yaxis=dict(autorange='reversed'),
+    )
+    return fig.to_html(include_plotlyjs=False, full_html=False, div_id=f"bar_{title.replace(' ','_')}")
+
+
+def make_trajectory_plot(angles, coefficients_relative, mode_groups, sensitivity, title=""):
+    """Create trajectory visualization: Zernike modes vs elevation angle."""
+    fig = go.Figure()
+
+    # Plot each mode group
+    colors = ['#f59e0b', '#ef4444', '#10b981', '#6366f1', '#ec4899', '#14b8a6', '#f97316']
+    color_idx = 0
+
+    for group_name, noll_indices in mode_groups.items():
+        indices = [n - 5 for n in noll_indices if 5 <= n < 5 + coefficients_relative.shape[1]]
+        if not indices:
+            continue
+
+        # RSS of modes in this group at each angle
+        rss = np.sqrt(np.sum(coefficients_relative[:, indices]**2, axis=1))
+        color = colors[color_idx % len(colors)]
+
+        fig.add_trace(go.Scatter(
+            x=angles, y=rss,
+            mode='lines+markers',
+            name=MODE_NAMES.get(group_name, group_name),
+            line=dict(color=color, width=2),
+            marker=dict(size=8),
+            hovertemplate=f"{group_name}<br>%{{x}}°: %{{y:.2f}} nm<extra></extra>"
+        ))
+        color_idx += 1
+
+    fig.update_layout(
+        height=400,
+        margin=dict(t=30, b=40, l=60, r=20),
+        paper_bgcolor='rgba(0,0,0,0)',
+        plot_bgcolor='rgba(17,24,39,0.8)',
+        font=dict(color='#e0e0e0'),
+        xaxis=dict(title="Elevation Angle (°)", gridcolor='rgba(128,128,128,0.2)',
+                   tickvals=angles, dtick=10),
+        yaxis=dict(title="RMS (nm)", gridcolor='rgba(128,128,128,0.2)'),
+        legend=dict(x=0.02, y=0.98, bgcolor='rgba(17,24,39,0.7)'),
+    )
+    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="trajectory_plot")
+
+
+def make_sensitivity_bar(sensitivity_dict):
+    """Create stacked bar chart of axial vs lateral sensitivity per mode."""
+    modes = list(sensitivity_dict.keys())
+    axial = [sensitivity_dict[m]['axial'] for m in modes]
+    lateral = [sensitivity_dict[m]['lateral'] for m in modes]
+    labels = [MODE_NAMES.get(m, m) for m in modes]
+
+    fig = go.Figure()
+    fig.add_trace(go.Bar(
+        y=labels, x=axial, orientation='h',
+        name='Axial (sin θ)', marker_color='#f59e0b',
+        hovertemplate="%{y}<br>Axial: %{x:.3f} nm/unit<extra></extra>"
+    ))
+    fig.add_trace(go.Bar(
+        y=labels, x=lateral, orientation='h',
+        name='Lateral (cos θ)', marker_color='#6366f1',
+        hovertemplate="%{y}<br>Lateral: %{x:.3f} nm/unit<extra></extra>"
+    ))
+    fig.update_layout(
+        barmode='group',
+        height=max(300, len(modes)*40),
+        margin=dict(t=30, b=40, l=200, r=20),
+        paper_bgcolor='rgba(0,0,0,0)',
+        plot_bgcolor='rgba(17,24,39,0.8)',
+        font=dict(color='#e0e0e0', size=11),
+        xaxis=dict(title="Sensitivity (nm per load fraction)", gridcolor='rgba(128,128,128,0.2)'),
+        yaxis=dict(autorange='reversed'),
+        legend=dict(x=0.6, y=0.98, bgcolor='rgba(17,24,39,0.7)'),
+    )
+    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="sensitivity_bar")
+
+
+def make_per_angle_rms_plot(angle_rms_data, ref_angle=20):
+    """Create bar chart of per-angle RMS relative to reference."""
+    angles = sorted(angle_rms_data.keys())
+    rms_vals = [angle_rms_data[a] for a in angles]
+    labels = [f"{a}° vs {ref_angle}°" for a in angles]
+
+    fig = go.Figure(go.Bar(
+        x=labels, y=rms_vals,
+        marker_color=['#10b981' if v < 10 else '#f59e0b' if v < 20 else '#ef4444' for v in rms_vals],
+        text=[f"{v:.2f} nm" for v in rms_vals],
+        textposition='outside',
+        hovertemplate="%{x}: %{y:.2f} nm<extra></extra>"
+    ))
+    fig.update_layout(
+        height=350,
+        margin=dict(t=30, b=40, l=60, r=20),
+        paper_bgcolor='rgba(0,0,0,0)',
+        plot_bgcolor='rgba(17,24,39,0.8)',
+        font=dict(color='#e0e0e0'),
+        yaxis=dict(title="Filtered RMS WFE (nm)", gridcolor='rgba(128,128,128,0.2)'),
+    )
+    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="per_angle_rms")
+
+
+# ============================================================================
+# Main Report Builder
+# ============================================================================
+
+def generate_report(
+    op2_path: Path,
+    inner_radius: float = None,
+    targets: dict = None,
+    study_db: str = None,
+    trial_id: int = None,
+    title: str = "M1 Mirror Optical Performance Report",
+    study_name: str = None,
+) -> Path:
+    """Generate comprehensive optical performance HTML report."""
+
+    targets = targets or DEFAULT_TARGETS
+    timestamp = datetime.now().strftime("%Y-%m-%d %H:%M")
+    ts_file = datetime.now().strftime("%Y%m%d_%H%M%S")
+
+    print("=" * 70)
+    print("  ATOMIZER OPTICAL PERFORMANCE REPORT GENERATOR")
+    print("=" * 70)
+    print(f"\nOP2 File: {op2_path.name}")
+    print(f"Inner Radius: {inner_radius} mm" if inner_radius else "Aperture: Full disk")
+
+    # ------------------------------------------------------------------
+    # 1. Load geometry & displacement data
+    # ------------------------------------------------------------------
+    print("\n[1/5] Loading data...")
+    geo_path = find_geometry_file(op2_path)
+    node_geo = read_node_geometry(geo_path)
+    print(f"  Geometry: {geo_path.name} ({len(node_geo)} nodes)")
+
+    op2 = OP2()
+    op2.read_op2(str(op2_path))
+    displacements = extract_displacements_by_subcase(op2_path)
+    print(f"  Subcases: {list(displacements.keys())}")
+
+    # Map subcases to angles
+    subcase_map = {}
+    for angle in ['90', '20', '40', '60']:
+        if angle in displacements:
+            subcase_map[angle] = angle
+    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"  Subcase map: {subcase_map}")
+
+    # Also detect intermediate angles (30, 50) if present
+    extra_angles = []
+    for a in ['30', '50']:
+        if a in displacements:
+            extra_angles.append(a)
+    if extra_angles:
+        print(f"  Extra angles detected: {extra_angles}")
+
+    # ------------------------------------------------------------------
+    # 2. Per-angle Zernike analysis
+    # ------------------------------------------------------------------
+    print("\n[2/5] Per-angle Zernike analysis...")
+
+    ref_label = subcase_map['20']
+    ref_data = displacements[ref_label]
+    X_ref, Y_ref, WFE_ref = build_wfe_arrays(ref_data['node_ids'], ref_data['disp'], node_geo)
+
+    # Analysis results storage
+    angle_results = {}  # angle -> {rms_data, X, Y, WFE, ...}
+
+    for angle_name, label in subcase_map.items():
+        data = displacements[label]
+        X, Y, WFE = build_wfe_arrays(data['node_ids'], data['disp'], node_geo)
+
+        # Absolute fit
+        rms_abs = zernike_fit(X, Y, WFE, inner_radius=inner_radius)
+
+        # Relative fit (vs 20 deg reference)
+        if angle_name != '20':
+            Xr, Yr, Wr = compute_relative_wfe(
+                X, Y, WFE, data['node_ids'],
+                X_ref, Y_ref, WFE_ref, ref_data['node_ids']
+            )
+            rms_rel = zernike_fit(Xr, Yr, Wr, inner_radius=inner_radius)
+        else:
+            Xr, Yr, Wr = X, Y, np.zeros_like(WFE)
+            rms_rel = {'filtered_rms': 0.0, 'rms_j1to3': 0.0, 'coefficients': np.zeros(N_MODES)}
+
+        angle_results[int(angle_name)] = {
+            'X': X, 'Y': Y, 'WFE': WFE,
+            'X_rel': Xr, 'Y_rel': Yr, 'WFE_rel': Wr,
+            'rms_abs': rms_abs,
+            'rms_rel': rms_rel,
+            'aberrations_abs': aberration_magnitudes(rms_abs['coefficients']),
+            'aberrations_rel': aberration_magnitudes(rms_rel['coefficients']) if angle_name != '20' else None,
+        }
+        print(f"  {angle_name}° - Abs Filt: {rms_abs['filtered_rms']:.2f} nm, "
+              f"Rel Filt: {rms_rel['filtered_rms']:.2f} nm")
+
+    # Extra angles (30, 50)
+    for ea in extra_angles:
+        data = displacements[ea]
+        X, Y, WFE = build_wfe_arrays(data['node_ids'], data['disp'], node_geo)
+        Xr, Yr, Wr = compute_relative_wfe(
+            X, Y, WFE, data['node_ids'],
+            X_ref, Y_ref, WFE_ref, ref_data['node_ids']
+        )
+        rms_abs = zernike_fit(X, Y, WFE, inner_radius=inner_radius)
+        rms_rel = zernike_fit(Xr, Yr, Wr, inner_radius=inner_radius)
+        angle_results[int(ea)] = {
+            'X': X, 'Y': Y, 'WFE': WFE,
+            'X_rel': Xr, 'Y_rel': Yr, 'WFE_rel': Wr,
+            'rms_abs': rms_abs,
+            'rms_rel': rms_rel,
+            'aberrations_abs': aberration_magnitudes(rms_abs['coefficients']),
+            'aberrations_rel': aberration_magnitudes(rms_rel['coefficients']),
+        }
+        print(f"  {ea}° - Abs Filt: {rms_abs['filtered_rms']:.2f} nm, "
+              f"Rel Filt: {rms_rel['filtered_rms']:.2f} nm")
+
+    # ------------------------------------------------------------------
+    # 3. Trajectory analysis
+    # ------------------------------------------------------------------
+    print("\n[3/5] Trajectory analysis...")
+    traj_result = None
+    try:
+        traj_extractor = ZernikeTrajectoryExtractor(
+            op2_file=op2_path,
+            bdf_file=geo_path,
+            reference_angle=20.0,
+            unit=DISP_UNIT,
+            n_modes=N_MODES,
+            inner_radius=inner_radius,
+        )
+        traj_result = traj_extractor.extract_trajectory(exclude_angles=[90.0])
+        print(f"  R² fit: {traj_result['linear_fit_r2']:.4f}")
+        print(f"  Dominant mode: {traj_result['dominant_mode']}")
+        print(f"  Total filtered RMS: {traj_result['total_filtered_rms_nm']:.2f} nm")
+    except Exception as e:
+        print(f"  [WARN] Trajectory analysis failed: {e}")
+
+    # ------------------------------------------------------------------
+    # 4. Manufacturing analysis
+    # ------------------------------------------------------------------
+    print("\n[4/5] Manufacturing analysis...")
+    r90 = angle_results[90]
+    mfg_abs_aberr = r90['aberrations_abs']
+    mfg_correction = aberration_magnitudes(angle_results[90]['rms_rel']['coefficients'])
+    mfg_rms_j1to3 = r90['rms_rel']['rms_j1to3']
+    print(f"  MFG 90 (J1-J3 filtered): {mfg_rms_j1to3:.2f} nm")
+    print(f"  Correction - Astigmatism: {mfg_correction['astigmatism']:.2f} nm, "
+          f"Coma: {mfg_correction['coma']:.2f} nm")
+
+    # ------------------------------------------------------------------
+    # 5. Load study params (optional)
+    # ------------------------------------------------------------------
+    study_params = None
+    if study_db:
+        print("\n[5/5] Loading study parameters...")
+        try:
+            study_params = load_study_params(study_db, trial_id)
+            print(f"  Trial #{study_params.get('trial_id', '?')}")
+        except Exception as e:
+            print(f"  [WARN] Could not load study params: {e}")
+    else:
+        print("\n[5/5] No study database provided (skipping design parameters)")
+
+    # ------------------------------------------------------------------
+    # Key metrics for executive summary
+    # ------------------------------------------------------------------
+    wfe_40_20 = angle_results[40]['rms_rel']['filtered_rms']
+    wfe_60_20 = angle_results[60]['rms_rel']['filtered_rms']
+    mfg_90 = mfg_rms_j1to3
+
+    # Weighted sum
+    ws = 6*wfe_40_20 + 5*wfe_60_20 + 3*mfg_90
+
+    # ------------------------------------------------------------------
+    # Generate HTML
+    # ------------------------------------------------------------------
+    print("\nGenerating HTML report...")
+
+    # Surface plots
+    surf_40 = make_surface_plot(
+        angle_results[40]['X_rel'], angle_results[40]['Y_rel'],
+        angle_results[40]['rms_rel']['W_res_filt'], angle_results[40]['rms_rel']['mask'],
+        inner_radius=inner_radius, title="40 vs 20"
+    )
+    surf_60 = make_surface_plot(
+        angle_results[60]['X_rel'], angle_results[60]['Y_rel'],
+        angle_results[60]['rms_rel']['W_res_filt'], angle_results[60]['rms_rel']['mask'],
+        inner_radius=inner_radius, title="60 vs 20"
+    )
+    surf_90 = make_surface_plot(
+        angle_results[90]['X'], angle_results[90]['Y'],
+        angle_results[90]['rms_abs']['W_res_filt'], angle_results[90]['rms_abs']['mask'],
+        inner_radius=inner_radius, title="90 abs"
+    )
+
+    # Bar charts
+    bar_40 = make_bar_chart(angle_results[40]['rms_rel']['coefficients'], title="40v20 coeffs")
+    bar_60 = make_bar_chart(angle_results[60]['rms_rel']['coefficients'], title="60v20 coeffs")
+    bar_90 = make_bar_chart(angle_results[90]['rms_abs']['coefficients'], title="90abs coeffs")
+
+    # Per-angle RMS plot
+    angle_rms_data = {}
+    for ang in sorted(angle_results.keys()):
+        if ang != 20:
+            angle_rms_data[ang] = angle_results[ang]['rms_rel']['filtered_rms']
+    per_angle_plot = make_per_angle_rms_plot(angle_rms_data)
+
+    # Trajectory & sensitivity plots
+    traj_plot_html = ""
+    sens_plot_html = ""
+    if traj_result:
+        coeffs_rel = np.array(traj_result['coefficients_relative'])
+        traj_plot_html = make_trajectory_plot(
+            traj_result['angles_deg'], coeffs_rel, MODE_GROUPS,
+            traj_result['sensitivity_matrix']
+        )
+        sens_plot_html = make_sensitivity_bar(traj_result['sensitivity_matrix'])
+
+    # Design parameters table
+    params_html = ""
+    if study_params and study_params.get('parameters'):
+        params = study_params['parameters']
+        rows = ""
+        for k, v in sorted(params.items()):
+            unit = "°" if "angle" in k else "mm"
+            rows += f"<tr><td>{k}</td><td>{v:.4f} {unit}</td></tr>\n"
+        params_html = f"""
+        <div class="section">
+            <h2>🔧 Design Parameters (Trial #{study_params.get('trial_id', '?')})</h2>
+            <table class="data-table"><thead><tr><th>Parameter</th><th>Value</th></tr></thead>
+            <tbody>{rows}</tbody></table>
+        </div>
+        """
+
+    # Per-angle detail table
+    angle_detail_rows = ""
+    for ang in sorted(angle_results.keys()):
+        r = angle_results[ang]
+        rel_filt = r['rms_rel']['filtered_rms']
+        abs_filt = r['rms_abs']['filtered_rms']
+        abs_glob = r['rms_abs']['global_rms']
+        ab = r['aberrations_abs']
+        angle_detail_rows += f"""<tr>
+            <td><b>{ang}°</b></td>
+            <td>{abs_glob:.2f}</td><td>{abs_filt:.2f}</td>
+            <td>{rel_filt:.2f}</td>
+            <td>{ab['astigmatism']:.2f}</td><td>{ab['coma']:.2f}</td>
+            <td>{ab['trefoil']:.2f}</td><td>{ab['spherical']:.2f}</td>
+        </tr>"""
+
+    # Trajectory metrics table
+    traj_metrics_html = ""
+    if traj_result:
+        traj_metrics_html = f"""
+        <div class="metrics-grid">
+            <div class="metric-card">
+                <div class="metric-label">Coma RMS</div>
+                <div class="metric-value">{traj_result['coma_rms_nm']:.2f} nm</div>
+            </div>
+            <div class="metric-card">
+                <div class="metric-label">Astigmatism RMS</div>
+                <div class="metric-value">{traj_result['astigmatism_rms_nm']:.2f} nm</div>
+            </div>
+            <div class="metric-card">
+                <div class="metric-label">Trefoil RMS</div>
+                <div class="metric-value">{traj_result['trefoil_rms_nm']:.2f} nm</div>
+            </div>
+            <div class="metric-card">
+                <div class="metric-label">Spherical RMS</div>
+                <div class="metric-value">{traj_result['spherical_rms_nm']:.2f} nm</div>
+            </div>
+            <div class="metric-card">
+                <div class="metric-label">Total Filtered RMS</div>
+                <div class="metric-value">{traj_result['total_filtered_rms_nm']:.2f} nm</div>
+            </div>
+            <div class="metric-card">
+                <div class="metric-label">Linear Fit R²</div>
+                <div class="metric-value">{traj_result['linear_fit_r2']:.4f}</div>
+            </div>
+        </div>
+        <p class="note">Dominant aberration mode: <b>{MODE_NAMES.get(traj_result['dominant_mode'], traj_result['dominant_mode'])}</b></p>
+        <p class="note">Mode ranking: {' → '.join(traj_result['mode_ranking'][:5])}</p>
+        """
+
+    # Manufacturing details
+    mfg_html = f"""
+    <table class="data-table">
+        <thead><tr><th>Metric</th><th>Absolute 90°</th><th>Correction (90°−20°)</th></tr></thead>
+        <tbody>
+            <tr><td>Defocus (J4)</td><td>{mfg_abs_aberr['defocus']:.2f} nm</td><td>{mfg_correction['defocus']:.2f} nm</td></tr>
+            <tr><td>Astigmatism (J5+J6)</td><td>{mfg_abs_aberr['astigmatism']:.2f} nm</td><td>{mfg_correction['astigmatism']:.2f} nm</td></tr>
+            <tr><td>Coma (J7+J8)</td><td>{mfg_abs_aberr['coma']:.2f} nm</td><td>{mfg_correction['coma']:.2f} nm</td></tr>
+            <tr><td>Trefoil (J9+J10)</td><td>{mfg_abs_aberr['trefoil']:.2f} nm</td><td>{mfg_correction['trefoil']:.2f} nm</td></tr>
+            <tr><td>Spherical (J11)</td><td>{mfg_abs_aberr['spherical']:.2f} nm</td><td>{mfg_correction['spherical']:.2f} nm</td></tr>
+            <tr class="highlight"><td><b>J1−J3 Filtered RMS</b></td><td>{r90['rms_abs']['rms_j1to3']:.2f} nm</td><td><b>{mfg_rms_j1to3:.2f} nm</b></td></tr>
+        </tbody>
+    </table>
+    """
+
+    # Assemble full HTML
+    html = f"""<!DOCTYPE html>
+<html lang="en">
+<head>
+<meta charset="UTF-8">
+<meta name="viewport" content="width=device-width, initial-scale=1.0">
+<title>{title}</title>
+<script src="https://cdn.plot.ly/plotly-2.27.0.min.js"></script>
+<style>
+    :root {{
+        --bg-primary: #0f172a;
+        --bg-secondary: #1e293b;
+        --bg-card: #1e293b;
+        --text-primary: #f1f5f9;
+        --text-secondary: #94a3b8;
+        --accent: #6366f1;
+        --accent-hover: #818cf8;
+        --success: #10b981;
+        --warning: #f59e0b;
+        --danger: #ef4444;
+        --border: #334155;
+    }}
+    * {{ margin: 0; padding: 0; box-sizing: border-box; }}
+    body {{
+        font-family: 'Inter', -apple-system, BlinkMacSystemFont, 'Segoe UI', system-ui, sans-serif;
+        background: var(--bg-primary);
+        color: var(--text-primary);
+        line-height: 1.6;
+    }}
+    .container {{ max-width: 1400px; margin: 0 auto; padding: 2rem; }}
+
+    /* Header */
+    .header {{
+        background: linear-gradient(135deg, #1e293b 0%, #0f172a 100%);
+        border: 1px solid var(--border);
+        border-radius: 12px;
+        padding: 2rem 3rem;
+        margin-bottom: 2rem;
+        display: flex;
+        justify-content: space-between;
+        align-items: center;
+    }}
+    .header h1 {{ font-size: 1.8rem; font-weight: 700; }}
+    .header .subtitle {{ color: var(--text-secondary); font-size: 0.95rem; margin-top: 0.3rem; }}
+    .header .branding {{
+        text-align: right;
+        font-size: 0.85rem;
+        color: var(--text-secondary);
+    }}
+    .header .branding .logo {{ font-size: 1.4rem; font-weight: 700; color: var(--accent); }}
+
+    /* Sections */
+    .section {{
+        background: var(--bg-card);
+        border: 1px solid var(--border);
+        border-radius: 12px;
+        padding: 1.5rem 2rem;
+        margin-bottom: 1.5rem;
+    }}
+    .section h2 {{
+        font-size: 1.3rem;
+        margin-bottom: 1rem;
+        padding-bottom: 0.5rem;
+        border-bottom: 1px solid var(--border);
+    }}
+
+    /* Executive Summary */
+    .exec-grid {{
+        display: grid;
+        grid-template-columns: repeat(auto-fit, minmax(300px, 1fr));
+        gap: 1rem;
+        margin: 1rem 0;
+    }}
+    .exec-card {{
+        background: var(--bg-primary);
+        border: 1px solid var(--border);
+        border-radius: 8px;
+        padding: 1.2rem;
+    }}
+    .exec-card .label {{ font-size: 0.85rem; color: var(--text-secondary); margin-bottom: 0.3rem; }}
+    .exec-card .value {{ font-size: 1.8rem; font-weight: 700; }}
+    .exec-card .target {{ font-size: 0.8rem; color: var(--text-secondary); margin-top: 0.3rem; }}
+
+    /* Badges */
+    .badge {{
+        display: inline-block;
+        padding: 0.2rem 0.6rem;
+        border-radius: 4px;
+        font-size: 0.85rem;
+        font-weight: 600;
+    }}
+    .badge.pass {{ background: rgba(16,185,129,0.15); color: var(--success); }}
+    .badge.warn {{ background: rgba(245,158,11,0.15); color: var(--warning); }}
+    .badge.fail {{ background: rgba(239,68,68,0.15); color: var(--danger); }}
+
+    /* Tables */
+    .data-table {{
+        width: 100%;
+        border-collapse: collapse;
+        margin: 0.5rem 0;
+    }}
+    .data-table th, .data-table td {{
+        padding: 0.6rem 1rem;
+        text-align: left;
+        border-bottom: 1px solid var(--border);
+    }}
+    .data-table th {{
+        background: var(--bg-primary);
+        font-weight: 600;
+        font-size: 0.85rem;
+        text-transform: uppercase;
+        letter-spacing: 0.05em;
+        color: var(--text-secondary);
+    }}
+    .data-table tr.highlight td {{
+        background: rgba(99,102,241,0.08);
+        font-weight: 600;
+    }}
+
+    /* Metrics Grid */
+    .metrics-grid {{
+        display: grid;
+        grid-template-columns: repeat(auto-fit, minmax(180px, 1fr));
+        gap: 0.8rem;
+        margin: 1rem 0;
+    }}
+    .metric-card {{
+        background: var(--bg-primary);
+        border: 1px solid var(--border);
+        border-radius: 8px;
+        padding: 1rem;
+        text-align: center;
+    }}
+    .metric-label {{ font-size: 0.8rem; color: var(--text-secondary); margin-bottom: 0.3rem; }}
+    .metric-value {{ font-size: 1.3rem; font-weight: 700; color: var(--accent); }}
+
+    /* Plots */
+    .plot-grid {{
+        display: grid;
+        grid-template-columns: repeat(auto-fit, minmax(650px, 1fr));
+        gap: 1rem;
+    }}
+    .plot-container {{
+        background: var(--bg-primary);
+        border: 1px solid var(--border);
+        border-radius: 8px;
+        padding: 1rem;
+    }}
+    .plot-container h3 {{
+        font-size: 1rem;
+        margin-bottom: 0.5rem;
+        color: var(--text-secondary);
+    }}
+
+    /* Collapsible */
+    details {{ margin: 0.5rem 0; }}
+    summary {{
+        cursor: pointer;
+        font-weight: 600;
+        padding: 0.5rem;
+        background: var(--bg-primary);
+        border-radius: 6px;
+        border: 1px solid var(--border);
+    }}
+    summary:hover {{ background: rgba(99,102,241,0.1); }}
+    details > div {{ padding: 1rem; }}
+
+    .note {{ color: var(--text-secondary); font-size: 0.9rem; margin: 0.5rem 0; }}
+
+    /* Print styles */
+    @media print {{
+        body {{ background: white; color: black; }}
+        .section {{ border: 1px solid #ccc; page-break-inside: avoid; }}
+    }}
+
+    .two-col {{ display: grid; grid-template-columns: 1fr 1fr; gap: 1.5rem; }}
+    @media (max-width: 900px) {{ .two-col {{ grid-template-columns: 1fr; }} }}
+</style>
+</head>
+<body>
+<div class="container">
+
+<!-- Header -->
+<div class="header">
+    <div>
+        <h1>🔭 {title}</h1>
+        <div class="subtitle">Generated {timestamp} &nbsp;|&nbsp; OP2: {op2_path.name}</div>
+        {'<div class="subtitle">Study: ' + study_name + '</div>' if study_name else ''}
+    </div>
+    <div class="branding">
+        <div class="logo">ATOMIZER</div>
+        <div>by Atomaste</div>
+        <div style="margin-top:0.3rem">FEA Optimization Platform</div>
+    </div>
+</div>
+
+<!-- Executive Summary -->
+<div class="section">
+    <h2>📋 Executive Summary</h2>
+    <div class="exec-grid">
+        <div class="exec-card">
+            <div class="label">WFE 40° vs 20° (Tracking)</div>
+            <div class="value">{wfe_40_20:.2f} <small>nm</small></div>
+            <div class="target">{status_badge(wfe_40_20, targets['wfe_40_20'])}</div>
+        </div>
+        <div class="exec-card">
+            <div class="label">WFE 60° vs 20° (Tracking)</div>
+            <div class="value">{wfe_60_20:.2f} <small>nm</small></div>
+            <div class="target">{status_badge(wfe_60_20, targets['wfe_60_20'])}</div>
+        </div>
+        <div class="exec-card">
+            <div class="label">MFG 90° (J1−J3 Filtered)</div>
+            <div class="value">{mfg_90:.2f} <small>nm</small></div>
+            <div class="target">{status_badge(mfg_90, targets['mfg_90'])}</div>
+        </div>
+        <div class="exec-card">
+            <div class="label">Weighted Sum (6·W40 + 5·W60 + 3·MFG)</div>
+            <div class="value" style="color: var(--accent)">{ws:.1f}</div>
+            <div class="target">Lower is better</div>
+        </div>
+    </div>
+    {'<p class="note">Annular aperture: inner radius = ' + f'{inner_radius:.1f} mm (ø{2*inner_radius:.1f} mm central hole)' + '</p>' if inner_radius else ''}
+</div>
+
+<!-- Per-Angle Summary -->
+<div class="section">
+    <h2>📊 Per-Angle RMS Summary</h2>
+    {per_angle_plot}
+    <table class="data-table" style="margin-top:1rem">
+        <thead>
+            <tr>
+                <th>Angle</th><th>Abs Global RMS</th><th>Abs Filtered RMS</th>
+                <th>Rel Filtered RMS</th>
+                <th>Astigmatism</th><th>Coma</th><th>Trefoil</th><th>Spherical</th>
+            </tr>
+        </thead>
+        <tbody>{angle_detail_rows}</tbody>
+    </table>
+    <p class="note">All values in nm. Filtered = J1−J4 removed. Relative = vs 20° reference. Aberrations are absolute.</p>
+</div>
+
+<!-- Surface Plots -->
+<div class="section">
+    <h2>🌊 Wavefront Error Surface Maps</h2>
+    <p class="note">3D residual surfaces after removing piston, tip, tilt, and defocus (J1−J4). Interactive — drag to rotate.</p>
+    <div class="plot-grid">
+        <div class="plot-container">
+            <h3>40° vs 20° (Relative)</h3>
+            {surf_40}
+        </div>
+        <div class="plot-container">
+            <h3>60° vs 20° (Relative)</h3>
+            {surf_60}
+        </div>
+    </div>
+    <div class="plot-container" style="margin-top:1rem">
+        <h3>90° Manufacturing (Absolute)</h3>
+        {surf_90}
+    </div>
+</div>
+
+<!-- Trajectory Analysis -->
+{'<div class="section"><h2>📈 Zernike Trajectory Analysis</h2>' +
+ '<p class="note">Mode-specific integrated RMS across the operating elevation range. ' +
+ 'The linear model c<sub>j</sub>(θ) = a<sub>j</sub>·Δsinθ + b<sub>j</sub>·Δcosθ decomposes gravity into axial and lateral components.</p>' +
+ traj_metrics_html +
+ '<div class="two-col" style="margin-top:1rem">' +
+ '<div class="plot-container"><h3>Mode RMS vs Elevation Angle</h3>' + traj_plot_html + '</div>' +
+ '<div class="plot-container"><h3>Axial vs Lateral Sensitivity</h3>' + sens_plot_html + '</div>' +
+ '</div></div>' if traj_result else ''}
+
+<!-- Manufacturing Analysis -->
+<div class="section">
+    <h2>🏭 Manufacturing Analysis (90° Orientation)</h2>
+    <p class="note">
+        The mirror is manufactured (polished) at 90° orientation. The "Correction" column shows the
+        aberrations that must be polished out to achieve the 20° operational figure.
+    </p>
+    {mfg_html}
+</div>
+
+<!-- Design Parameters -->
+{params_html}
+
+<!-- Zernike Coefficient Details -->
+<div class="section">
+    <h2>🔬 Zernike Coefficient Details</h2>
+    <details>
+        <summary>40° vs 20° — Relative Coefficients</summary>
+        <div>{bar_40}</div>
+    </details>
+    <details>
+        <summary>60° vs 20° — Relative Coefficients</summary>
+        <div>{bar_60}</div>
+    </details>
+    <details>
+        <summary>90° — Absolute Coefficients</summary>
+        <div>{bar_90}</div>
+    </details>
+</div>
+
+<!-- Methodology -->
+<div class="section" style="opacity:0.85">
+    <h2>📝 Methodology</h2>
+    <table class="data-table">
+        <tbody>
+            <tr><td><b>Zernike Modes</b></td><td>{N_MODES} (Noll convention)</td></tr>
+            <tr><td><b>Filtered Modes</b></td><td>J1−J4 (Piston, Tip, Tilt, Defocus)</td></tr>
+            <tr><td><b>WFE Calculation</b></td><td>WFE = 2 × Surface Error (reflective)</td></tr>
+            <tr><td><b>Displacement Unit</b></td><td>{DISP_UNIT} → nm ({NM_SCALE:.0e}×)</td></tr>
+            <tr><td><b>Aperture</b></td><td>{'Annular (inner R = ' + f'{inner_radius:.1f} mm)' if inner_radius else 'Full disk'}</td></tr>
+            <tr><td><b>Reference Angle</b></td><td>20° (polishing/measurement orientation)</td></tr>
+            <tr><td><b>MFG Objective</b></td><td>90°−20° relative, J1−J3 filtered (optician workload)</td></tr>
+            <tr><td><b>Weighted Sum</b></td><td>6×WFE(40−20) + 5×WFE(60−20) + 3×MFG(90)</td></tr>
+            {'<tr><td><b>Trajectory R²</b></td><td>' + f'{traj_result["linear_fit_r2"]:.6f}' + '</td></tr>' if traj_result else ''}
+        </tbody>
+    </table>
+</div>
+
+<!-- Footer -->
+<div style="text-align:center; padding:2rem; color:var(--text-secondary); font-size:0.8rem;">
+    Generated by <b>Atomizer</b> Optical Report Generator &nbsp;|&nbsp; {timestamp}<br>
+    © Atomaste &nbsp;|&nbsp; atomaste.ca
+</div>
+
+</div>
+</body>
+</html>"""
+
+    # Write output
+    output_path = op2_path.parent / f"{op2_path.stem}_OPTICAL_REPORT_{ts_file}.html"
+    output_path.write_text(html, encoding='utf-8')
+
+    print(f"\n{'=' * 70}")
+    print(f"REPORT GENERATED: {output_path.name}")
+    print(f"{'=' * 70}")
+    print(f"\nLocation: {output_path}")
+    print(f"Size: {output_path.stat().st_size / 1024:.0f} KB")
+
+    return output_path
+
+
+# ============================================================================
+# CLI
+# ============================================================================
+
+def main():
+    parser = argparse.ArgumentParser(
+        description='Atomizer Optical Performance Report Generator',
+        epilog='Generates a comprehensive CDR-ready HTML report from FEA results.'
+    )
+    parser.add_argument('op2_file', nargs='?', help='Path to OP2 results file')
+    parser.add_argument('--inner-radius', '-r', type=float, default=None,
+                       help=f'Inner radius of central hole in mm (default: {DEFAULT_INNER_RADIUS}mm for M1)')
+    parser.add_argument('--inner-diameter', '-d', type=float, default=None,
+                       help='Inner diameter of central hole in mm')
+    parser.add_argument('--no-annular', action='store_true',
+                       help='Disable annular aperture (treat as full disk)')
+    parser.add_argument('--target-40', type=float, default=DEFAULT_TARGETS['wfe_40_20'],
+                       help=f'WFE 40-20 target in nm (default: {DEFAULT_TARGETS["wfe_40_20"]})')
+    parser.add_argument('--target-60', type=float, default=DEFAULT_TARGETS['wfe_60_20'],
+                       help=f'WFE 60-20 target in nm (default: {DEFAULT_TARGETS["wfe_60_20"]})')
+    parser.add_argument('--target-mfg', type=float, default=DEFAULT_TARGETS['mfg_90'],
+                       help=f'MFG 90 target in nm (default: {DEFAULT_TARGETS["mfg_90"]})')
+    parser.add_argument('--study-db', type=str, default=None,
+                       help='Path to study.db for design parameters')
+    parser.add_argument('--trial', type=int, default=None,
+                       help='Trial ID (default: best trial)')
+    parser.add_argument('--title', type=str, default="M1 Mirror Optical Performance Report",
+                       help='Report title')
+    parser.add_argument('--study-name', type=str, default=None,
+                       help='Study name for report header')
+
+    args = parser.parse_args()
+
+    # 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:
+        # Search current directory
+        cwd = Path.cwd()
+        candidates = list(cwd.glob("*solution*.op2")) + list(cwd.glob("*.op2"))
+        if not candidates:
+            print("ERROR: No OP2 file found. Specify path as argument.")
+            sys.exit(1)
+        op2_path = max(candidates, key=lambda p: p.stat().st_mtime)
+        print(f"Found: {op2_path}")
+
+    # Handle inner radius
+    inner_radius = DEFAULT_INNER_RADIUS  # Default to M1 annular
+    if args.no_annular:
+        inner_radius = None
+    elif args.inner_diameter is not None:
+        inner_radius = args.inner_diameter / 2.0
+    elif args.inner_radius is not None:
+        inner_radius = args.inner_radius
+
+    targets = {
+        'wfe_40_20': args.target_40,
+        'wfe_60_20': args.target_60,
+        'mfg_90': args.target_mfg,
+    }
+
+    try:
+        generate_report(
+            op2_path=op2_path,
+            inner_radius=inner_radius,
+            targets=targets,
+            study_db=args.study_db,
+            trial_id=args.trial,
+            title=args.title,
+            study_name=args.study_name,
+        )
+    except Exception as e:
+        print(f"\nERROR: {e}")
+        import traceback
+        traceback.print_exc()
+        sys.exit(1)
+
+
+if __name__ == '__main__':
+    main()