#!/usr/bin/env python3
"""Validate and export the ABACUS R8 harmonic tripod intake shell.

This mirrors the browser geometry:
  k(q,u) = 2/W(u) [1 + eps(u) cos(6 pi q)]
  z0(q)  = h_a [1 - cos(6 pi q)] / 2

The shell is the original analytic surface restricted to u >= z0(q)/H.
The arch rise is solved so the filtered effective intake equals the swept disk
area of one 140 mm fan.
"""
from __future__ import annotations

import json
import math
from pathlib import Path

import numpy as np
import trimesh

OUT = Path(__file__).resolve().parent

W_B = 275.0
W_T = 180.0
H = 285.0
EPS0 = 0.20
WALL = 4.0
MIN_DRAFT = 0.005
N = 480
NZ = 101
FAN_D = 140.0
FILTER_FREE = 0.80
TARGET_EFFECTIVE_RATIO = 1.0
DENSITY = 7.2


def smoothstep(u):
    return 3*u*u - 2*u*u*u


def width_at(u):
    m = -2 * MIN_DRAFT * H
    return W_B + m*u + (W_T-W_B-m)*smoothstep(u)


def eps_at(u):
    return EPS0 * (1 - 6*u*u + 4*u*u*u)


def signed_area(p):
    return 0.5*np.sum(p[:,0]*np.roll(p[:,1],-1) - p[:,1]*np.roll(p[:,0],-1))


def centroid(p):
    cross = p[:,0]*np.roll(p[:,1],-1) - np.roll(p[:,0],-1)*p[:,1]
    a = 0.5*np.sum(cross)
    if abs(a) < 1e-12:
        return p.mean(axis=0)
    cx = np.sum((p[:,0] + np.roll(p[:,0],-1))*cross)/(6*a)
    cy = np.sum((p[:,1] + np.roll(p[:,1],-1))*cross)/(6*a)
    return np.array([cx,cy])


def section(W, eps, n=N):
    qmid = (np.arange(n)+0.5)/n
    psi = 2*np.pi*qmid + eps/3*np.sin(6*np.pi*qmid)
    ds = np.pi*W/n
    increments = ds*np.column_stack([np.cos(psi),np.sin(psi)])
    cumulative = np.vstack([np.zeros(2), np.cumsum(increments,axis=0)])
    closure = cumulative[-1]
    q = np.arange(n)/n
    p = cumulative[:-1] - q[:,None]*closure
    p -= centroid(p)
    if signed_area(p)<0:
        p=p[::-1]
    return p


def offset_inward(p, wall):
    prev=np.roll(p,1,axis=0); nxt=np.roll(p,-1,axis=0)
    t=nxt-prev
    t/=np.linalg.norm(t,axis=1)[:,None]
    n=np.column_stack([-t[:,1],t[:,0]])
    return p+wall*n


def norm_arch(q):
    return 0.5*(1-np.cos(6*np.pi*q))


us=np.linspace(0,1,NZ)
outer=np.array([section(width_at(u),eps_at(u)) for u in us])
inner=np.array([offset_inward(s,WALL) for s in outer])


def interp(arr,u):
    f=np.clip(u,0,1)*(NZ-1)
    j0=np.floor(f).astype(int)
    j0=np.minimum(j0,NZ-2)
    j1=j0+1
    t=f-j0
    # one u per perimeter point
    idx=np.arange(N)
    return arr[j0,idx]*(1-t)[:,None]+arr[j1,idx]*t[:,None]


def curtain(rise):
    q=np.arange(N)/N
    z=rise*norm_arch(q)
    u=z/H
    rim2=interp(outer,u)
    base=np.column_stack([outer[0],np.zeros(N)])
    rim=np.column_stack([rim2,z])
    area=0.0
    for i in range(N):
        j=(i+1)%N
        for a,b,c in [(base[i],base[j],rim[i]),(base[j],rim[j],rim[i])]:
            area += 0.5*np.linalg.norm(np.cross(b-a,c-a))
    return area,base,rim

fan_area=math.pi*(FAN_D/2)**2
target_effective=TARGET_EFFECTIVE_RATIO*fan_area
target_geo=target_effective/FILTER_FREE
lo=0.0; hi=80.0
while curtain(hi)[0] < target_geo:
    hi*=1.25
for _ in range(60):
    mid=(lo+hi)/2
    if curtain(mid)[0] < target_geo: lo=mid
    else: hi=mid
rise=(lo+hi)/2
intake_area,base_boundary,rim_boundary=curtain(rise)

# Build trimmed-domain shell mesh.
verts=[]
q=np.arange(N)/N
u0=rise*norm_arch(q)/H
for jt in range(NZ):
    t=jt/(NZ-1)
    u=u0+(1-u0)*t
    xy=interp(outer,u)
    z=H*u
    verts.extend(np.column_stack([xy,z]))
outer_offset=0
inner_offset=len(verts)
for jt in range(NZ):
    t=jt/(NZ-1)
    u=u0+(1-u0)*t
    xy=interp(inner,u)
    z=H*u
    verts.extend(np.column_stack([xy,z]))
verts=np.asarray(verts,float)
faces=[]
for k in range(NZ-1):
    for i in range(N):
        j=(i+1)%N
        a=k*N+i; b=k*N+j; c=(k+1)*N+j; d=(k+1)*N+i
        faces.extend([(a,b,c),(a,c,d)])
for k in range(NZ-1):
    for i in range(N):
        j=(i+1)%N
        a=inner_offset+k*N+i; b=inner_offset+(k+1)*N+i
        c=inner_offset+(k+1)*N+j; d=inner_offset+k*N+j
        faces.extend([(a,b,c),(a,c,d)])
# variable bottom rim
for i in range(N):
    j=(i+1)%N
    o0=i;o1=j;i0=inner_offset+i;i1=inner_offset+j
    faces.extend([(o0,i1,o1),(o0,i0,i1)])
# top rim
ob=(NZ-1)*N; ib=inner_offset+(NZ-1)*N
for i in range(N):
    j=(i+1)%N
    o0=ob+i;o1=ob+j;i0=ib+i;i1=ib+j
    faces.extend([(o0,o1,i1),(o0,i1,i0)])
mesh=trimesh.Trimesh(vertices=verts,faces=np.asarray(faces),process=True)
mesh.update_faces(mesh.unique_faces())
mesh.update_faces(mesh.nondegenerate_faces())
mesh.remove_unreferenced_vertices()
trimesh.repair.fix_normals(mesh,multibody=False)
if mesh.volume<0: mesh.invert()
mesh.export(OUT/'ABACUS_R8_harmonic_tripod_shell.stl')
mesh.export(OUT/'ABACUS_R8_harmonic_tripod_shell.obj')

# Intake curtain mesh for visualization.
cv=np.vstack([base_boundary,rim_boundary])
cf=[]
for i in range(N):
    j=(i+1)%N
    cf.extend([(i,j,N+i),(j,N+j,N+i)])
curtain_mesh=trimesh.Trimesh(vertices=cv,faces=np.asarray(cf),process=False)

# GLB scene.
scene=trimesh.Scene()
sm=mesh.copy();sm.visual.face_colors=[42,39,36,235]
cm=curtain_mesh.copy();cm.visual.face_colors=[64,173,201,75]
scene.add_geometry(sm,node_name='cast_shell')
scene.add_geometry(cm,node_name='intake_area_curtain')
scene.export(OUT/'ABACUS_R8_harmonic_tripod_shell.glb')

# Preview.
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

fig=plt.figure(figsize=(11,9),dpi=180)
ax=fig.add_subplot(111,projection='3d')
tri=mesh.vertices[mesh.faces[::22]]
ax.add_collection3d(Poly3DCollection(tri,facecolor=(.16,.145,.13,.88),edgecolor=(.06,.055,.05,.12),linewidth=.1))
ctri=curtain_mesh.vertices[curtain_mesh.faces[::4]]
ax.add_collection3d(Poly3DCollection(ctri,facecolor=(.20,.70,.82,.20),edgecolor=(.20,.75,.88,.22),linewidth=.2))
# lower rim
ax.plot(np.r_[rim_boundary[:,0],rim_boundary[0,0]],np.r_[rim_boundary[:,1],rim_boundary[0,1]],np.r_[rim_boundary[:,2],rim_boundary[0,2]],color='#d8bd77',linewidth=1.3)
# desk plane
bounds=mesh.bounds
xx=np.linspace(bounds[0,0]-20,bounds[1,0]+20,2); yy=np.linspace(bounds[0,1]-20,bounds[1,1]+20,2)
X,Y=np.meshgrid(xx,yy);Z=np.zeros_like(X)
ax.plot_surface(X,Y,Z,color=(.07,.075,.08),alpha=.55,shade=False)
ax.set_xlim(bounds[0,0]-15,bounds[1,0]+15);ax.set_ylim(bounds[0,1]-15,bounds[1,1]+15);ax.set_zlim(0,H)
ax.set_box_aspect([bounds[1,0]-bounds[0,0],bounds[1,1]-bounds[0,1],H])
ax.view_init(elev=18,azim=-52)
ax.set_axis_off()
ax.set_title('ABACUS R8 — analytic tripod intake boundary',pad=18)
fig.tight_layout()
fig.savefig(OUT/'ABACUS_R8_harmonic_tripod_preview.png',transparent=True)
plt.close(fig)

# Lower-edge profile plot.
fig,ax=plt.subplots(figsize=(11,4.5),dpi=180)
ax.plot(q,rise*norm_arch(q),linewidth=2.2,color='#b68f37')
for k in range(3):
    ax.axvline(k/3,color='black',alpha=.15,linewidth=.8)
    ax.axvline(k/3+1/6,color='#4fa9be',alpha=.18,linewidth=.8)
ax.set_xlabel('normalized bottom arclength q')
ax.set_ylabel('rim height above desk (mm)')
ax.set_title('Single-harmonic lower edge: three feet and three intake arches')
ax.grid(True,linewidth=.25,alpha=.4)
fig.tight_layout();fig.savefig(OUT/'ABACUS_R8_harmonic_tripod_rim_law.png');plt.close(fig)

metrics={
    'bottom_mean_width_mm':W_B,
    'height_mm':H,
    'fan_diameter_mm':FAN_D,
    'fan_disk_area_mm2':fan_area,
    'filter_free_area':FILTER_FREE,
    'target_geometric_intake_area_mm2':target_geo,
    'solved_arch_rise_mm':rise,
    'actual_intake_curtain_area_mm2':intake_area,
    'effective_intake_area_mm2':intake_area*FILTER_FREE,
    'effective_to_fan_ratio':intake_area*FILTER_FREE/fan_area,
    'per_arch_geometric_area_mm2':intake_area/3,
    'shell_watertight':bool(mesh.is_watertight),
    'shell_volume_cm3':float(mesh.volume/1000),
    'shell_mass_kg':float(mesh.volume/1000*DENSITY/1000),
    'euler_number':int(mesh.euler_number),
    'vertices':int(len(mesh.vertices)),
    'faces':int(len(mesh.faces)),
}
(OUT/'TRIPOD_GEOMETRY_METRICS.json').write_text(json.dumps(metrics,indent=2))
report=f'''# Harmonic tripod intake validation\n\n- Solved arch rise: **{rise:.3f} mm**\n- Total geometric intake curtain: **{intake_area/100:.2f} cm²**\n- Per arch: **{intake_area/300:.2f} cm²**\n- Filtered effective area / fan disk: **{metrics['effective_to_fan_ratio']:.6f}:1**\n- Shell watertight: **{mesh.is_watertight}**\n- Shell material volume: **{metrics['shell_volume_cm3']:.2f} cm³**\n- Estimated gray-iron mass: **{metrics['shell_mass_kg']:.2f} kg**\n\nThe lower edge is `z0(q)=h_a[1-cos(6πq)]/2`.  The solver evaluates the actual sloped curtain between the bottom footprint and the tapered analytic rim, rather than relying only on the zero-taper Barbier estimate.\n'''
(OUT/'TRIPOD_GEOMETRY_REPORT.md').write_text(report)
print(json.dumps(metrics,indent=2))
