diff --git a/tools/generate_optical_report.py b/tools/generate_optical_report.py index 0ec14eba..0727d5c5 100644 --- a/tools/generate_optical_report.py +++ b/tools/generate_optical_report.py @@ -280,27 +280,27 @@ def aberration_magnitudes(coeffs): def compute_surface_psd(X, Y, Z, aperture_radius): - """Compute power spectral density of surface height data (Tony Hull methodology). + """Compute PSD and band RMS of surface height data (Tony Hull methodology). - 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. + 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: Surface x-coordinates [mm] - Y: Surface y-coordinates [mm] - Z: Surface height values [nm] - aperture_radius: Physical radius of the mirror [mm] + X, Y: coordinates [mm] + Z: surface height [nm] + aperture_radius: mirror radius [mm] Returns: - (freqs, psd_values) — spatial frequencies in cycles/aperture and PSD amplitude + dict with 'freqs', 'psd' (for plotting), 'bands' (gravity/support/hf/total RMS in nm) """ Xc = X - np.mean(X) Yc = Y - np.mean(Y) - grid_size = 256 if len(X) >= 1000 else min(128, int(np.sqrt(len(X)))) + N = 256 if len(X) >= 1000 else min(128, int(np.sqrt(len(X)))) - x_range = np.linspace(Xc.min(), Xc.max(), grid_size) - y_range = np.linspace(Yc.min(), Yc.max(), grid_size) + 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') @@ -308,52 +308,66 @@ def compute_surface_psd(X, Y, Z, aperture_radius): # Circular aperture mask R_grid = np.sqrt(X_grid**2 + Y_grid**2) - Z_grid[R_grid > aperture_radius] = 0.0 + apt_mask = R_grid <= aperture_radius + Z_grid[~apt_mask] = 0.0 - # Hann window to reduce spectral leakage - hann = np.outer(np.hanning(grid_size), np.hanning(grid_size)) + # 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(grid_size, dx) - freqs_y = fftfreq(grid_size, dy) + freqs_x = fftfreq(N, dx) + freqs_y = fftfreq(N, dy) fy, fx = np.meshgrid(freqs_y, freqs_x) - freqs_radial = np.sqrt(fx**2 + fy**2) * 2 * aperture_radius # cycles/aperture + D = 2 * aperture_radius + freqs_radial = np.sqrt(fx**2 + fy**2) * D # cycles/aperture - power_spectrum = np.abs(fft_result)**2 + # ── 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 - bin_edges = np.logspace(-1.5, np.log10(grid_size / 2), 60) + 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(power_spectrum[mask]) * dx * dy) + 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 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 { - 'gravity_rms': gravity, - 'support_rms': support, - 'hf_rms': hf, - 'total_rms': total, + 'freqs': np.array(freqs_center), + 'psd': np.array(psd_values), + 'bands': bands, } @@ -938,21 +952,29 @@ def generate_report( 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'), + # 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] - 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'] + 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: - 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") + 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}")