From 59a435f1193cb2a78e93d889f7f1f325aaec879e Mon Sep 17 00:00:00 2001 From: Anto01 Date: Tue, 23 Dec 2025 15:03:32 -0500 Subject: [PATCH] feat: Add debug script for lateral displacement analysis MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds tests/debug_lateral_discrepancy.py to investigate differences between Zernike OPD lateral displacement reporting and Simcenter post-processing. Key findings documented: - OPD reports sqrt(dx² + dy²) - combined XY magnitude - Simcenter shows individual components (dx or dy) - Both are correct, OPD magnitude is more meaningful for optics 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 --- tests/debug_lateral_discrepancy.py | 181 +++++++++++++++++++++++++++++ 1 file changed, 181 insertions(+) create mode 100644 tests/debug_lateral_discrepancy.py diff --git a/tests/debug_lateral_discrepancy.py b/tests/debug_lateral_discrepancy.py new file mode 100644 index 00000000..514c3979 --- /dev/null +++ b/tests/debug_lateral_discrepancy.py @@ -0,0 +1,181 @@ +""" +Debug script to investigate the lateral displacement discrepancy between +Zernike OPD output and Simcenter post-processing. + +User observation: +- Zernike OPD shows max lateral XY displacement: ~0.2238 µm +- Simcenter shows XX displacement at 20deg: 7.457e-05 mm = 0.0746 µm + +Hypothesis: The Zernike "max_lateral_disp_um" is sqrt(dx² + dy²), not just dx or dy. +""" + +import sys +import os +sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + +import numpy as np +from pathlib import Path + +# Find a mirror study OP2 to analyze +STUDIES_PATH = Path(r"c:\Users\antoi\Atomizer\studies\M1_Mirror") + +def find_latest_op2(): + """Find a recent OP2 file to analyze.""" + op2_files = list(STUDIES_PATH.rglob("*.op2")) + if not op2_files: + print("No OP2 files found!") + return None + # Get most recent + return max(op2_files, key=lambda p: p.stat().st_mtime) + +def analyze_displacements(op2_path: Path): + """Analyze displacement components for all subcases.""" + from pyNastran.op2.op2 import OP2 + from pyNastran.bdf.bdf import BDF + + print(f"\n{'='*70}") + print(f"Analyzing: {op2_path.name}") + print(f"Path: {op2_path}") + print(f"{'='*70}") + + # Find BDF + bdf_path = None + for ext in ['.dat', '.bdf']: + candidate = op2_path.with_suffix(ext) + if candidate.exists(): + bdf_path = candidate + break + + if not bdf_path: + print("No BDF found, searching parent...") + for f in op2_path.parent.iterdir(): + if f.suffix.lower() in ['.dat', '.bdf']: + bdf_path = f + break + + if not bdf_path: + print("ERROR: No geometry file found!") + return + + print(f"BDF: {bdf_path.name}") + + # Read data + print("\nLoading OP2...") + op2 = OP2() + op2.read_op2(str(op2_path)) + + print("Loading BDF...") + bdf = BDF() + bdf.read_bdf(str(bdf_path)) + + node_geo = {int(nid): node.get_position() for nid, node in bdf.nodes.items()} + + SUBCASE_MAP = {'1': 90, '2': 20, '3': 40, '4': 60} + + print(f"\nNode count in BDF: {len(node_geo)}") + print(f"\n{'='*70}") + print("DISPLACEMENT ANALYSIS BY SUBCASE") + print(f"{'='*70}") + + 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) + angle = SUBCASE_MAP.get(label, label) + + print(f"\n--- Subcase {label} ({angle} deg) ---") + print(f"Nodes with displacement: {len(node_ids)}") + + # Extract displacement components + dx = dmat[:, 0] # X displacement (mm) + dy = dmat[:, 1] # Y displacement (mm) + dz = dmat[:, 2] # Z displacement (mm) + + # Compute statistics + dx_um = dx * 1000.0 # Convert mm to µm + dy_um = dy * 1000.0 + dz_um = dz * 1000.0 + + # Lateral magnitude (XY combined) + lateral_um = np.sqrt(dx_um**2 + dy_um**2) + + print(f"\nComponent Statistics (in um):") + print(f" X displacement (dx):") + print(f" Min: {np.min(dx_um):+.4f} um") + print(f" Max: {np.max(dx_um):+.4f} um") + print(f" RMS: {np.sqrt(np.mean(dx_um**2)):.4f} um") + + print(f" Y displacement (dy):") + print(f" Min: {np.min(dy_um):+.4f} um") + print(f" Max: {np.max(dy_um):+.4f} um") + print(f" RMS: {np.sqrt(np.mean(dy_um**2)):.4f} um") + + print(f" Z displacement (dz):") + print(f" Min: {np.min(dz_um):+.4f} um") + print(f" Max: {np.max(dz_um):+.4f} um") + print(f" RMS: {np.sqrt(np.mean(dz_um**2)):.4f} um") + + print(f"\n Lateral Magnitude (sqrt(dx^2 + dy^2)):") + print(f" Max: {np.max(lateral_um):.4f} um <-- This is what Zernike OPD reports!") + print(f" RMS: {np.sqrt(np.mean(lateral_um**2)):.4f} um") + + # Find the node with max lateral displacement + max_lat_idx = np.argmax(lateral_um) + max_lat_nid = int(node_ids[max_lat_idx]) + max_dx = dx_um[max_lat_idx] + max_dy = dy_um[max_lat_idx] + + print(f"\n Node with max lateral displacement: Node {max_lat_nid}") + print(f" dx = {max_dx:+.4f} um") + print(f" dy = {max_dy:+.4f} um") + print(f" sqrt(dx^2+dy^2) = {lateral_um[max_lat_idx]:.4f} um") + + # Compare with just max(|dx|) + max_abs_dx = np.max(np.abs(dx_um)) + max_abs_dy = np.max(np.abs(dy_um)) + print(f"\n For comparison (what you see in Simcenter):") + print(f" max(|dx|) = {max_abs_dx:.4f} um") + print(f" max(|dy|) = {max_abs_dy:.4f} um") + + print(f"\n{'='*70}") + print("EXPLANATION OF DISCREPANCY") + print(f"{'='*70}") + print(""" +The Zernike OPD insight reports "max_lateral_disp_um" as: + lateral = sqrt(dx^2 + dy^2) -- combined XY magnitude at each node + +Simcenter's "Displacement - Nodal, X" shows: + Just the X component (dx) at each node + +These are different metrics: +- If dx_max = 0.0746 um and dy is significant, then: + lateral = sqrt(0.0746^2 + dy^2) > 0.0746 um + +To match Simcenter exactly, look at the individual dx/dy/dz stats above. +The "max_lateral_um" in Zernike OPD is the MAGNITUDE of the XY vector, +not the individual X or Y components. + +For a node where both dx and dy are non-zero: + dx = 0.0746 um, dy = 0.21 um + lateral = sqrt(0.0746^2 + 0.21^2) = sqrt(0.0056 + 0.044) = sqrt(0.0496) = 0.22 um +""") + +if __name__ == '__main__': + # Find and analyze an OP2 file + if len(sys.argv) > 1: + op2_path = Path(sys.argv[1]) + else: + op2_path = find_latest_op2() + + if op2_path and op2_path.exists(): + analyze_displacements(op2_path) + else: + print("Please provide an OP2 file path as argument, or place OP2 files in the studies directory.") + print(f"\nUsage: python {sys.argv[0]} ")