"""
Zernike Wavefront Error (WFE) Insight

Provides 3D surface visualization of mirror wavefront errors with
Zernike polynomial decomposition. Generates three views:
- 40 deg vs 20 deg (operational tilt comparison)
- 60 deg vs 20 deg (operational tilt comparison)
- 90 deg Manufacturing (absolute with optician workload metrics)

Applicable to: Mirror optimization studies with multi-subcase gravity loads.
"""

from pathlib import Path
from datetime import datetime
from typing import Dict, Any, List, Optional, Tuple
import numpy as np
from math import factorial
from numpy.linalg import LinAlgError

from .base import StudyInsight, InsightConfig, InsightResult, register_insight

# Lazy imports to avoid startup overhead
_plotly_loaded = False
_go = None
_make_subplots = None
_Triangulation = None
_OP2 = None
_BDF = None


def _load_dependencies():
    """Lazy load heavy dependencies."""
    global _plotly_loaded, _go, _make_subplots, _Triangulation, _OP2, _BDF
    if not _plotly_loaded:
        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
        _go = go
        _make_subplots = make_subplots
        _Triangulation = Triangulation
        _OP2 = OP2
        _BDF = BDF
        _plotly_loaded = True


# ============================================================================
# Zernike Mathematics
# ============================================================================
def noll_indices(j: int) -> Tuple[int, int]:
    """Convert Noll index to (n, m) radial/azimuthal orders."""
    if j < 1:
        raise ValueError("Noll index j must be >= 1")
    count = 0
    n = 0
    while True:
        if n == 0:
            ms = [0]
        elif n % 2 == 0:
            ms = [0] + [m for k in range(1, n//2 + 1) for m in (-2*k, 2*k)]
        else:
            ms = [m for k in range(0, (n+1)//2) for m in (-(2*k+1), (2*k+1))]
        for m in ms:
            count += 1
            if count == j:
                return n, m
        n += 1


def zernike_noll(j: int, r: np.ndarray, th: np.ndarray) -> np.ndarray:
    """Evaluate Zernike polynomial j at (r, theta)."""
    n, m = noll_indices(j)
    R = np.zeros_like(r)
    for s in range((n - abs(m)) // 2 + 1):
        c = ((-1)**s * factorial(n - s) /
             (factorial(s) *
              factorial((n + abs(m)) // 2 - s) *
              factorial((n - abs(m)) // 2 - s)))
        R += c * r**(n - 2*s)
    if m == 0:
        return R
    return R * (np.cos(m * th) if m > 0 else np.sin(-m * th))


def zernike_common_name(n: int, m: int) -> str:
    """Get common name for Zernike mode."""
    names = {
        (0, 0): "Piston", (1, -1): "Tilt X", (1, 1): "Tilt Y",
        (2, 0): "Defocus", (2, -2): "Astig 45°", (2, 2): "Astig 0°",
        (3, -1): "Coma X", (3, 1): "Coma Y", (3, -3): "Trefoil X", (3, 3): "Trefoil Y",
        (4, 0): "Primary Spherical", (4, -2): "Sec Astig X", (4, 2): "Sec Astig Y",
        (4, -4): "Quadrafoil X", (4, 4): "Quadrafoil Y",
        (5, -1): "Sec Coma X", (5, 1): "Sec Coma Y",
        (5, -3): "Sec Trefoil X", (5, 3): "Sec Trefoil Y",
        (5, -5): "Pentafoil X", (5, 5): "Pentafoil Y",
        (6, 0): "Sec Spherical",
    }
    return names.get((n, m), f"Z(n={n}, m={m})")


def zernike_label(j: int) -> str:
    """Get label for Zernike coefficient J{j}."""
    n, m = noll_indices(j)
    return f"J{j:02d} - {zernike_common_name(n, m)} (n={n}, m={m})"


def compute_zernike_coeffs(
    X: np.ndarray,
    Y: np.ndarray,
    vals: np.ndarray,
    n_modes: int,
    chunk_size: int = 100000
) -> Tuple[np.ndarray, float]:
    """Fit Zernike coefficients to WFE data."""
    Xc, Yc = X - np.mean(X), Y - np.mean(Y)
    R = float(np.max(np.hypot(Xc, Yc)))
    r = np.hypot(Xc / R, Yc / R).astype(np.float32)
    th = np.arctan2(Yc, Xc).astype(np.float32)
    mask = (r <= 1.0) & ~np.isnan(vals)
    if not np.any(mask):
        raise RuntimeError("No valid points inside unit disk.")
    idx = np.nonzero(mask)[0]
    m = int(n_modes)
    G = np.zeros((m, m), dtype=np.float64)
    h = np.zeros((m,), dtype=np.float64)
    v = vals.astype(np.float64)

    for start in range(0, len(idx), chunk_size):
        sl = idx[start:start + chunk_size]
        r_b, th_b, v_b = r[sl], th[sl], v[sl]
        Zb = np.column_stack([zernike_noll(j, r_b, th_b).astype(np.float32)
                              for j in range(1, m + 1)])
        G += (Zb.T @ Zb).astype(np.float64)
        h += (Zb.T @ v_b).astype(np.float64)
    try:
        coeffs = np.linalg.solve(G, h)
    except LinAlgError:
        coeffs = np.linalg.lstsq(G, h, rcond=None)[0]
    return coeffs, R


# ============================================================================
# Configuration Defaults
# ============================================================================
DEFAULT_CONFIG = {
    'n_modes': 50,
    'amp': 0.5,                    # Visual deformation scale
    'pancake': 3.0,                # Z-axis range multiplier
    'plot_downsample': 10000,
    'filter_low_orders': 4,        # Piston, tip, tilt, defocus
    'colorscale': 'Turbo',
    'disp_unit': 'mm',
    'show_bar_chart': True,
}


@register_insight
class ZernikeWFEInsight(StudyInsight):
    """
    Zernike Wavefront Error visualization for mirror optimization.

    Generates interactive 3D surface plots showing:
    - Residual WFE after Zernike fit
    - Coefficient bar charts
    - RMS metrics tables
    - Manufacturing orientation analysis
    """

    insight_type = "zernike_wfe"
    name = "Zernike WFE Analysis"
    description = "3D wavefront error surface with Zernike decomposition"
    applicable_to = ["mirror", "optics", "wfe"]
    required_files = ["*.op2"]

    def __init__(self, study_path: Path):
        super().__init__(study_path)
        self.op2_path: Optional[Path] = None
        self.geo_path: Optional[Path] = None
        self._node_geo: Optional[Dict] = None
        self._displacements: Optional[Dict] = None

    def can_generate(self) -> bool:
        """Check if OP2 and geometry files exist."""
        # Look for OP2 in results or iterations
        search_paths = [
            self.results_path,
            self.study_path / "2_iterations",
            self.setup_path / "model",
        ]

        for search_path in search_paths:
            if not search_path.exists():
                continue
            op2_files = list(search_path.glob("**/*solution*.op2"))
            if not op2_files:
                op2_files = list(search_path.glob("**/*.op2"))
            if op2_files:
                self.op2_path = max(op2_files, key=lambda p: p.stat().st_mtime)
                break

        if self.op2_path is None:
            return False

        # Find geometry
        try:
            self.geo_path = self._find_geometry_file(self.op2_path)
            return True
        except FileNotFoundError:
            return False

    def _find_geometry_file(self, op2_path: Path) -> Path:
        """Find BDF/DAT geometry file for OP2."""
        folder = op2_path.parent
        base = op2_path.stem

        for ext in ['.dat', '.bdf']:
            cand = folder / (base + ext)
            if cand.exists():
                return cand

        for f in folder.iterdir():
            if f.suffix.lower() in ['.dat', '.bdf']:
                return f

        raise FileNotFoundError(f"No geometry file found for {op2_path}")

    def _load_data(self):
        """Load geometry and displacement data."""
        if self._node_geo is not None:
            return  # Already loaded

        _load_dependencies()

        # Read geometry
        bdf = _BDF()
        bdf.read_bdf(str(self.geo_path))
        self._node_geo = {int(nid): node.get_position()
                         for nid, node in bdf.nodes.items()}

        # Read displacements
        op2 = _OP2()
        op2.read_op2(str(self.op2_path))

        if not op2.displacements:
            raise RuntimeError("No displacement data in OP2")

        self._displacements = {}
        for key, darr in op2.displacements.items():
            data = darr.data
            dmat = data[0] if data.ndim == 3 else (data if data.ndim == 2 else None)
            if dmat is None:
                continue

            ngt = darr.node_gridtype.astype(int)
            node_ids = ngt if ngt.ndim == 1 else ngt[:, 0]
            isubcase = getattr(darr, 'isubcase', None)
            label = str(isubcase) if isubcase else str(key)

            self._displacements[label] = {
                'node_ids': node_ids.astype(int),
                'disp': dmat.copy()
            }

    def _build_wfe_arrays(
        self,
        label: str,
        disp_unit: str = 'mm'
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
        """Build X, Y, WFE arrays for a subcase."""
        nm_per_unit = 1e6 if disp_unit == 'mm' else 1e9

        data = self._displacements[label]
        node_ids = data['node_ids']
        dmat = data['disp']

        X, Y, WFE = [], [], []
        valid_nids = []
        for nid, vec in zip(node_ids, dmat):
            geo = self._node_geo.get(int(nid))
            if geo is None:
                continue
            X.append(geo[0])
            Y.append(geo[1])
            wfe = vec[2] * 2.0 * nm_per_unit  # Z-disp to WFE
            WFE.append(wfe)
            valid_nids.append(nid)

        return (np.array(X), np.array(Y), np.array(WFE), np.array(valid_nids))

    def _compute_relative_wfe(
        self,
        X1, Y1, WFE1, nids1,
        X2, Y2, WFE2, nids2
    ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Compute WFE1 - WFE2 for common nodes."""
        ref_map = {int(nid): (x, y, w) for nid, x, y, w in zip(nids2, X2, Y2, WFE2)}

        X_rel, Y_rel, WFE_rel = [], [], []
        for nid, x, y, w in zip(nids1, X1, Y1, WFE1):
            nid = int(nid)
            if nid in ref_map:
                _, _, w_ref = ref_map[nid]
                X_rel.append(x)
                Y_rel.append(y)
                WFE_rel.append(w - w_ref)

        return np.array(X_rel), np.array(Y_rel), np.array(WFE_rel)

    def _compute_metrics(
        self,
        X: np.ndarray,
        Y: np.ndarray,
        W_nm: np.ndarray,
        n_modes: int,
        filter_orders: int
    ) -> Dict[str, Any]:
        """Compute RMS metrics and Zernike coefficients."""
        coeffs, R = compute_zernike_coeffs(X, Y, W_nm, n_modes)

        Xc = X - np.mean(X)
        Yc = Y - np.mean(Y)
        r = np.hypot(Xc / R, Yc / R)
        th = np.arctan2(Yc, Xc)

        Z = np.column_stack([zernike_noll(j, r, th) for j in range(1, n_modes + 1)])
        W_res_filt = W_nm - Z[:, :filter_orders].dot(coeffs[:filter_orders])
        W_res_filt_j1to3 = W_nm - Z[:, :3].dot(coeffs[:3])

        return {
            'coefficients': coeffs,
            'R': R,
            'global_rms': float(np.sqrt(np.mean(W_nm**2))),
            'filtered_rms': float(np.sqrt(np.mean(W_res_filt**2))),
            'rms_filter_j1to3': float(np.sqrt(np.mean(W_res_filt_j1to3**2))),
            'W_res_filt': W_res_filt,
        }

    def _compute_aberration_magnitudes(self, coeffs: np.ndarray) -> Dict[str, float]:
        """Compute magnitude of specific aberration modes."""
        return {
            'defocus_nm': float(abs(coeffs[3])) if len(coeffs) > 3 else 0.0,
            'astigmatism_rms': float(np.sqrt(coeffs[4]**2 + coeffs[5]**2)) if len(coeffs) > 5 else 0.0,
            'coma_rms': float(np.sqrt(coeffs[6]**2 + coeffs[7]**2)) if len(coeffs) > 7 else 0.0,
            'trefoil_rms': float(np.sqrt(coeffs[8]**2 + coeffs[9]**2)) if len(coeffs) > 9 else 0.0,
            'spherical_nm': float(abs(coeffs[10])) if len(coeffs) > 10 else 0.0,
        }

    def _generate_view_html(
        self,
        title: str,
        X: np.ndarray,
        Y: np.ndarray,
        W_nm: np.ndarray,
        rms_data: Dict,
        config: Dict,
        is_relative: bool = False,
        ref_title: str = "20 deg",
        abs_pair: Optional[Tuple[float, float]] = None,
        is_manufacturing: bool = False,
        mfg_metrics: Optional[Dict] = None,
        correction_metrics: Optional[Dict] = None,
    ) -> str:
        """Generate HTML for a single view."""
        _load_dependencies()

        n_modes = config.get('n_modes', 50)
        amp = config.get('amp', 0.5)
        pancake = config.get('pancake', 3.0)
        downsample = config.get('plot_downsample', 10000)
        colorscale = config.get('colorscale', 'Turbo')
        show_bar = config.get('show_bar_chart', True)

        coeffs = rms_data['coefficients']
        global_rms = rms_data['global_rms']
        filtered_rms = rms_data['filtered_rms']
        W_res_filt = rms_data['W_res_filt']

        labels = [zernike_label(j) for j in range(1, n_modes + 1)]
        coeff_abs = np.abs(coeffs)

        # Downsample
        n = len(X)
        if n > downsample:
            rng = np.random.default_rng(42)
            sel = rng.choice(n, size=downsample, replace=False)
            Xp, Yp, Wp = X[sel], Y[sel], W_res_filt[sel]
        else:
            Xp, Yp, Wp = X, Y, W_res_filt

        res_amp = amp * Wp
        max_amp = float(np.max(np.abs(res_amp))) if res_amp.size else 1.0

        # Build mesh
        mesh_traces = []
        try:
            tri = _Triangulation(Xp, Yp)
            if tri.triangles is not None and len(tri.triangles) > 0:
                i, j, k = tri.triangles.T
                mesh_traces.append(_go.Mesh3d(
                    x=Xp, y=Yp, z=res_amp,
                    i=i, j=j, k=k,
                    intensity=res_amp,
                    colorscale=colorscale,
                    opacity=1.0,
                    flatshading=False,
                    lighting=dict(ambient=0.4, diffuse=0.8, specular=0.3,
                                  roughness=0.5, fresnel=0.2),
                    lightposition=dict(x=100, y=200, z=300),
                    showscale=True,
                    colorbar=dict(title=dict(text="Residual (nm)", side="right"),
                                  thickness=15, len=0.6, tickformat=".1f"),
                    hovertemplate="X: %{x:.1f}<br>Y: %{y:.1f}<br>Residual: %{z:.2f} nm<extra></extra>"
                ))
        except Exception:
            pass

        if not mesh_traces:
            mesh_traces.append(_go.Scatter3d(
                x=Xp, y=Yp, z=res_amp,
                mode='markers',
                marker=dict(size=2, color=res_amp, colorscale=colorscale, showscale=True),
                showlegend=False
            ))

        title_suffix = f" (relative to {ref_title})" if is_relative else " (absolute)"

        # Build subplots
        if is_manufacturing and mfg_metrics and correction_metrics:
            fig = _make_subplots(
                rows=5, cols=1,
                specs=[[{"type": "scene"}], [{"type": "table"}], [{"type": "table"}],
                       [{"type": "table"}], [{"type": "xy"}]],
                row_heights=[0.38, 0.12, 0.12, 0.18, 0.20],
                vertical_spacing=0.025,
                subplot_titles=[
                    f"<b>Surface Residual - {title}{title_suffix}</b>",
                    "<b>RMS Metrics (Absolute 90 deg)</b>",
                    "<b>Mode Magnitudes at 90 deg</b>",
                    "<b>Pre-Correction (90 deg - 20 deg)</b>",
                    "<b>|Zernike Coefficients| (nm)</b>"
                ]
            )
        elif show_bar:
            fig = _make_subplots(
                rows=4, cols=1,
                specs=[[{"type": "scene"}], [{"type": "table"}],
                       [{"type": "table"}], [{"type": "xy"}]],
                row_heights=[0.45, 0.12, 0.25, 0.18],
                vertical_spacing=0.03,
                subplot_titles=[
                    f"<b>Surface Residual - {title}{title_suffix}</b>",
                    "<b>RMS Metrics</b>",
                    f"<b>Zernike Coefficients ({n_modes} modes)</b>",
                    "<b>|Zernike Coefficients| (nm)</b>"
                ]
            )
        else:
            fig = _make_subplots(
                rows=3, cols=1,
                specs=[[{"type": "scene"}], [{"type": "table"}], [{"type": "table"}]],
                row_heights=[0.55, 0.15, 0.30],
                vertical_spacing=0.03,
                subplot_titles=[
                    f"<b>Surface Residual - {title}{title_suffix}</b>",
                    "<b>RMS Metrics</b>",
                    f"<b>Zernike Coefficients ({n_modes} modes)</b>"
                ]
            )

        # Add mesh
        for tr in mesh_traces:
            fig.add_trace(tr, row=1, col=1)

        # Configure 3D scene
        fig.update_scenes(
            camera=dict(eye=dict(x=1.2, y=1.2, z=0.8), up=dict(x=0, y=0, z=1)),
            xaxis=dict(title="X (mm)", showgrid=True,
                       gridcolor='rgba(128,128,128,0.3)',
                       showbackground=True, backgroundcolor='rgba(240,240,240,0.9)'),
            yaxis=dict(title="Y (mm)", showgrid=True,
                       gridcolor='rgba(128,128,128,0.3)',
                       showbackground=True, backgroundcolor='rgba(240,240,240,0.9)'),
            zaxis=dict(title="Residual (nm)",
                       range=[-max_amp * pancake, max_amp * pancake],
                       showgrid=True, gridcolor='rgba(128,128,128,0.3)',
                       showbackground=True, backgroundcolor='rgba(230,230,250,0.9)'),
            aspectmode='manual',
            aspectratio=dict(x=1, y=1, z=0.4)
        )

        # Add tables
        if is_relative and abs_pair:
            abs_global, abs_filtered = abs_pair
            fig.add_trace(_go.Table(
                header=dict(values=["<b>Metric</b>", "<b>Relative (nm)</b>", "<b>Absolute (nm)</b>"],
                            align="left", fill_color='#1f2937', font=dict(color='white')),
                cells=dict(values=[
                    ["Global RMS", "Filtered RMS (J1-J4 removed)"],
                    [f"{global_rms:.2f}", f"{filtered_rms:.2f}"],
                    [f"{abs_global:.2f}", f"{abs_filtered:.2f}"],
                ], align="left", fill_color='#374151', font=dict(color='white'))
            ), row=2, col=1)
        elif is_manufacturing and mfg_metrics and correction_metrics:
            fig.add_trace(_go.Table(
                header=dict(values=["<b>Metric</b>", "<b>Value (nm)</b>"],
                            align="left", fill_color='#1f2937', font=dict(color='white')),
                cells=dict(values=[
                    ["Global RMS", "Filtered RMS (J1-J4)"],
                    [f"{global_rms:.2f}", f"{filtered_rms:.2f}"]
                ], align="left", fill_color='#374151', font=dict(color='white'))
            ), row=2, col=1)

            fig.add_trace(_go.Table(
                header=dict(values=["<b>Mode</b>", "<b>Value (nm)</b>"],
                            align="left", fill_color='#1f2937', font=dict(color='white')),
                cells=dict(values=[
                    ["Filtered RMS (J1-J3, with defocus)", "Astigmatism (J5+J6)",
                     "Coma (J7+J8)", "Trefoil (J9+J10)", "Spherical (J11)"],
                    [f"{rms_data['rms_filter_j1to3']:.2f}",
                     f"{mfg_metrics['astigmatism_rms']:.2f}",
                     f"{mfg_metrics['coma_rms']:.2f}",
                     f"{mfg_metrics['trefoil_rms']:.2f}",
                     f"{mfg_metrics['spherical_nm']:.2f}"]
                ], align="left", fill_color='#374151', font=dict(color='white'))
            ), row=3, col=1)

            fig.add_trace(_go.Table(
                header=dict(values=["<b>Mode</b>", "<b>Correction (nm)</b>"],
                            align="left", fill_color='#1f2937', font=dict(color='white')),
                cells=dict(values=[
                    ["Total RMS (J1-J3 filter)", "Defocus (J4)",
                     "Astigmatism (J5+J6)", "Coma (J7+J8)"],
                    [f"{correction_metrics['rms_filter_j1to3']:.2f}",
                     f"{correction_metrics['defocus_nm']:.2f}",
                     f"{correction_metrics['astigmatism_rms']:.2f}",
                     f"{correction_metrics['coma_rms']:.2f}"]
                ], align="left", fill_color='#374151', font=dict(color='white'))
            ), row=4, col=1)
        else:
            fig.add_trace(_go.Table(
                header=dict(values=["<b>Metric</b>", "<b>Value (nm)</b>"],
                            align="left", fill_color='#1f2937', font=dict(color='white')),
                cells=dict(values=[
                    ["Global RMS", "Filtered RMS (J1-J4 removed)"],
                    [f"{global_rms:.2f}", f"{filtered_rms:.2f}"]
                ], align="left", fill_color='#374151', font=dict(color='white'))
            ), row=2, col=1)

        # Coefficients table
        if not (is_manufacturing and mfg_metrics and correction_metrics):
            fig.add_trace(_go.Table(
                header=dict(values=["<b>Noll j</b>", "<b>Label</b>", "<b>|Coeff| (nm)</b>"],
                            align="left", fill_color='#1f2937', font=dict(color='white')),
                cells=dict(values=[
                    list(range(1, n_modes + 1)),
                    labels,
                    [f"{c:.3f}" for c in coeff_abs]
                ], align="left", fill_color='#374151', font=dict(color='white'))
            ), row=3, col=1)

        # Bar chart
        if show_bar:
            bar_row = 5 if (is_manufacturing and mfg_metrics and correction_metrics) else 4
            fig.add_trace(
                _go.Bar(
                    x=coeff_abs.tolist(), y=labels,
                    orientation='h', marker_color='#6366f1',
                    hovertemplate="%{y}<br>|Coeff| = %{x:.3f} nm<extra></extra>",
                    showlegend=False
                ),
                row=bar_row, col=1
            )

        # Layout
        height = 1500 if (is_manufacturing and mfg_metrics and correction_metrics) else 1300
        fig.update_layout(
            width=1400, height=height,
            margin=dict(t=60, b=20, l=20, r=20),
            paper_bgcolor='#111827', plot_bgcolor='#1f2937',
            font=dict(color='white'),
            title=dict(text=f"<b>Atomizer Zernike Analysis - {title}</b>",
                       x=0.5, font=dict(size=18))
        )

        return fig.to_html(include_plotlyjs='cdn', full_html=True)

    def _generate(self, config: InsightConfig) -> InsightResult:
        """Generate all Zernike WFE views."""
        self._load_data()

        # Merge config
        cfg = {**DEFAULT_CONFIG, **config.extra}
        cfg['colorscale'] = config.extra.get('colorscale', cfg['colorscale'])
        cfg['amp'] = config.amplification if config.amplification != 1.0 else cfg['amp']

        n_modes = cfg['n_modes']
        filter_orders = cfg['filter_low_orders']
        disp_unit = cfg['disp_unit']

        # Map subcases
        disps = self._displacements
        if '1' in disps and '2' in disps:
            sc_map = {'90': '1', '20': '2', '40': '3', '60': '4'}
        elif '90' in disps and '20' in disps:
            sc_map = {'90': '90', '20': '20', '40': '40', '60': '60'}
        else:
            available = sorted(disps.keys(), key=lambda x: int(x) if x.isdigit() else 0)
            if len(available) >= 4:
                sc_map = {'90': available[0], '20': available[1],
                          '40': available[2], '60': available[3]}
            else:
                return InsightResult(success=False,
                                     error=f"Need 4 subcases, found: {available}")

        # Check subcases
        for angle, label in sc_map.items():
            if label not in disps:
                return InsightResult(success=False,
                                     error=f"Subcase '{label}' (angle {angle}) not found")

        timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
        output_dir = config.output_dir or self.insights_path
        output_dir.mkdir(parents=True, exist_ok=True)

        html_files = []
        summary = {}

        # Reference: 20 deg
        X_ref, Y_ref, WFE_ref, nids_ref = self._build_wfe_arrays(sc_map['20'], disp_unit)
        rms_ref = self._compute_metrics(X_ref, Y_ref, WFE_ref, n_modes, filter_orders)

        # 90 deg
        X_90, Y_90, WFE_90, nids_90 = self._build_wfe_arrays(sc_map['90'], disp_unit)
        rms_90 = self._compute_metrics(X_90, Y_90, WFE_90, n_modes, filter_orders)
        mfg_metrics = self._compute_aberration_magnitudes(rms_90['coefficients'])

        # 40 deg vs 20 deg
        X_40, Y_40, WFE_40, nids_40 = self._build_wfe_arrays(sc_map['40'], disp_unit)
        X_40_rel, Y_40_rel, WFE_40_rel = self._compute_relative_wfe(
            X_40, Y_40, WFE_40, nids_40, X_ref, Y_ref, WFE_ref, nids_ref)
        rms_40_abs = self._compute_metrics(X_40, Y_40, WFE_40, n_modes, filter_orders)
        rms_40_rel = self._compute_metrics(X_40_rel, Y_40_rel, WFE_40_rel, n_modes, filter_orders)

        html_40 = self._generate_view_html(
            "40 deg", X_40_rel, Y_40_rel, WFE_40_rel, rms_40_rel, cfg,
            is_relative=True, ref_title="20 deg",
            abs_pair=(rms_40_abs['global_rms'], rms_40_abs['filtered_rms']))
        path_40 = output_dir / f"zernike_{timestamp}_40_vs_20.html"
        path_40.write_text(html_40, encoding='utf-8')
        html_files.append(path_40)
        summary['40_vs_20_filtered_rms'] = rms_40_rel['filtered_rms']

        # 60 deg vs 20 deg
        X_60, Y_60, WFE_60, nids_60 = self._build_wfe_arrays(sc_map['60'], disp_unit)
        X_60_rel, Y_60_rel, WFE_60_rel = self._compute_relative_wfe(
            X_60, Y_60, WFE_60, nids_60, X_ref, Y_ref, WFE_ref, nids_ref)
        rms_60_abs = self._compute_metrics(X_60, Y_60, WFE_60, n_modes, filter_orders)
        rms_60_rel = self._compute_metrics(X_60_rel, Y_60_rel, WFE_60_rel, n_modes, filter_orders)

        html_60 = self._generate_view_html(
            "60 deg", X_60_rel, Y_60_rel, WFE_60_rel, rms_60_rel, cfg,
            is_relative=True, ref_title="20 deg",
            abs_pair=(rms_60_abs['global_rms'], rms_60_abs['filtered_rms']))
        path_60 = output_dir / f"zernike_{timestamp}_60_vs_20.html"
        path_60.write_text(html_60, encoding='utf-8')
        html_files.append(path_60)
        summary['60_vs_20_filtered_rms'] = rms_60_rel['filtered_rms']

        # 90 deg Manufacturing
        X_90_rel, Y_90_rel, WFE_90_rel = self._compute_relative_wfe(
            X_90, Y_90, WFE_90, nids_90, X_ref, Y_ref, WFE_ref, nids_ref)
        rms_90_rel = self._compute_metrics(X_90_rel, Y_90_rel, WFE_90_rel, n_modes, filter_orders)
        corr_abr = self._compute_aberration_magnitudes(rms_90_rel['coefficients'])
        correction_metrics = {
            'rms_filter_j1to3': rms_90_rel['rms_filter_j1to3'],
            **corr_abr
        }

        html_90 = self._generate_view_html(
            "90 deg (Manufacturing)", X_90, Y_90, WFE_90, rms_90, cfg,
            is_relative=False, is_manufacturing=True,
            mfg_metrics=mfg_metrics, correction_metrics=correction_metrics)
        path_90 = output_dir / f"zernike_{timestamp}_90_mfg.html"
        path_90.write_text(html_90, encoding='utf-8')
        html_files.append(path_90)
        summary['90_mfg_filtered_rms'] = rms_90['filtered_rms']
        summary['90_optician_workload'] = rms_90['rms_filter_j1to3']

        return InsightResult(
            success=True,
            html_path=html_files[0],  # Return first as primary
            summary={
                'html_files': [str(p) for p in html_files],
                **summary
            }
        )