"""
Check actual units from pyNastran OP2 reader
"""
from pyNastran.op2.op2 import OP2
import numpy as np

print("="*60)
print("CHECKING OP2 FILE UNITS")
print("="*60)

# Load OP2 file
op2_file = 'test_case_beam/output/model.op2'
print(f"\nLoading: {op2_file}")

op2 = OP2()
op2.read_op2(op2_file)

print("\n1. DISPLACEMENT DATA:")
print("-"*60)
if hasattr(op2, 'displacements') and op2.displacements:
    for subcase_id, disp in op2.displacements.items():
        print(f"  Subcase {subcase_id}:")
        print(f"    Data shape: {disp.data.shape}")
        print(f"    Max displacement (all DOF): {np.max(np.abs(disp.data)):.6f}")
        print(f"    Translation max (DOF 1-3): {np.max(np.abs(disp.data[0, :, :3])):.6f}")

        # Check if pyNastran has unit info
        if hasattr(disp, 'units'):
            print(f"    Units: {disp.units}")
        else:
            print(f"    Units: Not specified by pyNastran")

print("\n2. STRESS DATA:")
print("-"*60)
# Try new API first
stress_dict = None
if hasattr(op2, 'op2_results') and hasattr(op2.op2_results, 'stress'):
    if hasattr(op2.op2_results.stress, 'cquad4_stress'):
        stress_dict = op2.op2_results.stress.cquad4_stress
elif hasattr(op2, 'cquad4_stress'):
    try:
        stress_dict = op2.cquad4_stress
    except:
        stress_dict = None

if stress_dict:
    for subcase_id, stress in stress_dict.items():
        print(f"  Subcase {subcase_id}:")
        print(f"    Data shape: {stress.data.shape}")

        # Get von Mises stress (last column)
        von_mises = stress.data[0, :, -1]
        print(f"    Von Mises min: {np.min(von_mises):.2e}")
        print(f"    Von Mises max: {np.max(von_mises):.2e}")
        print(f"    Von Mises mean: {np.mean(von_mises):.2e}")

        # Check if pyNastran has unit info
        if hasattr(stress, 'units'):
            print(f"    Units: {stress.units}")
        else:
            print(f"    Units: Not specified by pyNastran")

        # Check data type names
        if hasattr(stress, 'data_names'):
            print(f"    Data names: {stress.data_names}")
else:
    print("  No CQUAD4 stress data found")

print("\n3. REACTION FORCES:")
print("-"*60)
if hasattr(op2, 'grid_point_forces') and op2.grid_point_forces:
    for subcase_id, forces in op2.grid_point_forces.items():
        print(f"  Subcase {subcase_id}:")
        print(f"    Data shape: {forces.data.shape}")
        print(f"    Max force: {np.max(np.abs(forces.data)):.2e}")

        if hasattr(forces, 'units'):
            print(f"    Units: {forces.units}")
        else:
            print(f"    Units: Not specified by pyNastran")

print("\n4. BDF FILE UNITS:")
print("-"*60)
# Try to read BDF to check units
from pyNastran.bdf.bdf import BDF
bdf_file = 'test_case_beam/input/model.bdf'
print(f"Loading: {bdf_file}")

bdf = BDF()
bdf.read_bdf(bdf_file)

# Check for PARAM cards that define units
if hasattr(bdf, 'params') and bdf.params:
    print("  PARAM cards found:")
    for param_name, param in bdf.params.items():
        if 'UNIT' in param_name.upper() or 'WTMASS' in param_name.upper():
            print(f"    {param_name}: {param}")

# Check material properties (can infer units from magnitude)
if hasattr(bdf, 'materials') and bdf.materials:
    print("\n  Material properties (first material):")
    mat = list(bdf.materials.values())[0]
    print(f"    Type: {mat.type}")
    if hasattr(mat, 'e') and mat.e:
        print(f"    Young's modulus E: {mat.e:.2e}")
        print(f"      If E~200,000 then units are MPa (steel)")
        print(f"      If E~200,000,000 then units are Pa (steel)")
        print(f"      If E~29,000,000 then units are psi (steel)")

print("\n5. NASTRAN UNIT SYSTEM ANALYSIS:")
print("-"*60)
print("  Common Nastran unit systems:")
print("    - SI: m, kg, N, Pa, s")
print("    - mm-tonne-N: mm, tonne, N, MPa, s")
print("    - mm-kg-N: mm, kg, N, MPa, s")
print("    - in-lb-s: in, lb, lbf, psi, s")
print("")
print("  Our metadata claims: mm, kg, N, MPa")
print("  But pyNastran might return stress in Pa (base SI)")
print("  This would explain the 1000× error!")

print("\n" + "="*60)