Files
Atomizer/tools/generate_optical_report.py
Antoine 34b52f9543 Add comprehensive optical performance report generator
New tool: tools/generate_optical_report.py
- CDR-ready single HTML report from OP2 results
- Executive summary with pass/fail vs targets
- Per-angle WFE analysis with 3D surface plots
- Zernike trajectory analysis (mode-specific RMS)
- Axial vs lateral sensitivity matrix
- Manufacturing correction metrics
- Collapsible Zernike coefficient bar charts
- Optional study DB integration for design params
- Annular aperture support (default M1 inner R=135.75mm)
- Dark theme, interactive Plotly charts, print-friendly
2026-01-29 18:28:10 +00:00

1224 lines
45 KiB
Python
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
#!/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
# 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,
}
# ============================================================================
# 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,
}
# ============================================================================
# HTML Report Generation
# ============================================================================
def status_badge(value, target, unit='nm'):
"""Return pass/fail badge HTML."""
if value <= target:
return f'<span class="badge pass">✅ {value:.2f} {unit}{target:.1f}</span>'
ratio = value / target
if ratio < 1.5:
return f'<span class="badge warn">⚠️ {value:.2f} {unit} ({ratio:.1f}× target)</span>'
return f'<span class="badge fail">❌ {value:.2f} {unit} ({ratio:.1f}× target)</span>'
def make_surface_plot(X, Y, W_res, mask, inner_radius=None, title="", amp=0.5, downsample=PLOT_DOWNSAMPLE):
"""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"), thickness=12, len=0.5, tickformat=".1f"),
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='white', 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="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*3.0, max_amp*3.0],
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),
),
margin=dict(t=30, b=10, l=10, r=10),
paper_bgcolor='rgba(0,0,0,0)',
plot_bgcolor='rgba(0,0,0,0)',
font=dict(color='#e0e0e0'),
height=500,
width=700,
)
return fig.to_html(include_plotlyjs=False, full_html=False, div_id=f"surface_{title.replace(' ','_')}")
def make_bar_chart(coeffs, title="Zernike Coefficients", max_modes=30):
"""Create horizontal bar chart of Zernike coefficient magnitudes."""
n = min(len(coeffs), max_modes)
labels = [zernike_label(j) for j in range(1, n+1)]
vals = np.abs(coeffs[:n])
fig = go.Figure(go.Bar(
x=vals, y=labels, orientation='h',
marker_color='#6366f1',
hovertemplate="%{y}<br>|c| = %{x:.3f} nm<extra></extra>",
))
fig.update_layout(
height=max(400, n*22),
margin=dict(t=30, b=10, l=200, r=20),
paper_bgcolor='rgba(0,0,0,0)',
plot_bgcolor='rgba(17,24,39,0.8)',
font=dict(color='#e0e0e0', size=10),
xaxis=dict(title="|Coefficient| (nm)", gridcolor='rgba(128,128,128,0.2)'),
yaxis=dict(autorange='reversed'),
)
return fig.to_html(include_plotlyjs=False, full_html=False, div_id=f"bar_{title.replace(' ','_')}")
def make_trajectory_plot(angles, coefficients_relative, mode_groups, sensitivity, title=""):
"""Create trajectory visualization: Zernike modes vs elevation angle."""
fig = go.Figure()
# Plot each mode group
colors = ['#f59e0b', '#ef4444', '#10b981', '#6366f1', '#ec4899', '#14b8a6', '#f97316']
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
# RSS of modes in this group at each angle
rss = np.sqrt(np.sum(coefficients_relative[:, indices]**2, axis=1))
color = colors[color_idx % len(colors)]
fig.add_trace(go.Scatter(
x=angles, y=rss,
mode='lines+markers',
name=MODE_NAMES.get(group_name, group_name),
line=dict(color=color, width=2),
marker=dict(size=8),
hovertemplate=f"{group_name}<br>%{{x}}°: %{{y:.2f}} nm<extra></extra>"
))
color_idx += 1
fig.update_layout(
height=400,
margin=dict(t=30, b=40, l=60, r=20),
paper_bgcolor='rgba(0,0,0,0)',
plot_bgcolor='rgba(17,24,39,0.8)',
font=dict(color='#e0e0e0'),
xaxis=dict(title="Elevation Angle (°)", gridcolor='rgba(128,128,128,0.2)',
tickvals=angles, dtick=10),
yaxis=dict(title="RMS (nm)", gridcolor='rgba(128,128,128,0.2)'),
legend=dict(x=0.02, y=0.98, bgcolor='rgba(17,24,39,0.7)'),
)
return fig.to_html(include_plotlyjs=False, full_html=False, div_id="trajectory_plot")
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 = [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 θ)', marker_color='#f59e0b',
hovertemplate="%{y}<br>Axial: %{x:.3f} nm/unit<extra></extra>"
))
fig.add_trace(go.Bar(
y=labels, x=lateral, orientation='h',
name='Lateral (cos θ)', marker_color='#6366f1',
hovertemplate="%{y}<br>Lateral: %{x:.3f} nm/unit<extra></extra>"
))
fig.update_layout(
barmode='group',
height=max(300, len(modes)*40),
margin=dict(t=30, b=40, l=200, r=20),
paper_bgcolor='rgba(0,0,0,0)',
plot_bgcolor='rgba(17,24,39,0.8)',
font=dict(color='#e0e0e0', size=11),
xaxis=dict(title="Sensitivity (nm per load fraction)", gridcolor='rgba(128,128,128,0.2)'),
yaxis=dict(autorange='reversed'),
legend=dict(x=0.6, y=0.98, bgcolor='rgba(17,24,39,0.7)'),
)
return fig.to_html(include_plotlyjs=False, full_html=False, div_id="sensitivity_bar")
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}° vs {ref_angle}°" for a in angles]
fig = go.Figure(go.Bar(
x=labels, y=rms_vals,
marker_color=['#10b981' if v < 10 else '#f59e0b' if v < 20 else '#ef4444' for v in rms_vals],
text=[f"{v:.2f} nm" for v in rms_vals],
textposition='outside',
hovertemplate="%{x}: %{y:.2f} nm<extra></extra>"
))
fig.update_layout(
height=350,
margin=dict(t=30, b=40, l=60, r=20),
paper_bgcolor='rgba(0,0,0,0)',
plot_bgcolor='rgba(17,24,39,0.8)',
font=dict(color='#e0e0e0'),
yaxis=dict(title="Filtered RMS WFE (nm)", gridcolor='rgba(128,128,128,0.2)'),
)
return fig.to_html(include_plotlyjs=False, full_html=False, div_id="per_angle_rms")
# ============================================================================
# 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
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"
)
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")
# 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 = "°" 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>🔧 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}°</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²</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: {''.join(traj_result['mode_ranking'][:5])}</p>
"""
# Manufacturing details
mfg_html = f"""
<table class="data-table">
<thead><tr><th>Metric</th><th>Absolute 90°</th><th>Correction (90°20°)</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>J1J3 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>
"""
# 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>
<script src="https://cdn.plot.ly/plotly-2.27.0.min.js"></script>
<style>
:root {{
--bg-primary: #0f172a;
--bg-secondary: #1e293b;
--bg-card: #1e293b;
--text-primary: #f1f5f9;
--text-secondary: #94a3b8;
--accent: #6366f1;
--accent-hover: #818cf8;
--success: #10b981;
--warning: #f59e0b;
--danger: #ef4444;
--border: #334155;
}}
* {{ margin: 0; padding: 0; box-sizing: border-box; }}
body {{
font-family: 'Inter', -apple-system, BlinkMacSystemFont, 'Segoe UI', system-ui, sans-serif;
background: var(--bg-primary);
color: var(--text-primary);
line-height: 1.6;
}}
.container {{ max-width: 1400px; margin: 0 auto; padding: 2rem; }}
/* Header */
.header {{
background: linear-gradient(135deg, #1e293b 0%, #0f172a 100%);
border: 1px solid var(--border);
border-radius: 12px;
padding: 2rem 3rem;
margin-bottom: 2rem;
display: flex;
justify-content: space-between;
align-items: center;
}}
.header h1 {{ font-size: 1.8rem; font-weight: 700; }}
.header .subtitle {{ color: var(--text-secondary); font-size: 0.95rem; margin-top: 0.3rem; }}
.header .branding {{
text-align: right;
font-size: 0.85rem;
color: var(--text-secondary);
}}
.header .branding .logo {{ font-size: 1.4rem; font-weight: 700; color: var(--accent); }}
/* Sections */
.section {{
background: var(--bg-card);
border: 1px solid var(--border);
border-radius: 12px;
padding: 1.5rem 2rem;
margin-bottom: 1.5rem;
}}
.section h2 {{
font-size: 1.3rem;
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(300px, 1fr));
gap: 1rem;
margin: 1rem 0;
}}
.exec-card {{
background: var(--bg-primary);
border: 1px solid var(--border);
border-radius: 8px;
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 .target {{ font-size: 0.8rem; color: var(--text-secondary); margin-top: 0.3rem; }}
/* Badges */
.badge {{
display: inline-block;
padding: 0.2rem 0.6rem;
border-radius: 4px;
font-size: 0.85rem;
font-weight: 600;
}}
.badge.pass {{ background: rgba(16,185,129,0.15); color: var(--success); }}
.badge.warn {{ background: rgba(245,158,11,0.15); color: var(--warning); }}
.badge.fail {{ background: rgba(239,68,68,0.15); color: var(--danger); }}
/* 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-primary);
font-weight: 600;
font-size: 0.85rem;
text-transform: uppercase;
letter-spacing: 0.05em;
color: var(--text-secondary);
}}
.data-table tr.highlight td {{
background: rgba(99,102,241,0.08);
font-weight: 600;
}}
/* Metrics Grid */
.metrics-grid {{
display: grid;
grid-template-columns: repeat(auto-fit, minmax(180px, 1fr));
gap: 0.8rem;
margin: 1rem 0;
}}
.metric-card {{
background: var(--bg-primary);
border: 1px solid var(--border);
border-radius: 8px;
padding: 1rem;
text-align: center;
}}
.metric-label {{ font-size: 0.8rem; color: var(--text-secondary); margin-bottom: 0.3rem; }}
.metric-value {{ font-size: 1.3rem; font-weight: 700; color: var(--accent); }}
/* Plots */
.plot-grid {{
display: grid;
grid-template-columns: repeat(auto-fit, minmax(650px, 1fr));
gap: 1rem;
}}
.plot-container {{
background: var(--bg-primary);
border: 1px solid var(--border);
border-radius: 8px;
padding: 1rem;
}}
.plot-container h3 {{
font-size: 1rem;
margin-bottom: 0.5rem;
color: var(--text-secondary);
}}
/* Collapsible */
details {{ margin: 0.5rem 0; }}
summary {{
cursor: pointer;
font-weight: 600;
padding: 0.5rem;
background: var(--bg-primary);
border-radius: 6px;
border: 1px solid var(--border);
}}
summary:hover {{ background: rgba(99,102,241,0.1); }}
details > div {{ padding: 1rem; }}
.note {{ color: var(--text-secondary); font-size: 0.9rem; margin: 0.5rem 0; }}
/* Print styles */
@media print {{
body {{ background: white; color: black; }}
.section {{ border: 1px solid #ccc; page-break-inside: avoid; }}
}}
.two-col {{ display: grid; grid-template-columns: 1fr 1fr; gap: 1.5rem; }}
@media (max-width: 900px) {{ .two-col {{ grid-template-columns: 1fr; }} }}
</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">
<div class="logo">ATOMIZER</div>
<div>by Atomaste</div>
<div style="margin-top:0.3rem">FEA Optimization Platform</div>
</div>
</div>
<!-- Executive Summary -->
<div class="section">
<h2>📋 Executive Summary</h2>
<div class="exec-grid">
<div class="exec-card">
<div class="label">WFE 40° vs 20° (Tracking)</div>
<div class="value">{wfe_40_20:.2f} <small>nm</small></div>
<div class="target">{status_badge(wfe_40_20, targets['wfe_40_20'])}</div>
</div>
<div class="exec-card">
<div class="label">WFE 60° vs 20° (Tracking)</div>
<div class="value">{wfe_60_20:.2f} <small>nm</small></div>
<div class="target">{status_badge(wfe_60_20, targets['wfe_60_20'])}</div>
</div>
<div class="exec-card">
<div class="label">MFG 90° (J1J3 Filtered)</div>
<div class="value">{mfg_90:.2f} <small>nm</small></div>
<div class="target">{status_badge(mfg_90, targets['mfg_90'])}</div>
</div>
<div class="exec-card">
<div class="label">Weighted Sum (6·W40 + 5·W60 + 3·MFG)</div>
<div class="value" style="color: var(--accent)">{ws:.1f}</div>
<div class="target">Lower is better</div>
</div>
</div>
{'<p class="note">Annular aperture: inner radius = ' + f'{inner_radius:.1f} mm (ø{2*inner_radius:.1f} mm central hole)' + '</p>' if inner_radius else ''}
</div>
<!-- Per-Angle Summary -->
<div class="section">
<h2>📊 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 = J1J4 removed. Relative = vs 20° reference. Aberrations are absolute.</p>
</div>
<!-- Surface Plots -->
<div class="section">
<h2>🌊 Wavefront Error Surface Maps</h2>
<p class="note">3D residual surfaces after removing piston, tip, tilt, and defocus (J1J4). Interactive — drag to rotate.</p>
<div class="plot-grid">
<div class="plot-container">
<h3>40° vs 20° (Relative)</h3>
{surf_40}
</div>
<div class="plot-container">
<h3>60° vs 20° (Relative)</h3>
{surf_60}
</div>
</div>
<div class="plot-container" style="margin-top:1rem">
<h3>90° Manufacturing (Absolute)</h3>
{surf_90}
</div>
</div>
<!-- Trajectory Analysis -->
{'<div class="section"><h2>📈 Zernike Trajectory Analysis</h2>' +
'<p class="note">Mode-specific integrated RMS across the operating elevation range. ' +
'The linear model c<sub>j</sub>(θ) = a<sub>j</sub>·Δsinθ + b<sub>j</sub>·Δcosθ decomposes gravity into axial and lateral components.</p>' +
traj_metrics_html +
'<div class="two-col" style="margin-top:1rem">' +
'<div class="plot-container"><h3>Mode RMS vs Elevation Angle</h3>' + traj_plot_html + '</div>' +
'<div class="plot-container"><h3>Axial vs Lateral Sensitivity</h3>' + sens_plot_html + '</div>' +
'</div></div>' if traj_result else ''}
<!-- Manufacturing Analysis -->
<div class="section">
<h2>🏭 Manufacturing Analysis (90° Orientation)</h2>
<p class="note">
The mirror is manufactured (polished) at 90° orientation. The "Correction" column shows the
aberrations that must be polished out to achieve the 20° operational figure.
</p>
{mfg_html}
</div>
<!-- Design Parameters -->
{params_html}
<!-- Zernike Coefficient Details -->
<div class="section">
<h2>🔬 Zernike Coefficient Details</h2>
<details>
<summary>40° vs 20° — Relative Coefficients</summary>
<div>{bar_40}</div>
</details>
<details>
<summary>60° vs 20° — Relative Coefficients</summary>
<div>{bar_60}</div>
</details>
<details>
<summary>90° — Absolute Coefficients</summary>
<div>{bar_90}</div>
</details>
</div>
<!-- Methodology -->
<div class="section" style="opacity:0.85">
<h2>📝 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>J1J4 (Piston, Tip, Tilt, Defocus)</td></tr>
<tr><td><b>WFE Calculation</b></td><td>WFE = 2 × Surface Error (reflective)</td></tr>
<tr><td><b>Displacement Unit</b></td><td>{DISP_UNIT} → nm ({NM_SCALE:.0e}×)</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° (polishing/measurement orientation)</td></tr>
<tr><td><b>MFG Objective</b></td><td>90°20° relative, J1J3 filtered (optician workload)</td></tr>
<tr><td><b>Weighted Sum</b></td><td>6×WFE(4020) + 5×WFE(6020) + 3×MFG(90)</td></tr>
{'<tr><td><b>Trajectory R²</b></td><td>' + f'{traj_result["linear_fit_r2"]:.6f}' + '</td></tr>' if traj_result else ''}
</tbody>
</table>
</div>
<!-- Footer -->
<div style="text-align:center; padding:2rem; color:var(--text-secondary); font-size:0.8rem;">
Generated by <b>Atomizer</b> Optical Report Generator &nbsp;|&nbsp; {timestamp}<br>
© 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()