feat(report): replace LSF/MSF with Tony Hull PSD analysis

- Remove compute_spatial_freq_metrics() and _spatial_freq_html() - Add compute_surface_psd(): 2D FFT + Hann window + radial averaging - Add compute_psd_band_rms(): gravity/support/HF band decomposition - Add make_psd_plot(): interactive log-log PSD plot with band annotations - Add _psd_summary_html(): band RMS cards with % breakdown - New section in report: Power Spectral Density Analysis - Zernike details now show only coefficient bar charts (cleaner) - Methodology: Tony Hull JWST approach for WFE spatial frequency analysis
2026-01-29 22:15:42 +00:00
parent 487ecf67dc
commit eeacfbe41a
1 changed files with 170 additions and 74 deletions
--- a/tools/generate_optical_report.py
+++ b/tools/generate_optical_report.py
@@ -44,6 +44,8 @@ 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
@@ -277,101 +279,159 @@ def aberration_magnitudes(coeffs):
    }


-def compute_spatial_freq_metrics(coefficients):
-    """Compute spatial frequency band metrics from Zernike coefficients.
+def compute_surface_psd(X, Y, Z, aperture_radius):
+    """Compute power spectral density of surface height data (Tony Hull methodology).

-    Decomposes the Zernike spectrum into Low Spatial Frequency (form/figure)
-    and Mid Spatial Frequency (ripple) bands based on radial order.
+    Interpolates scattered FEA node data onto a uniform grid, applies a Hann
+    window, computes the 2D FFT, and radially averages to produce a 1D PSD.

    Args:
-        coefficients: array of Zernike coefficients (50 modes, Noll convention)
+        X: Surface x-coordinates [mm]
+        Y: Surface y-coordinates [mm]
+        Z: Surface height values [nm]
+        aperture_radius: Physical radius of the mirror [mm]

    Returns:
-        dict with band RMS values and per-radial-order breakdown
+        (freqs, psd_values) — spatial frequencies in cycles/aperture and PSD amplitude
    """
-    n_modes = min(len(coefficients), 50)
-    coeffs = np.array(coefficients[:n_modes])
+    Xc = X - np.mean(X)
+    Yc = Y - np.mean(Y)

-    # Build radial order mapping using noll_indices
-    order_map = {}  # n -> list of coeff values
-    for j in range(1, n_modes + 1):
-        n, m = noll_indices(j)
-        if n not in order_map:
-            order_map[n] = []
-        order_map[n].append(coeffs[j - 1])
+    grid_size = 256 if len(X) >= 1000 else min(128, int(np.sqrt(len(X))))

-    # Band definitions (by Noll index, 1-based)
-    # LSF: J4-J15 (n=2..4) — Low Spatial Frequency (form/figure)
-    # MSF: J16-J50 (n=5..9) — Mid Spatial Frequency (ripple)
-    lsf_coeffs = coeffs[3:15]              # indices 3..14 → J4..J15
-    msf_coeffs = coeffs[15:n_modes]        # indices 15..49 → J16..J50
-    total_coeffs = coeffs[3:n_modes]       # J4..J50 (excluding piston/tilt)
+    x_range = np.linspace(Xc.min(), Xc.max(), grid_size)
+    y_range = np.linspace(Yc.min(), Yc.max(), grid_size)
+    X_grid, Y_grid = np.meshgrid(x_range, y_range)

-    lsf_rms = float(np.sqrt(np.sum(lsf_coeffs**2)))
-    msf_rms = float(np.sqrt(np.sum(msf_coeffs**2)))
-    total_filtered_rms = float(np.sqrt(np.sum(total_coeffs**2)))
-    lsf_msf_ratio = lsf_rms / msf_rms if msf_rms > 0 else float('inf')
+    Z_grid = griddata((Xc, Yc), Z, (X_grid, Y_grid), method='cubic')
+    Z_grid = np.nan_to_num(Z_grid, nan=0.0)

-    # Per-radial-order RMS
-    per_order = {}
-    for n in sorted(order_map.keys()):
-        order_coeffs = np.array(order_map[n])
-        per_order[n] = float(np.sqrt(np.sum(order_coeffs**2)))
+    # Circular aperture mask
+    R_grid = np.sqrt(X_grid**2 + Y_grid**2)
+    Z_grid[R_grid > aperture_radius] = 0.0

+    # Hann window to reduce spectral leakage
+    hann = np.outer(np.hanning(grid_size), np.hanning(grid_size))
+    Z_windowed = Z_grid * hann
+
+    fft_result = fft2(Z_windowed)
+    dx = x_range[1] - x_range[0]
+    dy = y_range[1] - y_range[0]
+    freqs_x = fftfreq(grid_size, dx)
+    freqs_y = fftfreq(grid_size, dy)
+    fy, fx = np.meshgrid(freqs_y, freqs_x)
+    freqs_radial = np.sqrt(fx**2 + fy**2) * 2 * aperture_radius  # cycles/aperture
+
+    power_spectrum = np.abs(fft_result)**2
+
+    bin_edges = np.logspace(-1.5, np.log10(grid_size / 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(power_spectrum[mask]) * dx * dy)
+            if avg > 0:
+                psd_values.append(avg)
+                freqs_center.append(float(np.sqrt(bin_edges[i] * bin_edges[i + 1])))
+
+    return np.array(freqs_center), np.array(psd_values)
+
+
+def compute_psd_band_rms(freqs, psd):
+    """Compute integrated RMS for gravity, support, and HF bands."""
+    def _band_rms(lo, hi):
+        m = (freqs >= lo) & (freqs <= hi)
+        if not np.any(m):
+            return 0.0
+        return float(np.sqrt(np.trapz(psd[m], freqs[m])))
+
+    gravity = _band_rms(0.1, 2.0)
+    support = _band_rms(2.0, 20.0)
+    hf = _band_rms(20.0, freqs.max())
+    total = float(np.sqrt(np.trapz(psd, freqs)))
    return {
-        'lsf_rms': lsf_rms,
-        'msf_rms': msf_rms,
-        'total_filtered_rms': total_filtered_rms,
-        'lsf_msf_ratio': lsf_msf_ratio,
-        'per_order': per_order,
+        'gravity_rms': gravity,
+        'support_rms': support,
+        'hf_rms': hf,
+        'total_rms': total,
    }


-def _spatial_freq_html(metrics):
-    """Generate HTML card for spatial frequency band metrics."""
-    m = metrics
+def make_psd_plot(psd_data_dict, title="Power Spectral Density"):
+    """Create a log-log PSD plot for multiple angles.

-    # Per-order breakdown for n=2..9
-    order_items = ""
-    for n in range(2, 10):
-        val = m['per_order'].get(n, 0.0)
-        order_items += (
-            f'<div class="sf-order-item">'
-            f'<span class="sf-order-n">n={n}:</span>'
-            f'<span class="sf-order-val">{val:.2f}</span>'
-            f'</div>\n'
-        )
+    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()

-    ratio_str = (
-        f"{m['lsf_msf_ratio']:.1f}\u00d7"
-        if m['lsf_msf_ratio'] < 1000
-        else "\u221e"
+    for label, (freqs, psd) in psd_data_dict.items():
+        valid = psd > 0
+        fig.add_trace(go.Scatter(
+            x=freqs[valid], y=psd[valid],
+            mode='lines', name=label,
+            line=dict(color=colors.get(label, '#6366f1'), width=2),
+            hovertemplate="%{x:.1f} cyc/apt: %{y:.2e}<extra>" + label + "</extra>",
+        ))
+
+    # Band annotations
+    fig.add_vrect(x0=0.1, x1=2.0, fillcolor='rgba(59,130,246,0.08)', line_width=0, layer='below',
+                  annotation_text="Gravity", annotation_position="top left",
+                  annotation=dict(font=dict(size=10, color='#3b82f6')))
+    fig.add_vrect(x0=2.0, x1=20.0, fillcolor='rgba(245,158,11,0.08)', line_width=0, layer='below',
+                  annotation_text="Support", annotation_position="top left",
+                  annotation=dict(font=dict(size=10, color='#f59e0b')))
+    fig.add_vrect(x0=20.0, x1=200.0, fillcolor='rgba(239,68,68,0.06)', line_width=0, layer='below',
+                  annotation_text="High Freq", annotation_position="top left",
+                  annotation=dict(font=dict(size=10, 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')),
+                   gridcolor='#e2e8f0', tickfont=dict(color='#475569')),
+        yaxis=dict(type='log', title=dict(text="PSD [nm²·mm²]", font=dict(color='#1e293b')),
+                   gridcolor='#e2e8f0', tickfont=dict(color='#475569')),
+        legend=dict(x=0.01, y=0.99, bgcolor='rgba(255,255,255,0.9)',
+                    bordercolor='#e2e8f0', borderwidth=1, font=dict(size=11, color='#1e293b')),
    )
+    return fig.to_html(include_plotlyjs=False, full_html=False, div_id="psd_plot")
+
+
+def _psd_summary_html(band_dict):
+    """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>Spatial Frequency Breakdown</h4>
+        <h4>PSD Band Decomposition (Tony Hull Methodology)</h4>
        <div class="sf-band-grid">
            <div class="sf-band-card sf-lsf">
-                <div class="sf-band-label">LSF (Form)</div>
-                <div class="sf-band-range">J4\u2013J15 &nbsp; n\u22644</div>
-                <div class="sf-band-value">{m['lsf_rms']:.2f} <span class="sf-unit">nm</span></div>
+                <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">MSF (Ripple)</div>
-                <div class="sf-band-range">J16\u2013J50 &nbsp; n\u22655</div>
-                <div class="sf-band-value">{m['msf_rms']:.2f} <span class="sf-unit">nm</span></div>
+                <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">LSF / MSF</div>
-                <div class="sf-band-range">Ratio</div>
-                <div class="sf-band-value">{ratio_str}</div>
+                <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 class="sf-order-grid">
-            <div class="sf-order-title">Per-Order RMS (nm)</div>
-            <div class="sf-order-items">{order_items}</div>
-        </div>
    </div>
    """

@@ -870,10 +930,36 @@ def generate_report(
    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")

-    # Spatial frequency band metrics
-    sfm_40 = _spatial_freq_html(compute_spatial_freq_metrics(angle_results[40]['rms_rel']['coefficients']))
-    sfm_60 = _spatial_freq_html(compute_spatial_freq_metrics(angle_results[60]['rms_rel']['coefficients']))
-    sfm_90 = _spatial_freq_html(compute_spatial_freq_metrics(angle_results[90]['rms_abs']['coefficients']))
+    # 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 = {}
+    for ang_label, ang_key, wfe_key in [
+        ('40° vs 20°', 40, 'WFE_rel'),
+        ('60° vs 20°', 60, 'WFE_rel'),
+        ('90° (Abs)', 90, 'WFE'),
+    ]:
+        r = angle_results[ang_key]
+        Xp = r['X_rel'] if 'rel' in wfe_key else r['X']
+        Yp = r['Y_rel'] if 'rel' in wfe_key else r['Y']
+        Zp = r[wfe_key] if wfe_key == 'WFE' else r['WFE_rel']
+        try:
+            freqs, psd = compute_surface_psd(Xp, Yp, Zp, R_outer_est)
+            psd_plot_data[ang_label] = (freqs, psd)
+            psd_bands[ang_label] = compute_psd_band_rms(freqs, psd)
+            print(f"  {ang_label}: gravity={psd_bands[ang_label]['gravity_rms']:.2f} nm, "
+                  f"support={psd_bands[ang_label]['support_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°']) if '40° vs 20°' in psd_bands else ""
+    psd_summary_60 = _psd_summary_html(psd_bands['60° vs 20°']) if '60° vs 20°' in psd_bands else ""
+    psd_summary_90 = _psd_summary_html(psd_bands['90° (Abs)']) if '90° (Abs)' in psd_bands else ""

    # Per-angle RMS plot
    angle_rms_data = {}
@@ -982,7 +1068,8 @@ def generate_report(
    # Section numbering: adjust if trajectory present
    sec_surface = 3
    sec_traj = 4
-    sec_mfg = 5 if traj_result else 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
@@ -1370,6 +1457,15 @@ def generate_report(
 '<div class="plot-container"><h3>Axial vs Lateral Sensitivity</h3>' + sens_plot_html + '</div>' +
 '</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>
@@ -1388,15 +1484,15 @@ def generate_report(
    <h2>{sec_zernike}. Zernike Coefficient Details</h2>
    <details>
        <summary>40\u00b0 vs 20\u00b0 \u2014 Relative Coefficients</summary>
-        <div>{sfm_40}{bar_40}</div>
+        <div>{bar_40}</div>
    </details>
    <details>
        <summary>60\u00b0 vs 20\u00b0 \u2014 Relative Coefficients</summary>
-        <div>{sfm_60}{bar_60}</div>
+        <div>{bar_60}</div>
    </details>
    <details>
        <summary>90\u00b0 \u2014 Absolute Coefficients</summary>
-        <div>{sfm_90}{bar_90}</div>
+        <div>{bar_90}</div>
    </details>
 </div>