Atomizer/tools/generate_optical_report.py

#!/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
from scipy.interpolate import griddata
from scipy.fft import fft2, fftfreq

# 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,
}

# ── Professional light-theme Plotly layout defaults ──
_PLOTLY_LIGHT_LAYOUT = dict(
    paper_bgcolor='#ffffff',
    plot_bgcolor='#f8fafc',
    font=dict(family="Inter, system-ui, sans-serif", color='#1e293b', size=12),
)

# Professional blue palette for charts
_BLUE_PALETTE = ['#2563eb', '#3b82f6', '#60a5fa', '#93c5fd', '#1d4ed8', '#1e40af']
_MODE_COLORS = ['#2563eb', '#dc2626', '#16a34a', '#9333ea', '#ea580c', '#0891b2', '#4f46e5']
_MODE_DASHES = ['solid', 'dash', 'dot', 'dashdot', 'longdash', 'longdashdot', 'solid']

# High-resolution PNG export settings
PNG_EXPORT_SCALE = 4  # 4x resolution (e.g., 700px width -> 2800px export)
PNG_EXPORT_FORMAT = 'png'

def get_plotly_config(filename_prefix="plot"):
    """Get Plotly config with high-resolution PNG export settings."""
    return {
        'toImageButtonOptions': {
            'format': PNG_EXPORT_FORMAT,
            'filename': filename_prefix,
            'height': None,  # Use current height
            'width': None,   # Use current width
            'scale': PNG_EXPORT_SCALE,  # 4x resolution multiplier
        },
        'displaylogo': False,
        'modeBarButtonsToAdd': ['hoverClosest3d'],
    }


# ============================================================================
# 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,
    }


def compute_surface_psd(X, Y, Z, aperture_radius):
    """Compute PSD and band RMS of surface height data (Tony Hull methodology).

    Interpolates scattered FEA data onto a uniform grid, applies Hann window,
    computes 2D FFT, and radially averages. Band RMS uses direct Parseval
    summation for correct dimensional results.

    Args:
        X, Y: coordinates [mm]
        Z: surface height [nm]
        aperture_radius: mirror radius [mm]

    Returns:
        dict with 'freqs', 'psd' (for plotting), 'bands' (gravity/support/hf/total RMS in nm)
    """
    Xc = X - np.mean(X)
    Yc = Y - np.mean(Y)

    N = 256 if len(X) >= 1000 else min(128, int(np.sqrt(len(X))))

    x_range = np.linspace(Xc.min(), Xc.max(), N)
    y_range = np.linspace(Yc.min(), Yc.max(), N)
    X_grid, Y_grid = np.meshgrid(x_range, y_range)

    Z_grid = griddata((Xc, Yc), Z, (X_grid, Y_grid), method='cubic')
    Z_grid = np.nan_to_num(Z_grid, nan=0.0)

    # Circular aperture mask
    R_grid = np.sqrt(X_grid**2 + Y_grid**2)
    apt_mask = R_grid <= aperture_radius
    Z_grid[~apt_mask] = 0.0

    # Count valid aperture pixels for RMS normalization
    n_apt = int(np.sum(apt_mask))
    if n_apt < 10:
        raise ValueError("Too few aperture pixels")

    # Hann window (reduces spectral leakage at edges)
    hann = np.outer(np.hanning(N), np.hanning(N))
    Z_windowed = Z_grid * hann

    fft_result = fft2(Z_windowed)
    raw_power = np.abs(fft_result)**2

    # Build radial frequency grid in cycles/aperture
    dx = x_range[1] - x_range[0]
    dy = y_range[1] - y_range[0]
    freqs_x = fftfreq(N, dx)
    freqs_y = fftfreq(N, dy)
    fy, fx = np.meshgrid(freqs_y, freqs_x)
    D = 2 * aperture_radius
    freqs_radial = np.sqrt(fx**2 + fy**2) * D  # cycles/aperture

    # ── Band RMS via Parseval ──
    # Parseval: mean(|z|²) = (1/N⁴) Σ|FFT|²
    # But z is windowed, so we correct by the Hann window power:
    #   hann_power = mean(hann²) ≈ 0.375 for 2D Hann
    hann_power = float(np.mean(hann[apt_mask]**2))
    norm = 1.0 / (N * N * N * N * hann_power) if hann_power > 0 else 0

    def _band_rms(lo_cpa, hi_cpa):
        m = (freqs_radial >= lo_cpa) & (freqs_radial < hi_cpa)
        if not np.any(m):
            return 0.0
        return float(np.sqrt(np.sum(raw_power[m]) * norm))

    bands = {
        'gravity_rms': _band_rms(0.1, 2.0),
        'support_rms': _band_rms(2.0, 20.0),
        'hf_rms': _band_rms(20.0, N / 2),
        'total_rms': float(np.sqrt(np.sum(raw_power) * norm)),
    }

    # ── Radial PSD curve for plotting (arbitrary-unit, shape matters) ──
    psd_norm = dx * dy / (N * N)
    bin_edges = np.logspace(-1.5, np.log10(N / 2), 60)
    freqs_center, psd_values = [], []
    for i in range(len(bin_edges) - 1):
        mask = (freqs_radial >= bin_edges[i]) & (freqs_radial < bin_edges[i + 1])
        if np.any(mask):
            avg = float(np.mean(raw_power[mask]) * psd_norm)
            if avg > 0:
                psd_values.append(avg)
                freqs_center.append(float(np.sqrt(bin_edges[i] * bin_edges[i + 1])))

    return {
        'freqs': np.array(freqs_center),
        'psd': np.array(psd_values),
        'bands': bands,
    }


def make_psd_plot(psd_data_dict, title="Power Spectral Density"):
    """Create a log-log PSD plot for multiple angles.

    Args:
        psd_data_dict: {label: (freqs, psd)} for each angle/condition
    """
    colors = {'40° vs 20°': '#2563eb', '60° vs 20°': '#dc2626',
              '90° (Abs)': '#16a34a', '20° (Abs)': '#64748b'}
    fig = go.Figure()

    # Determine axis range from data
    all_freqs = []
    all_psd = []
    for label, (freqs, psd) in psd_data_dict.items():
        valid = psd > 0
        f_valid, p_valid = freqs[valid], psd[valid]
        all_freqs.extend(f_valid.tolist())
        all_psd.extend(p_valid.tolist())
        fig.add_trace(go.Scatter(
            x=f_valid, y=p_valid,
            mode='lines', name=label,
            line=dict(color=colors.get(label, '#6366f1'), width=2.5),
            hovertemplate="%{x:.1f} cyc/apt: %{y:.2e}<extra>" + label + "</extra>",
        ))

    # Compute axis limits from actual data
    if all_freqs:
        f_min = max(0.05, min(all_freqs) * 0.5)
        f_max = max(all_freqs) * 2.0
    else:
        f_min, f_max = 0.05, 200.0

    # Band annotations — clamped to data range
    band_max = min(f_max, 200.0)
    fig.add_vrect(x0=0.1, x1=2.0, fillcolor='rgba(59,130,246,0.07)', line_width=0, layer='below',
                  annotation_text="Gravity", annotation_position="top left",
                  annotation=dict(font=dict(size=11, color='#3b82f6')))
    fig.add_vrect(x0=2.0, x1=20.0, fillcolor='rgba(245,158,11,0.07)', line_width=0, layer='below',
                  annotation_text="Support", annotation_position="top left",
                  annotation=dict(font=dict(size=11, color='#f59e0b')))
    if band_max > 20:
        fig.add_vrect(x0=20.0, x1=band_max, fillcolor='rgba(239,68,68,0.05)', line_width=0, layer='below',
                      annotation_text="High Freq", annotation_position="top left",
                      annotation=dict(font=dict(size=11, color='#ef4444')))

    fig.update_layout(
        height=500, width=1100,
        margin=dict(t=30, b=60, l=80, r=30),
        **_PLOTLY_LIGHT_LAYOUT,
        xaxis=dict(type='log',
                   title=dict(text="Spatial Frequency [cycles/aperture]", font=dict(color='#1e293b', size=13)),
                   range=[np.log10(f_min), np.log10(f_max)],
                   gridcolor='#e2e8f0', tickfont=dict(color='#475569', size=11),
                   dtick=1),
        yaxis=dict(type='log',
                   title=dict(text="PSD [nm\u00b2\u00b7mm\u00b2]", font=dict(color='#1e293b', size=13)),
                   gridcolor='#e2e8f0', tickfont=dict(color='#475569', size=11)),
        legend=dict(x=0.70, y=0.98, bgcolor='rgba(255,255,255,0.92)',
                    bordercolor='#e2e8f0', borderwidth=1, font=dict(size=12, color='#1e293b')),
    )
    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="psd_plot")


def _psd_summary_html(band_dict, label=""):
    """Generate an HTML summary card for PSD band RMS values."""
    m = band_dict
    total = m['total_rms']
    grav_pct = 100 * m['gravity_rms'] / total if total > 0 else 0
    supp_pct = 100 * m['support_rms'] / total if total > 0 else 0
    hf_pct = 100 * m['hf_rms'] / total if total > 0 else 0

    return f"""
    <div class="sf-breakdown">
        <h4>{'PSD Band Decomposition — ' + label if label else 'PSD Band Decomposition (Tony Hull Methodology)'}</h4>
        <div class="sf-band-grid">
            <div class="sf-band-card sf-lsf">
                <div class="sf-band-label">Gravity Signature</div>
                <div class="sf-band-range">0.1\u20132 cyc/apt</div>
                <div class="sf-band-value">{m['gravity_rms']:.2f} <span class="sf-unit">nm</span></div>
                <div class="sf-band-range">{grav_pct:.0f}% of total</div>
            </div>
            <div class="sf-band-card sf-msf">
                <div class="sf-band-label">Support Print-through</div>
                <div class="sf-band-range">2\u201320 cyc/apt</div>
                <div class="sf-band-value">{m['support_rms']:.2f} <span class="sf-unit">nm</span></div>
                <div class="sf-band-range">{supp_pct:.0f}% of total</div>
            </div>
            <div class="sf-band-card sf-ratio">
                <div class="sf-band-label">High Frequency</div>
                <div class="sf-band-range">&gt;20 cyc/apt</div>
                <div class="sf-band-value">{m['hf_rms']:.2f} <span class="sf-unit">nm</span></div>
                <div class="sf-band-range">{hf_pct:.0f}% of total</div>
            </div>
        </div>
    </div>
    """


# ============================================================================
# HTML Report Generation
# ============================================================================

def _metric_color(value, target):
    """Return a CSS color class based on value vs target."""
    if value <= target:
        return 'color: #16a34a; font-weight: 700;'  # green
    ratio = value / target
    if ratio < 1.5:
        return 'color: #d97706; font-weight: 700;'  # amber
    return 'color: #dc2626; font-weight: 700;'  # red


def make_surface_plot(X, Y, W_res, mask, inner_radius=None, title="", amp=0.5, downsample=PLOT_DOWNSAMPLE, include_plotlyjs=False):
    """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", font=dict(color='#1e293b')), thickness=12, len=0.5, tickformat=".1f",
                              tickfont=dict(color='#1e293b')),
                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='#64748b', 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=dict(text="X (mm)", font=dict(color='#1e293b')), showgrid=True, gridcolor='#e2e8f0',
                       showbackground=True, backgroundcolor='#f1f5f9',
                       tickfont=dict(color='#475569')),
            yaxis=dict(title=dict(text="Y (mm)", font=dict(color='#1e293b')), showgrid=True, gridcolor='#e2e8f0',
                       showbackground=True, backgroundcolor='#f1f5f9',
                       tickfont=dict(color='#475569')),
            zaxis=dict(title=dict(text="Residual (nm)", font=dict(color='#1e293b')), range=[-max_amp*3.0, max_amp*3.0],
                       showgrid=True, gridcolor='#e2e8f0',
                       showbackground=True, backgroundcolor='#eff6ff',
                       tickfont=dict(color='#475569')),
            aspectmode='manual',
            aspectratio=dict(x=1, y=1, z=0.4),
        ),
        margin=dict(t=30, b=10, l=10, r=10),
        **_PLOTLY_LIGHT_LAYOUT,
        height=650,
        width=1200,
    )
    div_id = f"surface_{title.replace(' ','_')}"
    config = get_plotly_config(f"WFE_Surface_{title.replace(' ','_')}")
    return fig.to_html(include_plotlyjs=include_plotlyjs, full_html=False, div_id=div_id, config=config)


def make_bar_chart(coeffs, title="Zernike Coefficients", max_modes=50):
    """Create horizontal bar chart of Zernike coefficient magnitudes."""
    n = min(len(coeffs), max_modes)
    labels = [str(zernike_label(j))[:40] for j in range(1, n+1)]
    vals = np.abs(coeffs[:n])
    text_vals = [f"{v:.3f}" for v in vals]

    fig = go.Figure(go.Bar(
        x=vals, y=labels, orientation='h',
        marker_color=_BLUE_PALETTE[0],
        text=text_vals,
        textposition='outside',
        textfont=dict(size=9, color='#475569'),
        hovertemplate="%{y}<br>|c| = %{x:.3f} nm<extra></extra>",
    ))
    fig.update_layout(
        height=max(400, n*22),
        margin=dict(t=30, b=10, l=220, r=60),
        **_PLOTLY_LIGHT_LAYOUT,
        xaxis=dict(title=dict(text="|Coefficient| (nm)", font=dict(color='#1e293b')),
                   gridcolor='#e2e8f0', zeroline=True,
                   zerolinecolor='#cbd5e1',
                   tickfont=dict(color='#475569')),
        yaxis=dict(autorange='reversed', tickfont=dict(size=10, color='#1e293b')),
    )
    div_id = f"bar_{title.replace(' ','_')}"
    config = get_plotly_config(f"Zernike_Coefficients_{title.replace(' ','_')}")
    return fig.to_html(include_plotlyjs=False, full_html=False, div_id=div_id, config=config)


def make_trajectory_plot(angles, coefficients_relative, mode_groups, sensitivity, title=""):
    """Create trajectory visualization: Zernike modes vs elevation angle."""
    fig = go.Figure()

    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

        color = _MODE_COLORS[color_idx % len(_MODE_COLORS)]
        dash = _MODE_DASHES[color_idx % len(_MODE_DASHES)]

        # RSS of modes in this group at each angle
        rss = np.sqrt(np.sum(coefficients_relative[:, indices]**2, axis=1))
        display_name = MODE_NAMES.get(group_name, group_name)

        # Group RSS trace (thick)
        fig.add_trace(go.Scatter(
            x=angles, y=rss,
            mode='lines+markers',
            name=f"{display_name} (RSS)",
            line=dict(color=color, width=3, dash=dash),
            marker=dict(size=10, symbol='circle'),
            hovertemplate=f"{display_name} RSS<br>%{{x:.1f}}°: %{{y:.3f}} nm<extra></extra>",
            legendgroup=group_name,
        ))

        # Individual mode traces (thin, same color, lighter)
        for idx_i, noll_idx in enumerate(noll_indices):
            col_idx = noll_idx - 5
            if col_idx < 0 or col_idx >= coefficients_relative.shape[1]:
                continue
            mode_vals = np.abs(coefficients_relative[:, col_idx])
            mode_label = f"J{noll_idx}"
            fig.add_trace(go.Scatter(
                x=angles, y=mode_vals,
                mode='lines+markers',
                name=f"  {mode_label}",
                line=dict(color=color, width=1, dash='dot'),
                marker=dict(size=5, symbol='diamond'),
                opacity=0.5,
                hovertemplate=f"{mode_label}<br>%{{x:.1f}}°: %{{y:.3f}} nm<extra></extra>",
                legendgroup=group_name,
                showlegend=True,
            ))

        color_idx += 1

    # Check if log scale would help (values span >2 orders of magnitude)
    all_y = []
    for trace in fig.data:
        all_y.extend([v for v in trace.y if v > 0])
    use_log = False
    if all_y:
        y_min, y_max = min(all_y), max(all_y)
        if y_max > 0 and y_min > 0 and (y_max / y_min) > 100:
            use_log = True

    yaxis_cfg = dict(
        title=dict(text="RMS (nm)", font=dict(color='#1e293b')),
        gridcolor='#e2e8f0',
        zeroline=True, zerolinecolor='#cbd5e1',
        tickfont=dict(color='#475569'),
    )
    if use_log:
        yaxis_cfg['type'] = 'log'
        yaxis_cfg['title'] = dict(text="RMS (nm) — log scale", font=dict(color='#1e293b'))

    fig.update_layout(
        height=500,
        width=1100,
        margin=dict(t=30, b=50, l=70, r=20),
        **_PLOTLY_LIGHT_LAYOUT,
        xaxis=dict(title=dict(text="Elevation Angle (°)", font=dict(color='#1e293b')),
                   gridcolor='#e2e8f0',
                   tickvals=list(angles), dtick=10,
                   tickfont=dict(color='#475569')),
        yaxis=yaxis_cfg,
        legend=dict(x=0.01, y=0.99, bgcolor='rgba(255,255,255,0.9)',
                    bordercolor='#e2e8f0', borderwidth=1,
                    font=dict(size=10, color='#1e293b')),
    )
    config = get_plotly_config("Zernike_Trajectory")
    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="trajectory_plot", config=config)


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 = [str(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 \u03b8)', marker_color='#2563eb',
        text=[f"{v:.3f}" for v in axial], textposition='outside',
        textfont=dict(size=9, color='#475569'),
        hovertemplate="%{y}<br>Axial: %{x:.3f} nm/unit<extra></extra>"
    ))
    fig.add_trace(go.Bar(
        y=labels, x=lateral, orientation='h',
        name='Lateral (cos \u03b8)', marker_color='#7c3aed',
        text=[f"{v:.3f}" for v in lateral], textposition='outside',
        textfont=dict(size=9, color='#475569'),
        hovertemplate="%{y}<br>Lateral: %{x:.3f} nm/unit<extra></extra>"
    ))
    fig.update_layout(
        barmode='group',
        height=max(300, len(modes)*45),
        margin=dict(t=30, b=40, l=180, r=60),
        **_PLOTLY_LIGHT_LAYOUT,
        xaxis=dict(title=dict(text="Sensitivity (nm per load fraction)", font=dict(color='#1e293b')),
                   gridcolor='#e2e8f0',
                   zeroline=True, zerolinecolor='#cbd5e1',
                   tickfont=dict(color='#475569')),
        yaxis=dict(autorange='reversed', tickfont=dict(size=11, color='#1e293b')),
        legend=dict(x=0.55, y=0.99, bgcolor='rgba(255,255,255,0.9)',
                    bordercolor='#e2e8f0', borderwidth=1,
                    font=dict(color='#1e293b')),
    )
    config = get_plotly_config("Sensitivity_Axial_vs_Lateral")
    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="sensitivity_bar", config=config)


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}\u00b0 vs {ref_angle}\u00b0" for a in angles]

    fig = go.Figure(go.Bar(
        x=labels, y=rms_vals,
        marker_color=_BLUE_PALETTE[0],
        text=[f"{v:.2f} nm" for v in rms_vals],
        textposition='outside',
        textfont=dict(size=11, color='#1e293b'),
        hovertemplate="%{x}: %{y:.2f} nm<extra></extra>"
    ))
    fig.update_layout(
        height=350,
        margin=dict(t=30, b=40, l=60, r=20),
        **_PLOTLY_LIGHT_LAYOUT,
        yaxis=dict(title=dict(text="Filtered RMS WFE (nm)", font=dict(color='#1e293b')),
                   gridcolor='#e2e8f0',
                   zeroline=True, zerolinecolor='#cbd5e1',
                   tickfont=dict(color='#475569')),
        xaxis=dict(tickfont=dict(color='#1e293b', size=12)),
    )
    config = get_plotly_config("Per_Angle_RMS_WFE")
    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="per_angle_rms", config=config)


# ============================================================================
# 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
    # First plot embeds the full Plotly library (~3.5MB) for offline viewing
    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", include_plotlyjs=True
    )
    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")

    # PSD analysis (Tony Hull methodology)
    print("\n[PSD] Computing power spectral density...")
    R_outer_est = float(np.max(np.hypot(
        angle_results[40]['X'] - np.mean(angle_results[40]['X']),
        angle_results[40]['Y'] - np.mean(angle_results[40]['Y']))))

    psd_plot_data = {}
    psd_bands = {}
    # Use Zernike-filtered residual surface (J1-J4 removed) for PSD —
    # avoids piston/tilt/defocus dominating and circular-to-square boundary artifacts
    for ang_label, ang_key, use_rel in [
        ('40° vs 20°', 40, True),
        ('60° vs 20°', 60, True),
        ('90° (Abs)', 90, False),
    ]:
        r = angle_results[ang_key]
        rms_data = r['rms_rel'] if use_rel else r['rms_abs']
        Xp = r['X_rel'] if use_rel else r['X']
        Yp = r['Y_rel'] if use_rel else r['Y']
        Zp = rms_data['W_res_filt']  # J1-J4 filtered residual
        mask = rms_data['mask']
        # Only pass masked (valid aperture) points
        Xm, Ym, Zm = Xp[mask], Yp[mask], Zp[mask]
        try:
            result = compute_surface_psd(Xm, Ym, Zm, R_outer_est)
            psd_plot_data[ang_label] = (result['freqs'], result['psd'])
            psd_bands[ang_label] = result['bands']
            b = result['bands']
            print(f"  {ang_label}: gravity={b['gravity_rms']:.2f} nm, "
                  f"support={b['support_rms']:.2f} nm, "
                  f"total={b['total_rms']:.2f} nm")
        except Exception as e:
            print(f"  [WARN] PSD for {ang_label} failed: {e}")

    psd_plot_html = make_psd_plot(psd_plot_data) if psd_plot_data else ""
    psd_summary_40 = _psd_summary_html(psd_bands['40° vs 20°'], '40° vs 20° (Relative)') if '40° vs 20°' in psd_bands else ""
    psd_summary_60 = _psd_summary_html(psd_bands['60° vs 20°'], '60° vs 20° (Relative)') if '60° vs 20°' in psd_bands else ""
    psd_summary_90 = _psd_summary_html(psd_bands['90° (Abs)'], '90° Manufacturing (Absolute)') if '90° (Abs)' in psd_bands else ""

    # 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 = "\u00b0" 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>7. 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}\u00b0</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\u00b2</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: {' \u2192 '.join(traj_result['mode_ranking'][:5])}</p>
        """

    # Manufacturing details
    mfg_html = f"""
    <table class="data-table">
        <thead><tr><th>Metric</th><th>Absolute 90\u00b0</th><th>Correction (90\u00b0\u221220\u00b0)</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\u2212J3 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>
    """

    # Executive summary metric styling
    style_40 = _metric_color(wfe_40_20, targets['wfe_40_20'])
    style_60 = _metric_color(wfe_60_20, targets['wfe_60_20'])
    style_mfg = _metric_color(mfg_90, targets['mfg_90'])

    # Section numbering: adjust if trajectory present
    sec_surface = 3
    sec_traj = 4
    sec_psd = 5 if traj_result else 4
    sec_mfg = sec_psd + 1
    sec_params = sec_mfg + 1
    sec_zernike = sec_params + 1 if (study_params and study_params.get('parameters')) else sec_mfg + 1
    sec_method = sec_zernike + 1

    # 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>
<!-- Plotly.js is embedded inline in the first surface plot for offline viewing -->
<style>
    :root {{
        --bg-primary: #ffffff;
        --bg-secondary: #f8fafc;
        --bg-card: #ffffff;
        --text-primary: #1e293b;
        --text-secondary: #64748b;
        --accent: #2563eb;
        --accent-light: #dbeafe;
        --success: #16a34a;
        --warning: #d97706;
        --danger: #dc2626;
        --border: #e2e8f0;
        --border-strong: #cbd5e1;
    }}
    * {{ margin: 0; padding: 0; box-sizing: border-box; }}
    body {{
        font-family: 'Inter', -apple-system, BlinkMacSystemFont, 'Segoe UI', system-ui, sans-serif;
        background: var(--bg-secondary);
        color: var(--text-primary);
        line-height: 1.6;
    }}
    .container {{ max-width: 1400px; margin: 0 auto; padding: 2rem; }}

    /* Header */
    .header {{
        background: var(--bg-card);
        border: 1px solid var(--border);
        border-radius: 8px;
        padding: 2rem 3rem;
        margin-bottom: 2rem;
        display: flex;
        justify-content: space-between;
        align-items: center;
    }}
    .header h1 {{ font-size: 1.6rem; font-weight: 700; color: var(--text-primary); }}
    .header .subtitle {{ color: var(--text-secondary); font-size: 0.9rem; margin-top: 0.3rem; }}
    .header .branding {{
        text-align: right;
        font-size: 0.85rem;
        color: var(--text-secondary);
    }}
    .header .branding .by-line {{ color: #94a3b8; font-size: 0.8rem; margin-top: 0.2rem; }}
    .header .branding .tagline {{ color: var(--text-secondary); font-size: 0.8rem; margin-top: 0.15rem; }}

    /* Sections */
    .section {{
        background: var(--bg-card);
        border: 1px solid var(--border);
        border-radius: 8px;
        padding: 1.5rem 2rem;
        margin-bottom: 1.5rem;
    }}
    .section h2 {{
        font-size: 1.2rem;
        font-weight: 600;
        color: var(--text-primary);
        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(280px, 1fr));
        gap: 1rem;
        margin: 1rem 0;
    }}
    .exec-card {{
        background: var(--bg-secondary);
        border: 1px solid var(--border);
        border-radius: 6px;
        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 .unit {{ font-size: 0.9rem; font-weight: 400; color: var(--text-secondary); }}
    .exec-footnote {{
        margin-top: 1rem;
        padding: 0.75rem 1rem;
        background: var(--bg-secondary);
        border-left: 3px solid var(--accent);
        border-radius: 0 4px 4px 0;
        font-size: 0.82rem;
        color: var(--text-secondary);
    }}

    /* 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-secondary);
        font-weight: 600;
        font-size: 0.82rem;
        text-transform: uppercase;
        letter-spacing: 0.05em;
        color: var(--text-secondary);
    }}
    .data-table tr:hover {{ background: #f1f5f9; }}
    .data-table tr.highlight td {{
        background: var(--accent-light);
        font-weight: 600;
    }}

    /* Metrics Grid */
    .metrics-grid {{
        display: grid;
        grid-template-columns: repeat(auto-fit, minmax(170px, 1fr));
        gap: 0.8rem;
        margin: 1rem 0;
    }}
    .metric-card {{
        background: var(--bg-secondary);
        border: 1px solid var(--border);
        border-radius: 6px;
        padding: 1rem;
        text-align: center;
    }}
    .metric-label {{ font-size: 0.8rem; color: var(--text-secondary); margin-bottom: 0.3rem; }}
    .metric-value {{ font-size: 1.2rem; font-weight: 700; color: var(--accent); }}

    /* Plots */
    .plot-grid {{
        display: grid;
        grid-template-columns: repeat(auto-fit, minmax(1100px, 1fr));
        gap: 1rem;
    }}
    .plot-container {{
        background: var(--bg-card);
        border: 1px solid var(--border);
        border-radius: 6px;
        padding: 1rem;
        min-height: 400px;
    }}
    .plot-container h3 {{
        font-size: 0.95rem;
        font-weight: 600;
        margin-bottom: 0.5rem;
        color: var(--text-secondary);
    }}

    /* Spatial Frequency Breakdown */
    .sf-breakdown {{
        margin: 1rem 0 1.5rem 0;
        padding: 1.2rem;
        background: var(--bg-secondary);
        border: 1px solid var(--border);
        border-radius: 8px;
    }}
    .sf-breakdown h4 {{
        font-size: 0.9rem;
        font-weight: 600;
        color: var(--text-secondary);
        text-transform: uppercase;
        letter-spacing: 0.05em;
        margin-bottom: 0.8rem;
    }}
    .sf-band-grid {{
        display: grid;
        grid-template-columns: 1fr 1fr auto;
        gap: 0.8rem;
        margin-bottom: 1rem;
    }}
    .sf-band-card {{
        background: var(--bg-card);
        border: 1px solid var(--border);
        border-radius: 6px;
        padding: 0.8rem 1rem;
        text-align: center;
    }}
    .sf-band-label {{
        font-size: 0.82rem;
        font-weight: 600;
        color: var(--text-primary);
    }}
    .sf-band-range {{
        font-size: 0.75rem;
        color: var(--text-secondary);
        margin: 0.15rem 0 0.4rem 0;
    }}
    .sf-band-value {{
        font-size: 1.3rem;
        font-weight: 700;
        color: var(--accent);
    }}
    .sf-unit {{
        font-size: 0.8rem;
        font-weight: 400;
        color: var(--text-secondary);
    }}
    .sf-lsf {{ border-top: 3px solid #2563eb; }}
    .sf-msf {{ border-top: 3px solid #7c3aed; }}
    .sf-ratio {{ border-top: 3px solid #16a34a; min-width: 100px; }}
    .sf-order-grid {{
        background: var(--bg-card);
        border: 1px solid var(--border);
        border-radius: 6px;
        padding: 0.8rem 1rem;
    }}
    .sf-order-title {{
        font-size: 0.8rem;
        font-weight: 600;
        color: var(--text-secondary);
        margin-bottom: 0.5rem;
    }}
    .sf-order-items {{
        display: flex;
        flex-wrap: wrap;
        gap: 0.3rem 1.2rem;
    }}
    .sf-order-item {{
        display: flex;
        gap: 0.3rem;
        font-size: 0.85rem;
    }}
    .sf-order-n {{
        color: var(--text-secondary);
        font-weight: 500;
    }}
    .sf-order-val {{
        color: var(--text-primary);
        font-weight: 600;
    }}

    /* Collapsible */
    details {{ margin: 0.5rem 0; }}
    summary {{
        cursor: pointer;
        font-weight: 600;
        padding: 0.6rem 0.8rem;
        background: var(--bg-secondary);
        border-radius: 6px;
        border: 1px solid var(--border);
        color: var(--text-primary);
        font-size: 0.95rem;
    }}
    summary:hover {{ background: var(--accent-light); }}
    details > div {{ padding: 1rem; }}

    .note {{ color: var(--text-secondary); font-size: 0.88rem; margin: 0.5rem 0; }}

    .two-col {{ display: grid; grid-template-columns: 1fr 1fr; gap: 1.5rem; }}
    @media (max-width: 900px) {{ .two-col {{ grid-template-columns: 1fr; }} }}

    /* Print styles */
    @media print {{
        body {{ background: white; color: black; font-size: 10pt; }}
        .container {{ max-width: 100%; padding: 0.5cm; }}
        .header {{ border: 1px solid #999; }}
        .section {{ border: 1px solid #ccc; page-break-inside: avoid; margin-bottom: 0.5cm; }}
        .plot-container {{ page-break-inside: avoid; }}
        details {{ display: block; }}
        details > summary {{ display: none; }}
        details > div {{ padding: 0; }}
        .exec-grid {{ grid-template-columns: repeat(2, 1fr); }}
        .no-print {{ display: none; }}
    }}

    /* Footer */
    .footer {{
        text-align: center;
        padding: 2rem;
        color: var(--text-secondary);
        font-size: 0.8rem;
        border-top: 1px solid var(--border);
        margin-top: 1rem;
    }}
</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">
        <svg width="120" height="32" viewBox="0 0 120 32" xmlns="http://www.w3.org/2000/svg">
            <text x="0" y="24" font-family="Inter, system-ui, sans-serif" font-size="22" font-weight="700" fill="#2563eb">ATOMIZER</text>
        </svg>
        <div class="by-line">by Atomaste</div>
        <div class="tagline">FEA Optimization Platform</div>
    </div>
</div>

<!-- 1. Executive Summary -->
<div class="section">
    <h2>1. Executive Summary</h2>
    <div class="exec-grid">
        <div class="exec-card">
            <div class="label">WFE 40\u00b0 vs 20\u00b0 (Tracking)</div>
            <div class="value" style="{style_40}">{wfe_40_20:.2f} <span class="unit">nm</span></div>
        </div>
        <div class="exec-card">
            <div class="label">WFE 60\u00b0 vs 20\u00b0 (Tracking)</div>
            <div class="value" style="{style_60}">{wfe_60_20:.2f} <span class="unit">nm</span></div>
        </div>
        <div class="exec-card">
            <div class="label">MFG 90\u00b0 (J1\u2212J3 Filtered)</div>
            <div class="value" style="{style_mfg}">{mfg_90:.2f} <span class="unit">nm</span></div>
        </div>
        <div class="exec-card">
            <div class="label">Weighted Sum (6\u00b7W40 + 5\u00b7W60 + 3\u00b7MFG)</div>
            <div class="value" style="color: var(--accent); font-weight: 700;">{ws:.1f}</div>
        </div>
    </div>
    <div class="exec-footnote">
        Design targets &mdash;
        WFE 40\u00b0\u221220\u00b0 \u2264 {targets['wfe_40_20']:.1f} nm &nbsp;|&nbsp;
        WFE 60\u00b0\u221220\u00b0 \u2264 {targets['wfe_60_20']:.1f} nm &nbsp;|&nbsp;
        MFG 90\u00b0 \u2264 {targets['mfg_90']:.1f} nm &nbsp;|&nbsp;
        Weighted sum: lower is better.
        {'<br>Annular aperture: inner radius = ' + f'{inner_radius:.1f} mm (\u00f8{2*inner_radius:.1f} mm central hole)' if inner_radius else ''}
    </div>
</div>

<!-- 2. Per-Angle Summary -->
<div class="section">
    <h2>2. 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\u2212J4 removed. Relative = vs 20\u00b0 reference. Aberrations are absolute.</p>
</div>

<!-- 3. Surface Plots -->
<div class="section">
    <h2>{sec_surface}. Wavefront Error Surface Maps</h2>
    <p class="note">3D residual surfaces after removing piston, tip, tilt, and defocus (J1\u2212J4). Interactive \u2014 drag to rotate.</p>
    <div class="plot-grid">
        <div class="plot-container">
            <h3>40\u00b0 vs 20\u00b0 (Relative)</h3>
            {surf_40}
        </div>
        <div class="plot-container">
            <h3>60\u00b0 vs 20\u00b0 (Relative)</h3>
            {surf_60}
        </div>
    </div>
    <div class="plot-container" style="margin-top:1rem">
        <h3>90\u00b0 Manufacturing (Absolute)</h3>
        {surf_90}
    </div>
</div>

<!-- Trajectory Analysis -->
{'<div class="section"><h2>' + str(sec_traj) + '. Zernike Trajectory Analysis</h2>' +
 '<p class="note">Mode-specific integrated RMS across the operating elevation range. ' +
 'The linear model c<sub>j</sub>(\u03b8) = a<sub>j</sub>\u00b7\u0394sin\u03b8 + b<sub>j</sub>\u00b7\u0394cos\u03b8 decomposes gravity into axial and lateral components.</p>' +
 traj_metrics_html +
 '<div class="plot-container" style="margin-top:1rem"><h3>Mode RMS vs Elevation Angle</h3>' + traj_plot_html + '</div>' +
 '<div class="plot-container" style="margin-top:1rem"><h3>Axial vs Lateral Sensitivity</h3>' + sens_plot_html + '</div>' +
 '</div>' if traj_result else ''}

<!-- PSD Analysis -->
{'<div class="section"><h2>' + str(sec_psd) + '. Power Spectral Density Analysis</h2>' +
 '<p class="note">Surface PSD computed via 2D FFT with Hann windowing and radial averaging ' +
 '(Tony Hull / JWST methodology). Frequency bands: Gravity signature (0.1\u20132 cyc/apt), ' +
 'Support print-through (2\u201320 cyc/apt), High frequency (&gt;20 cyc/apt).</p>' +
 '<div class="plot-container" style="margin:1rem 0"><h3>PSD \u2014 Log-Log</h3>' + psd_plot_html + '</div>' +
 psd_summary_40 + psd_summary_60 + psd_summary_90 +
 '</div>' if psd_plot_html else ''}

<!-- Manufacturing Analysis -->
<div class="section">
    <h2>{sec_mfg}. Manufacturing Analysis (90\u00b0 Orientation)</h2>
    <p class="note">
        The mirror is manufactured (polished) at 90\u00b0 orientation. The "Correction" column shows the
        aberrations that must be polished out to achieve the 20\u00b0 operational figure.
    </p>
    {mfg_html}
</div>

<!-- Design Parameters -->
{params_html}

<!-- Zernike Coefficient Details -->
<div class="section">
    <h2>{sec_zernike}. Zernike Coefficient Details</h2>
    <details>
        <summary>40\u00b0 vs 20\u00b0 \u2014 Relative Coefficients</summary>
        <div>{bar_40}</div>
    </details>
    <details>
        <summary>60\u00b0 vs 20\u00b0 \u2014 Relative Coefficients</summary>
        <div>{bar_60}</div>
    </details>
    <details>
        <summary>90\u00b0 \u2014 Absolute Coefficients</summary>
        <div>{bar_90}</div>
    </details>
</div>

<!-- Methodology -->
<div class="section">
    <h2>{sec_method}. 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\u2212J4 (Piston, Tip, Tilt, Defocus)</td></tr>
            <tr><td><b>WFE Calculation</b></td><td>WFE = 2 \u00d7 Surface Error (reflective)</td></tr>
            <tr><td><b>Displacement Unit</b></td><td>{DISP_UNIT} \u2192 nm ({NM_SCALE:.0e}\u00d7)</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\u00b0 (polishing/measurement orientation)</td></tr>
            <tr><td><b>MFG Objective</b></td><td>90\u00b0\u221220\u00b0 relative, J1\u2212J3 filtered (optician workload)</td></tr>
            <tr><td><b>Weighted Sum</b></td><td>6\u00d7WFE(40\u221220) + 5\u00d7WFE(60\u221220) + 3\u00d7MFG(90)</td></tr>
            {'<tr><td><b>Trajectory R\u00b2</b></td><td>' + f'{traj_result["linear_fit_r2"]:.6f}' + '</td></tr>' if traj_result else ''}
        </tbody>
    </table>
</div>

<!-- Footer -->
<div class="footer">
    Generated by <b>Atomizer</b> Optical Report Generator &nbsp;|&nbsp; {timestamp}<br>
    \u00a9 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()