import struct, numpy as np

path = '/app/projects/esp32-box/original/Window 20x60mm.stl'

with open(path, 'rb') as f:
    f.read(80)
    n = struct.unpack('<I', f.read(4))[0]
    verts = []
    for i in range(n):
        f.read(12)
        tri = []
        for j in range(3):
            x, y, z = struct.unpack('<fff', f.read(12))
            tri.append([x, y, z])
        f.read(2)
        verts.append(tri)

verts = np.array(verts)

# Show precise outline at key Z levels
z_levels = [0.0, 1.6, 3.2, 4.8]
for z_lev in z_levels:
    xy_points = []
    for tri in verts:
        zs = tri[:,2]
        if abs(zs.min() - z_lev) < 0.001 or abs(zs.max() - z_lev) < 0.001:
            for i in range(3):
                if abs(tri[i,2] - z_lev) < 0.001:
                    xy_points.append([tri[i,0], tri[i,1]])
    if xy_points:
        xy = np.array(xy_points)
        x_c = xy[:,0].copy()
        y_c = xy[:,1].copy()
        print('Z=%.1f:' % z_lev)
        print('  X: %.2f..%.2f (%.2f)' % (x_c.min(), x_c.max(), x_c.max()-x_c.min()))
        print('  Y: %.2f..%.2f (%.2f)' % (y_c.min(), y_c.max(), y_c.max()-y_c.min()))
        print('  X sorted:', np.sort(np.unique(np.round(x_c, 1))))
        print('  Y sorted:', np.sort(np.unique(np.round(y_c, 1))))

# Analyze the inner cutout features at Z=1.6 level (the keystone/DC cutouts)
print()
print('=== Interior cutout features at Z=1.6 ===')
for tri in verts:
    zs = tri[:,2]
    if zs.max() <= 1.61 and zs.min() >= 0.0:
        pass  # These are the flat base at Z=1.6 (no cutout)

# Let me look at what happens in the interior (not at edges)
# Check vertices that are NOT on the outer perimeter
print()
print('=== Interior vertices (not on outer XY perimeter) ===')
all_pts = verts.reshape(-1, 3)
x_min, x_max = all_pts[:,0].min(), all_pts[:,0].max()
y_min, y_max = all_pts[:,1].min(), all_pts[:,1].max()

# Look at center region
center_pts = all_pts[(all_pts[:,0] > x_min+2) & (all_pts[:,0] < x_max-2) & 
                     (all_pts[:,1] > y_min+2) & (all_pts[:,1] < y_max-2)]
print('Interior points:', len(center_pts))
z_u = np.unique(np.round(center_pts[:,2], 2))
for z in z_u:
    count = np.sum(np.round(center_pts[:,2], 2) == z)
    print('  Z=%.2f: %d pts' % (z, count))
    
# Look at the X centerline profile
print()
print('=== X centerline (Y=%s) at various Z ===')
y_mid = all_pts[:,1].mean()
print('  Y center:', y_mid)
for z_lev in [0, 0.5, 1.0, 1.6, 2.0, 2.5, 3.0, 3.2, 3.5, 4.0, 4.5, 4.8]:
    mask = np.abs(all_pts[:,2] - z_lev) < 0.15
    if mask.any():
        xs = all_pts[mask, 0]
        ys = all_pts[mask, 1]
        print('  Z=%.1f: X=[%.2f, %.2f] Y=[%.2f, %.2f]' % (z_lev, xs.min(), xs.max(), ys.min(), ys.max()))

# Let's try to find the curves - look at X profile across Y at mid Z
print()
print('=== X profile at mid-Y across the window ===')
for z_lev in [0, 1.6, 3.2, 4.8]:
    print('Z=%.1f:' % z_lev, end=' ')
    for y_pos in np.linspace(y_min+0.5, y_max-0.5, 8):
        mask = (np.abs(all_pts[:,2] - z_lev) < 0.15) & (np.abs(all_pts[:,1] - y_pos) < 0.15)
        if mask.any():
            xs = all_pts[mask, 0]
            if 'x_min_p' not in dir():
                pass
            # Find gaps in X (cutouts)
            x_sorted = np.unique(np.round(xs, 1))
            gaps = []
            for i in range(len(x_sorted)-1):
                if x_sorted[i+1] - x_sorted[i] > 2:
                    gaps.append((x_sorted[i], x_sorted[i+1]))
            if gaps:
                print('\n    Y=%.1f gaps:' % y_pos, gaps)