diff --git a/tools/generate_optical_report.py b/tools/generate_optical_report.py index 957fbce6..0ec14eba 100644 --- 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'
' - f'n={n}:' - f'{val:.2f}' - f'
\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}" + label + "", + )) + + # 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"""
-

Spatial Frequency Breakdown

+

PSD Band Decomposition (Tony Hull Methodology)

-
LSF (Form)
-
J4\u2013J15   n\u22644
-
{m['lsf_rms']:.2f} nm
+
Gravity Signature
+
0.1\u20132 cyc/apt
+
{m['gravity_rms']:.2f} nm
+
{grav_pct:.0f}% of total
-
MSF (Ripple)
-
J16\u2013J50   n\u22655
-
{m['msf_rms']:.2f} nm
+
Support Print-through
+
2\u201320 cyc/apt
+
{m['support_rms']:.2f} nm
+
{supp_pct:.0f}% of total
-
LSF / MSF
-
Ratio
-
{ratio_str}
+
High Frequency
+
>20 cyc/apt
+
{m['hf_rms']:.2f} nm
+
{hf_pct:.0f}% of total
-
-
Per-Order RMS (nm)
-
{order_items}
-
""" @@ -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( '

Axial vs Lateral Sensitivity

' + sens_plot_html + '
' + '' if traj_result else ''} + +{'

' + str(sec_psd) + '. Power Spectral Density Analysis

' + + '

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 (>20 cyc/apt).

' + + '

PSD \u2014 Log-Log

' + psd_plot_html + '
' + + psd_summary_40 + psd_summary_60 + psd_summary_90 + + '
' if psd_plot_html else ''} +

{sec_mfg}. Manufacturing Analysis (90\u00b0 Orientation)

@@ -1388,15 +1484,15 @@ def generate_report(

{sec_zernike}. Zernike Coefficient Details

40\u00b0 vs 20\u00b0 \u2014 Relative Coefficients -
{sfm_40}{bar_40}
+
{bar_40}
60\u00b0 vs 20\u00b0 \u2014 Relative Coefficients -
{sfm_60}{bar_60}
+
{bar_60}
90\u00b0 \u2014 Absolute Coefficients -
{sfm_90}{bar_90}
+
{bar_90}