"""
Script to investigate unit conversion issues in AtomizerField
"""
import json
import h5py
import numpy as np

print("="*60)
print("UNIT VALIDATION CHECK")
print("="*60)

# Load JSON metadata
with open('test_case_beam/neural_field_data.json', 'r') as f:
    data = json.load(f)

# Check units
print("\n1. UNITS FROM METADATA:")
print("-"*60)
units = data['metadata']['units']
for key, value in units.items():
    print(f"  {key}: {value}")

# Check loads
print("\n2. APPLIED LOADS:")
print("-"*60)
loads = data['loads']
print(f"  Available load types: {list(loads.keys())}")

if 'point_forces' in loads:
    point_forces = loads['point_forces']
    print(f"  Point forces: {len(point_forces)} applied")

    if point_forces:
        print("\n  Sample forces (first 3):")
        for i in range(min(3, len(point_forces))):
            force = point_forces[i]
            print(f"    Force {i+1}:")
            print(f"      Node: {force['node']}")
            print(f"      Magnitude: {force['magnitude']:.2e} N")
            print(f"      Direction: {force['direction']}")

# Load HDF5 data
print("\n3. REACTION FORCES FROM HDF5:")
print("-"*60)
with h5py.File('test_case_beam/neural_field_data.h5', 'r') as f:
    reactions = f['results/reactions'][:]

    # Get statistics
    non_zero_reactions = reactions[np.abs(reactions) > 1e-10]
    print(f"  Shape: {reactions.shape}")
    print(f"  Min: {np.min(reactions):.2e}")
    print(f"  Max: {np.max(reactions):.2e}")
    print(f"  Mean (non-zero): {np.mean(np.abs(non_zero_reactions)):.2e}")
    print(f"  Max absolute: {np.max(np.abs(reactions)):.2e}")

    # Check displacement for comparison
    displacement = f['results/displacement'][:]
    max_disp = np.max(np.abs(displacement[:, :3]))  # Translation only
    print(f"\n  Max displacement: {max_disp:.6f} mm")

# Check stress
print("\n4. STRESS VALUES FROM HDF5:")
print("-"*60)
with h5py.File('test_case_beam/neural_field_data.h5', 'r') as f:
    stress = f['results/stress/cquad4_stress/data'][:]

    # Von Mises stress is in column 7 (0-indexed)
    von_mises = stress[:, 7]

    print(f"  Shape: {stress.shape}")
    print(f"  Von Mises min: {np.min(von_mises):.2e} MPa")
    print(f"  Von Mises max: {np.max(von_mises):.2e} MPa")
    print(f"  Von Mises mean: {np.mean(von_mises):.2e} MPa")

# Check if values make sense
print("\n5. SANITY CHECK:")
print("-"*60)
print(f"  Units claimed: force=N, stress=MPa, length=mm")
print(f"  Max reaction force: {np.max(np.abs(reactions)):.2e} N")
print(f"  Max von Mises: {np.max(von_mises):.2e} MPa")
print(f"  Max displacement: {max_disp:.6f} mm")

# Typical beam: if force is 1000 N, stress should be ~10-100 MPa
# If reaction is 152 million N, that's 152,000 kN - VERY high!
max_reaction = np.max(np.abs(reactions))
max_stress_val = np.max(von_mises)

print(f"\n  If force unit is actually kN instead of N:")
print(f"    Max reaction: {max_reaction/1000:.2e} kN")
print(f"  If stress unit is actually Pa instead of MPa:")
print(f"    Max stress: {max_stress_val/1e6:.2e} MPa")

print("\n6. HYPOTHESIS:")
print("-"*60)
if max_reaction > 1e6:
    print("  [!] Reaction forces seem TOO LARGE (>1 MN)")
    print("  [!] Possible issue: pyNastran returns forces in different units")
    print("  [!] Check: Nastran may export in base units (N) while expecting kN")

if max_stress_val > 1e6:
    print("  [!] Stresses seem TOO LARGE (>1000 MPa)")
    print("  [!] Possible issue: pyNastran returns stress in Pa, not MPa")

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