#!/usr/bin/env python3
"""MuJoCo coffee-hopper discharge: model, 5 s simulation, and one report figure."""
from pathlib import Path
import os
import math
import numpy as np
import mujoco
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from matplotlib.lines import Line2D

# ----------------------------- declared model -----------------------------
SEED = 19
N_BEANS = 50
DT = 0.0005
DURATION = 5.0
GATE_TIME = 1.0
BEAN_MASS = 0.00018          # kg, each complete bean
BEAN_SIZE = (0.0100, 0.0072, 0.0052)  # overall nominal bounding size, m
THROAT_DIAMETER = 0.018
WALL_FRICTION = (0.28, 0.005, 0.0001)
BEAN_FRICTION = (0.45, 0.005, 0.0001)
R_OUT = 0.060
Z_CONE_BOTTOM, Z_CONE_TOP = 0.075, 0.200
Z_THROAT_BOTTOM = 0.057
Z_GATE = 0.0545
CROSS_Z = Z_THROAT_BOTTOM
SNAPSHOT_TIMES = (0.0, 1.0, 2.0, 5.0)
DIAG_TIME = 1.35
rng = np.random.default_rng(SEED)


def initial_beans():
    """Deterministic, non-overlapping bounding-sphere packing in the cone."""
    pts = []
    while len(pts) < N_BEANS:
        # Sample cone volume, biased naturally toward its broad upper end.
        z = rng.uniform(Z_CONE_BOTTOM + 0.010, Z_CONE_TOP - 0.014)
        rmax = 0.009 + (z - Z_CONE_BOTTOM) / (Z_CONE_TOP - Z_CONE_BOTTOM) * (R_OUT - 0.009)
        r = (rmax - 0.0065) * math.sqrt(rng.random())
        a = rng.uniform(0, 2 * np.pi)
        p = np.array([r * np.cos(a), r * np.sin(a), z])
        if all(np.linalg.norm(p - q) > 0.0108 for q in pts):
            pts.append(p)
    return np.array(pts)


def ring_mesh(r0, r1, z0, z1, thickness=0.0025, n=64):
    """Closed hollow tapered tube, with inner and outer walls and annular rims."""
    verts = []
    for r, z in [(r0, z0), (r1, z1), (r0 + thickness, z0), (r1 + thickness, z1)]:
        verts.extend([(r*np.cos(2*np.pi*i/n), r*np.sin(2*np.pi*i/n), z) for i in range(n)])
    faces = []
    for i in range(n):
        j = (i + 1) % n
        # inner, outer, lower rim, upper rim (triangles)
        faces += [(i,j,n+j),(i,n+j,n+i)]
        faces += [(2*n+i,3*n+j,2*n+j),(2*n+i,3*n+i,3*n+j)]
        faces += [(i,2*n+j,j),(i,2*n+i,2*n+j)]
        faces += [(n+i,n+j,3*n+j),(n+i,3*n+j,3*n+i)]
    # Inner normals point into the flow channel and outer normals point out.
    return np.array(verts), np.array(faces, dtype=int)[:, ::-1]


def wall_wedges(r0, r1, z0, z1, prefix, thickness=0.0025, n=48):
    """A hollow wall as convex collidable wedges (not one convex hull)."""
    assets, geoms = [], []
    for k in range(n):
        a, b = 2*np.pi*k/n, 2*np.pi*(k+1)/n
        verts=[]
        for r,z,ang in [(r0,z0,a),(r0,z0,b),(r1,z1,b),(r1,z1,a),
                        (r0+thickness,z0,a),(r0+thickness,z0,b),
                        (r1+thickness,z1,b),(r1+thickness,z1,a)]:
            verts.append((r*np.cos(ang),r*np.sin(ang),z))
        faces=[(0,2,1),(0,3,2), (4,5,6),(4,6,7),
               (0,5,4),(0,1,5), (3,6,2),(3,7,6),
               (0,3,7),(0,7,4), (1,5,6),(1,6,2)]
        name=f"{prefix}_{k:02d}"
        assets.append(mesh_xml(name,np.array(verts),np.array(faces)))
        geoms.append(f'<geom name="{name}" type="mesh" mesh="{name}" material="steel" '
                     f'friction="{WALL_FRICTION[0]} {WALL_FRICTION[1]} {WALL_FRICTION[2]}"/>')
    return "\n    ".join(assets), "\n      ".join(geoms)


def mesh_xml(name, verts, faces):
    v = " ".join(f"{x:.7f} {y:.7f} {z:.7f}" for x,y,z in verts)
    f = " ".join(f"{a} {b} {c}" for a,b,c in faces)
    return f'<mesh name="{name}" vertex="{v}" face="{f}"/>'


def make_xml():
    cone_assets, cone_geoms = wall_wedges(0.009, R_OUT, Z_CONE_BOTTOM, Z_CONE_TOP, "hopper_wall", n=48)
    tube_assets, tube_geoms = wall_wedges(0.009, 0.009, Z_THROAT_BOTTOM, Z_CONE_BOTTOM, "throat_wall", n=32)
    positions = initial_beans()
    beans = []
    for i, p in enumerate(positions):
        # A thin, broad ellipsoid forms the flatter face; the offset fuller
        # ellipsoid forms the domed face. Both belong to the SAME free body.
        beans.append(f'''
    <body name="bean_{i+1:02d}" pos="{p[0]:.7f} {p[1]:.7f} {p[2]:.7f}">
      <freejoint name="free_{i+1:02d}"/>
      <geom name="bean_{i+1:02d}_flat" type="ellipsoid" size="0.0050 0.0036 0.00125"
            pos="0 0 -0.00065" quat="1 0 0 0" mass="0.00009" rgba="0.25 0.085 0.025 1"/>
      <geom name="bean_{i+1:02d}_dome" type="ellipsoid" size="0.00475 0.0034 0.00215"
            pos="0 0 0.00045" quat="1 0 0 0" mass="0.00009" rgba="0.40 0.16 0.055 1"/>
    </body>''')
    xml = f'''<mujoco model="50_bean_grinder_hopper">
  <compiler angle="radian" coordinate="local" inertiafromgeom="true"/>
  <option timestep="{DT}" gravity="0 0 -9.81" integrator="implicitfast"
          solver="Newton" iterations="100" ls_iterations="50" tolerance="1e-10"
          cone="elliptic" jacobian="dense" noslip_iterations="5"/>
  <size memory="64M"/>
  <visual><headlight ambient="0.35 0.35 0.35" diffuse="0.65 0.65 0.65" specular="0.2 0.2 0.2"/></visual>
  <default>
    <geom condim="4" friction="{BEAN_FRICTION[0]} {BEAN_FRICTION[1]} {BEAN_FRICTION[2]}"
          solref="0.003 1" solimp="0.92 0.99 0.001" margin="0" gap="0"/>
  </default>
  <asset>
    {cone_assets}
    {tube_assets}
    <material name="steel" rgba="0.55 0.31 0.15 1"/>
  </asset>
  <worldbody>
    <body name="hopper" pos="0 0 0">
      {cone_geoms}
      {tube_geoms}
    </body>
    <body name="gate" pos="0 0 {Z_GATE}">
      <joint name="gate_slide" type="slide" axis="1 0 0" range="0 0.055" damping="0.08"/>
      <geom name="gate_plate" type="box" size="0.023 0.020 0.0015" mass="0.025"
            rgba="0.80 0.48 0.12 1" friction="0.30 0.005 0.0001"/>
    </body>
    <body name="tray" pos="0 0 0">
      <geom name="tray_floor" type="box" size="0.075 0.060 0.002" pos="0 0 -0.002"
            rgba="0.30 0.34 0.36 1" friction="0.60 0.008 0.0001"/>
      <geom name="tray_front" type="box" size="0.075 0.002 0.005" pos="0 -0.060 0.003" rgba="0.30 0.34 0.36 1"/>
      <geom name="tray_back" type="box" size="0.075 0.002 0.005" pos="0 0.060 0.003" rgba="0.30 0.34 0.36 1"/>
      <geom name="tray_left" type="box" size="0.002 0.060 0.005" pos="-0.075 0 0.003" rgba="0.30 0.34 0.36 1"/>
      <geom name="tray_right" type="box" size="0.002 0.060 0.005" pos="0.075 0 0.003" rgba="0.30 0.34 0.36 1"/>
    </body>
    {''.join(beans)}
  </worldbody>
  <contact>
    <exclude body1="hopper" body2="gate"/>
    <exclude body1="tray" body2="gate"/>
  </contact>
  <actuator>
    <position name="gate_open" joint="gate_slide" kp="3" ctrlrange="0 0.055" forcerange="-4 4"/>
  </actuator>
</mujoco>'''
    return xml


# ------------------------------- simulation -------------------------------
xml_path = Path(__file__).resolve().with_name("model.xml")
xml_path.write_text(make_xml())
model = mujoco.MjModel.from_xml_path(str(xml_path))
data = mujoco.MjData(model)
mujoco.mj_forward(model, data)
qadr = np.array([model.jnt_qposadr[model.joint(name=f"free_{i:02d}").id] for i in range(1, N_BEANS+1)])
steps = int(round(DURATION / DT))
snap_steps = {int(round(t / DT)): t for t in SNAPSHOT_TIMES}
diag_step = int(round(DIAG_TIME / DT))
states, contacts_diag = {}, None
crossed = np.zeros(N_BEANS, dtype=bool)
cross_step = np.full(N_BEANS, -1, dtype=int)
count_history = np.zeros(steps + 1, dtype=int)
max_penetration = 0.0


def capture(step):
    states[step] = {
        "qpos": data.qpos.copy(), "time": data.time,
        "gate": data.qpos[model.jnt_qposadr[model.joint("gate_slide").id]],
        "count": int(crossed.sum())
    }


def contact_snapshot():
    out = []
    for k in range(data.ncon):
        c = data.contact[k]
        force6 = np.zeros(6)
        mujoco.mj_contactForce(model, data, k, force6)
        out.append((c.pos.copy(), c.frame.reshape(3,3).copy(), force6.copy(),
                    model.geom(c.geom1).name, model.geom(c.geom2).name))
    return out


capture(0)
count_history[0] = 0
for step in range(1, steps + 1):
    data.ctrl[0] = 0.0 if data.time < GATE_TIME else 0.055
    mujoco.mj_step(model, data)
    centers = np.array([data.qpos[a:a+3] for a in qadr])
    newly = (~crossed) & (centers[:, 2] < CROSS_Z)
    cross_step[newly] = step
    crossed[newly] = True
    count_history[step] = int(crossed.sum())
    for k in range(data.ncon):
        c = data.contact[k]
        names = (model.geom(c.geom1).name, model.geom(c.geom2).name)
        if any(n.startswith("hopper_wall_") or n.startswith("throat_wall_") for n in names):
            max_penetration = max(max_penetration, max(0.0, -float(c.dist)))
    if step in snap_steps:
        capture(step)
    if step == diag_step:
        capture(step)
        contacts_diag = contact_snapshot()

# Ensure 5 s state is captured even if a timestep convention changes.
if steps not in states:
    capture(steps)
centers_final = np.array([data.qpos[a:a+3] for a in qadr])
left_ids = [i+1 for i,p in enumerate(centers_final)
            if (not np.all(np.isfinite(p))) or abs(p[0]) > 0.16 or abs(p[1]) > 0.14 or p[2] < -0.05 or p[2] > 0.30]

print(f"Simulated {data.time:.3f} s; {N_BEANS} beans; timestep {DT*1000:.2f} ms")
print(f"Maximum penetration into hopper/throat walls: {max_penetration*1000:.4f} mm")
for t in SNAPSHOT_TIMES:
    st = states[int(round(t/DT))]
    print(f"Through throat at t={t:.0f} s: {st['count']} / {N_BEANS}")
print("Beans that left the model: " + (", ".join(f"bean_{i:02d}" for i in left_ids) if left_ids else "none"))
print(f"MJCF written to {xml_path}")

# ----------------------------- drawing helpers ----------------------------
BG = "#f5f1e9"
INK = "#243138"
MUTED = "#6e7776"
COPPER = "#b77b45"
GOLD = "#e5a93d"
TEAL = "#167d83"


def quat_matrix(q):
    w,x,y,z = q
    return np.array([[1-2*(y*y+z*z),2*(x*y-z*w),2*(x*z+y*w)],
                     [2*(x*y+z*w),1-2*(x*x+z*z),2*(y*z-x*w)],
                     [2*(x*z-y*w),2*(y*z+x*w),1-2*(x*x+y*y)]])


def ellipsoid_faces(center, rot, radii, offset, base, all_faces, all_colors, nu=10, nv=6):
    # Local z separates flatter and domed lobes. Faces receive simple studio light.
    lat = np.linspace(-np.pi/2, np.pi/2, nv+1)
    lon = np.linspace(0, 2*np.pi, nu+1)
    local = []
    for b in lat:
        for a in lon:
            local.append(np.array([radii[0]*np.cos(b)*np.cos(a),
                                   radii[1]*np.cos(b)*np.sin(a), radii[2]*np.sin(b)]))
    local = np.array(local) + np.array(offset)
    vv = local @ rot.T + center
    light = np.array([-0.4,-0.5,0.75]); light /= np.linalg.norm(light)
    rgb = np.array(matplotlib.colors.to_rgb(base))
    for j in range(nv):
        for i in range(nu):
            ids = [j*(nu+1)+i,j*(nu+1)+i+1,(j+1)*(nu+1)+i+1,(j+1)*(nu+1)+i]
            face = vv[ids]
            n = np.cross(face[1]-face[0],face[-1]-face[0])
            nn = np.linalg.norm(n)
            shade = 0.68 if nn == 0 else 0.60 + 0.40*abs(np.dot(n/nn,light))
            all_faces.append(face)
            all_colors.append(np.clip(rgb*shade,0,1))


def line3(ax, pts, **kw):
    pts = np.asarray(pts)
    ax.plot(pts[:,0],pts[:,1],pts[:,2],**kw)


def ring3(ax, r, z, **kw):
    a = np.linspace(0,2*np.pi,100)
    line3(ax, np.c_[r*np.cos(a),r*np.sin(a),np.full_like(a,z)],**kw)


def draw_hopper(ax, cutaway=True, throat=True):
    # Back 220-degree shell is translucent; omitted front sector reveals beans.
    front = np.deg2rad(-58)
    angles = np.linspace(front+np.deg2rad(70),front+2*np.pi-np.deg2rad(70),45)
    verts, faces = [], []
    for z,r in [(Z_CONE_BOTTOM,0.009),(Z_CONE_TOP,R_OUT)]:
        verts.extend([(r*np.cos(a),r*np.sin(a),z) for a in angles])
    n = len(angles)
    for i in range(n-1):
        faces.append([verts[i],verts[i+1],verts[n+i+1],verts[n+i]])
    ax.add_collection3d(Poly3DCollection(faces,facecolors=COPPER,alpha=.16,edgecolors="none"))
    for z,r in [(Z_CONE_BOTTOM,.009),(Z_CONE_TOP,R_OUT)]:
        ring3(ax,r,z,color=COPPER,lw=1.25,alpha=.85)
    for a in np.linspace(0,2*np.pi,12,endpoint=False):
        line3(ax,[[.009*np.cos(a),.009*np.sin(a),Z_CONE_BOTTOM],
                  [R_OUT*np.cos(a),R_OUT*np.sin(a),Z_CONE_TOP]],color=COPPER,lw=.45,alpha=.28)
    if throat:
        a = np.linspace(front+np.pi/3,front+5*np.pi/3,30)
        v = []
        for z in [Z_THROAT_BOTTOM,Z_CONE_BOTTOM]:
            v.extend([[.009*np.cos(t),.009*np.sin(t),z] for t in a])
        nn=len(a)
        ff=[[v[i],v[i+1],v[nn+i+1],v[nn+i]] for i in range(nn-1)]
        ax.add_collection3d(Poly3DCollection(ff,facecolors=COPPER,alpha=.20,edgecolors="none"))
        ring3(ax,.009,Z_THROAT_BOTTOM,color=COPPER,lw=1.0)


def draw_tray(ax):
    x,y=.075,.060
    floor=[[-x,-y,0],[x,-y,0],[x,y,0],[-x,y,0]]
    ax.add_collection3d(Poly3DCollection([floor],facecolors="#b9c3c2",alpha=.30,edgecolors="none"))
    for z in [0,.008]:
        line3(ax,[[-x,-y,z],[x,-y,z],[x,y,z],[-x,y,z],[-x,-y,z]],color="#75817f",lw=.8,alpha=.75)
    for xx,yy in [(-x,-y),(x,-y),(x,y),(-x,y)]:
        line3(ax,[[xx,yy,0],[xx,yy,.008]],color="#75817f",lw=.7)


def draw_gate(ax, x):
    # Gate plate follows the simulated joint position.
    xx,yy,h=.023,.020,.0015
    v=np.array([[x+a,y+b,Z_GATE+c] for a,b,c in []]) if False else None
    corners=np.array([[x-xx,-yy,Z_GATE-h],[x+xx,-yy,Z_GATE-h],[x+xx,yy,Z_GATE-h],[x-xx,yy,Z_GATE-h]])
    ax.add_collection3d(Poly3DCollection([corners],facecolors=GOLD,alpha=.80,edgecolors="#956322",linewidths=.6))


def draw_beans(ax, qpos, highlight=None):
    ff,cc=[],[]
    for i,a in enumerate(qadr):
        p=qpos[a:a+3]
        rot=quat_matrix(qpos[a+3:a+7])
        ellipsoid_faces(p,rot,(.0050,.0036,.00125),(0,0,-.00065),"#6e301b",ff,cc)
        ellipsoid_faces(p,rot,(.00475,.0034,.00215),(0,0,.00045),"#a55428",ff,cc)
    coll=Poly3DCollection(ff,facecolors=cc,edgecolors="none",alpha=1.0)
    ax.add_collection3d(coll)


def setup_ax(ax, zoom=False):
    ax.set_facecolor(BG)
    ax.view_init(elev=19,azim=-58)
    ax.set_proj_type('ortho')
    if zoom:
        ax.set_xlim(-.041,.041); ax.set_ylim(-.036,.036); ax.set_zlim(.020,.142)
        ax.set_box_aspect((1,1,1.7))
    else:
        ax.set_xlim(-.079,.079); ax.set_ylim(-.069,.069); ax.set_zlim(-.008,.208)
        ax.set_box_aspect((1,1,1.45))
    ax.set_axis_off()


def draw_state(ax, step, title, subtitle, diagnostic=False):
    st=states[step]
    setup_ax(ax,diagnostic)
    draw_tray(ax)
    draw_hopper(ax)
    draw_gate(ax,st['gate'])
    draw_beans(ax,st['qpos'])
    ax.set_title(title,loc="left",fontsize=12,fontweight="bold",color=INK,pad=23)
    ax.text2D(0,1.005,subtitle,transform=ax.transAxes,fontsize=9.2,color=MUTED,va="bottom")
    return st


# --------------------------------- figure ---------------------------------
plt.rcParams.update({"font.family":"DejaVu Sans","axes.titlecolor":INK,"text.color":INK})
fig=plt.figure(figsize=(17.5,11.4),facecolor=BG)
gs=GridSpec(2,3,figure=fig,left=.025,right=.98,top=.855,bottom=.105,wspace=.08,hspace=.24)
fig.text(.035,.956,"COFFEE / GRAVITY DISCHARGE",fontsize=12,fontweight="bold",color=TEAL)
fig.text(.035,.914,"50 beans. One 18 mm throat. Five seconds.",fontsize=27,fontweight="bold",color=INK)
fig.text(.965,.944,"MUJOCO  /  RIGID-BODY STUDY",ha="right",fontsize=10,color=MUTED)
fig.text(.965,.920,"cutaway view · dimensions in mm",ha="right",fontsize=9,color=MUTED)

axes=[]
for j,t in enumerate(SNAPSHOT_TIMES):
    row,col=divmod(j,3)
    ax=fig.add_subplot(gs[row,col],projection="3d")
    step=int(round(t/DT)); st=states[step]
    if t==0:
        title="01  /  CHARGED"
        sub="t = 0.00 s   ·   gate closed   ·   0 through"
    elif t==1:
        title="02  /  RELEASE"
        sub=f"t = 1.00 s   ·   gate starts   ·   {st['count']} through"
    elif t==2:
        title="03  /  FLOW"
        sub=f"t = 2.00 s   ·   gate open   ·   {st['count']} through"
    else:
        title="04  /  COLLECTED"
        sub=f"t = 5.00 s   ·   {st['count']} / 50 through throat"
    draw_state(ax,step,title,sub)
    axes.append(ax)
    if t==0:
        # Throat dimension and bean anatomy callout.
        ax.text2D(.02,.10,"Ø 18 mm\nTHROAT",transform=ax.transAxes,fontsize=9,fontweight="bold",color=COPPER)
        ax.text2D(.70,.83,"2 ELLIPSOIDS\n1 FREE BODY",transform=ax.transAxes,fontsize=8.5,color=MUTED,ha="center")

# Diagnostic occupies the sixth cell.
axd=fig.add_subplot(gs[1,1],projection="3d")
axd.computed_zorder=False
st=draw_state(axd,diag_step,"05  /  CONTACT DIAGNOSTIC",f"t = {DIAG_TIME:.2f} s   ·   lower hopper detail",True)
if contacts_diag:
    pts=np.array([c[0] for c in contacts_diag])
    axd.scatter(pts[:,0],pts[:,1],pts[:,2],s=28,c="#ffd45e",edgecolors="#704817",
                linewidths=.45,depthshade=False,zorder=30)
    forces=np.array([max(0.0,c[2][0]) for c in contacts_diag])
    scale=1.5  # drawn metres per newton, arrows capped for legibility
    origins,vectors=[],[]
    for (pos,frame,f6,g1,g2),fn in zip(contacts_diag,forces):
        if fn < 1e-6: continue
        direction=frame[0,:]
        length=min(.022,fn*scale)
        origins.append(pos); vectors.append(direction*length)
    if origins:
        oo=np.array(origins); vv=np.array(vectors)
        axd.quiver(oo[:,0],oo[:,1],oo[:,2],vv[:,0],vv[:,1],vv[:,2],
                   color="#007f89",linewidths=1.4,arrow_length_ratio=.28,zorder=31)
    fmax=forces.max()
else:
    fmax=0.0
axd.text2D(.01,.015,f"●  {len(contacts_diag or [])} contact points    →  normal force\nmax |Fn| = {fmax:.4f} N   ·   vector scale 1.5 m/N (capped)",
           transform=axd.transAxes,fontsize=8.6,color=INK,zorder=40,
           bbox=dict(boxstyle="round,pad=0.35",facecolor=BG,edgecolor="none",alpha=.90))

# Cumulative and individual crossing plot.
axp=fig.add_subplot(gs[1,2])
axp.set_facecolor("#eee9df")
times=np.arange(steps+1)*DT
axp.fill_between(times,0,count_history,color=TEAL,alpha=.12,step="post")
axp.step(times,count_history,where="post",color=TEAL,lw=2.2)
valid=cross_step>=0
event_t=cross_step[valid]*DT
ids=np.nonzero(valid)[0]+1
axp.scatter(event_t,ids,s=16,c="#8d4b2b",edgecolors="white",linewidths=.3,zorder=4)
for t in SNAPSHOT_TIMES[1:]:
    axp.axvline(t,color="#9e9689",lw=.8,ls=(0,(3,3)))
    axp.text(t,49.0,f" {t:g} s",fontsize=8,color=MUTED,va="top")
axp.axvline(GATE_TIME,color=GOLD,lw=1.6)
axp.text(1.05,44,"gate opens",fontsize=8.5,color="#8b5a16")
axp.set_xlim(0,5); axp.set_ylim(0,52); axp.set_xticks(range(6)); axp.set_yticks([0,10,20,30,40,50])
axp.grid(axis="y",color="white",lw=.9)
for s in ["top","right","left"]: axp.spines[s].set_visible(False)
axp.spines["bottom"].set_color("#b9b1a4")
axp.tick_params(axis="both",labelsize=8.5,colors=MUTED,length=0)
axp.set_xlabel("TIME  /  s",fontsize=9,labelpad=8,color=MUTED)
axp.set_ylabel("CUMULATIVE BEANS + BEAN ID AT CROSSING",fontsize=8.5,color=MUTED,labelpad=9)
axp.set_title("06  /  THROUGH THE THROAT",loc="left",fontsize=12,fontweight="bold",pad=18)
axp.text(0,1.025,"step: cumulative count   ·   dots: individual bean crossing",transform=axp.transAxes,fontsize=9.2,color=MUTED)
axp.text(.98,.08,f"{int(crossed.sum())} / 50\nby 5 seconds",transform=axp.transAxes,ha="right",fontsize=15,fontweight="bold",color=TEAL)

fig.text(.035,.061,
         f"BEAN   10.0 × 7.2 × 5.2 mm  ·  0.18 g  ·  two overlapping ellipsoids / free joint     "
         f"FRICTION   bean 0.45  ·  wall 0.28  ·  tray 0.60",
         fontsize=9.5,color=INK)
fig.text(.035,.035,
         f"SOLVER   Newton / elliptic cone  ·  implicitfast  ·  Δt {DT*1000:.2f} ms  ·  100 iterations  ·  tolerance 10⁻¹⁰     "
         f"WALL PENETRATION   max {max_penetration*1000:.3f} mm     ESCAPED   {len(left_ids)}",
         fontsize=9.5,color=MUTED)
fig.text(.965,.012,"Downward center-plane crossing at throat outlet · deterministic seed 19",ha="right",fontsize=8,color=MUTED)

out=os.environ.get("OUT",str(Path(__file__).resolve().with_name("hopper.png")))
fig.savefig(out,dpi=175,facecolor=BG,bbox_inches="tight")
plt.close(fig)