"""Geometry helpers for the sliding pressure door (models/pressure-door-sliding): the hinged door's pd_lib.py
plus loop, clipping and loft helpers at the end. Run inside Blender (exec'd by sd_build.py).

Everything is built as quads: surfaces of revolution (lathe), tubes swept along paths, rounded boxes made like a
cube-sphere (no bevel corner patches), and slabs with a rounded front edge. Only caps are triangle fans; flat
n-gon caps are triangulated in finish(). No booleans anywhere.
Each face is oriented as it is made, against an 'outward' direction the builder knows, so open surfaces are right too.
"""
import bpy, bmesh, math
from mathutils import Vector, Matrix

TAU = 2.0 * math.pi
I4 = Matrix.Identity(4)


def V(*a):
    return Vector(a)


def face(bm, verts, mat=0, out=None, ref=None):
    """make a face and turn it so its normal agrees with `out` (a direction) or points away from `ref` (a point)"""
    uq = []
    for v in verts:
        if v not in uq:
            uq.append(v)
    if len(uq) < 3:
        return None
    try:
        f = bm.faces.new(uq)
    except ValueError:
        return None
    f.material_index = mat
    f.normal_update()
    if out is None and ref is not None:
        out = f.calc_center_median() - Vector(ref)
    if out is not None and f.normal.dot(out) < 0:
        f.normal_flip()
    return f


# ------------------------------------------------------------------ 2D profiles
def fillet2(pts, radii, seg=3, closed=False):
    """round the corners of a 2D polyline. radii: one number or one per point (0 = keep the corner sharp)"""
    P = [Vector(p) for p in pts]
    n = len(P)
    if not isinstance(radii, (list, tuple)):
        radii = [radii] * n
    out = []
    for i in range(n):
        if not closed and (i == 0 or i == n - 1):
            out.append(P[i].copy())
            continue
        r = radii[i]
        p0, p1, p2 = P[i - 1], P[i], P[(i + 1) % n]
        if r <= 0:
            out.append(p1.copy())
            continue
        d1, d2 = p0 - p1, p2 - p1
        l1, l2 = d1.length, d2.length
        d1.normalize()
        d2.normalize()
        ang = math.acos(max(-1.0, min(1.0, d1.dot(d2))))
        if ang < 1e-3 or ang > math.pi - 1e-3:
            out.append(p1.copy())
            continue
        t = min(r / math.tan(ang / 2), l1 * 0.5, l2 * 0.5)
        re = t * math.tan(ang / 2)
        a, b = p1 + d1 * t, p1 + d2 * t
        c = p1 + (d1 + d2).normalized() * (re / math.sin(ang / 2))
        va, vb = a - c, b - c
        sweep = math.atan2(va.x * vb.y - va.y * vb.x, va.dot(vb))
        for k in range(seg + 1):
            out.append(c + Matrix.Rotation(sweep * k / seg, 2) @ va)
    clean = [out[0]]
    for p in out[1:]:
        if (p - clean[-1]).length > 1e-6:
            clean.append(p)
    if closed and len(clean) > 2 and (clean[0] - clean[-1]).length < 1e-6:
        clean.pop()
    return clean


def prof(spec, seg=3):
    """spec: [(r, h, fillet)] -> filleted 2D profile"""
    return fillet2([(a, b) for a, b, _ in spec], [c for _, _, c in spec], seg)


# ------------------------------------------------------------------ lathe
AX = {'Y': (V(1, 0, 0), V(0, 0, 1), V(0, 1, 0)),     # angle 0 = +X, 90 deg = +Z  (seen from the front, -Y, counter-clockwise)
      'Z': (V(1, 0, 0), V(0, 1, 0), V(0, 0, 1)),
      'X': (V(0, 1, 0), V(0, 0, 1), V(1, 0, 0))}


def lathe(bm, profile, n, center=(0, 0, 0), axis='Y', a0=0.0, a1=TAU, mat=0, M=None):
    """revolve profile [(r, h)] about the axis through center. The outside must lie on the LEFT of the direction of
    travel in (r, h), i.e. along (-dh, dr). Points with r == 0 sit on the axis (fans)."""
    u, v, w = AX[axis]
    C = Vector(center)
    full = abs((a1 - a0) - TAU) < 1e-9
    ncol = n if full else n + 1
    angs = [a0 + (a1 - a0) * j / n for j in range(ncol)]
    R3 = M.to_3x3() if M else None
    rows = []
    for (r, h) in profile:
        if r < 1e-7:
            p = C + w * h
            vv = bm.verts.new(M @ p if M else p)
            rows.append([vv] * ncol)
        else:
            row = []
            for a in angs:
                p = C + (u * math.cos(a) + v * math.sin(a)) * r + w * h
                row.append(bm.verts.new(M @ p if M else p))
            rows.append(row)
    faces = []
    for i in range(len(profile) - 1):
        (r0, h0), (r1, h1) = profile[i], profile[i + 1]
        dr, dh = r1 - r0, h1 - h0
        if abs(dr) < 1e-9 and abs(dh) < 1e-9:
            continue
        for j in range(n):
            j2 = (j + 1) % ncol
            am = a0 + (a1 - a0) * (j + 0.5) / n
            rad = u * math.cos(am) + v * math.sin(am)
            out = rad * (-dh) + w * dr
            if R3:
                out = R3 @ out
            f = face(bm, [rows[i][j], rows[i + 1][j], rows[i + 1][j2], rows[i][j2]], mat, out=out)
            if f:
                faces.append(f)
    return faces, rows


# ------------------------------------------------------------------ tubes
def arc_pts(center, radius, a0, a1, n, axis='Y'):
    u, v, w = AX[axis]
    C = Vector(center)
    return [C + (u * math.cos(a0 + (a1 - a0) * k / n) + v * math.sin(a0 + (a1 - a0) * k / n)) * radius for k in range(n + 1)]


def fillet3(pts, r, seg=4):
    """round the corners of a 3D polyline with arcs of radius r"""
    P = [Vector(p) for p in pts]
    out = [P[0]]
    for i in range(1, len(P) - 1):
        p0, p1, p2 = P[i - 1], P[i], P[i + 1]
        d1, d2 = (p0 - p1), (p2 - p1)
        l1, l2 = d1.length, d2.length
        d1.normalize()
        d2.normalize()
        ang = math.acos(max(-1.0, min(1.0, d1.dot(d2))))
        if ang > math.radians(172.0):        # nearly straight: keep the point (arcs come already sampled)
            out.append(p1)
            continue
        t = min(r / math.tan(ang / 2), l1 * 0.5, l2 * 0.5)
        re = t * math.tan(ang / 2)
        a, b = p1 + d1 * t, p1 + d2 * t
        c = p1 + (d1 + d2).normalized() * (re / math.sin(ang / 2))
        va, vb = a - c, b - c
        ax = va.cross(vb)
        th = va.angle(vb)
        for k in range(seg + 1):
            out.append(c + Matrix.Rotation(th * k / seg, 3, ax.normalized()) @ va)
    out.append(P[-1])
    return out


def _frames(path, closed):
    n = len(path)
    T = []
    for i in range(n):
        if closed:
            t = path[(i + 1) % n] - path[i - 1]
        else:
            t = path[min(i + 1, n - 1)] - path[max(i - 1, 0)]
        T.append(t.normalized())
    t0 = T[0]
    ref = V(0, 0, 1) if abs(t0.z) < 0.9 else V(1, 0, 0)
    N = [(ref - t0 * ref.dot(t0)).normalized()]
    for i in range(1, n):
        b = T[i - 1].cross(T[i])
        nn = N[-1].copy()
        if b.length > 1e-9:
            nn = Matrix.Rotation(T[i - 1].angle(T[i]), 3, b.normalized()) @ nn
        nn = (nn - T[i] * nn.dot(T[i])).normalized()
        N.append(nn)
    return T, N


def tube(bm, path, radius, seg=12, closed=False, caps=('dome', 'dome'), mat=0, cap_mat=None):
    """sweep a circle along path. caps: 'dome', 'flat' or None for each end"""
    path = [Vector(p) for p in path]
    T, N = _frames(path, closed)
    rings = []
    for i, p in enumerate(path):
        B = T[i].cross(N[i])
        rings.append([bm.verts.new(p + (N[i] * math.cos(k * TAU / seg) + B * math.sin(k * TAU / seg)) * radius) for k in range(seg)])
    faces = []
    nseg = len(path) if closed else len(path) - 1
    for i in range(nseg):
        i2 = (i + 1) % len(path)
        c = (path[i] + path[i2]) * 0.5
        for k in range(seg):
            k2 = (k + 1) % seg
            f = face(bm, [rings[i][k], rings[i2][k], rings[i2][k2], rings[i][k2]], mat, ref=c)
            if f:
                faces.append(f)
    if not closed:
        cm = mat if cap_mat is None else cap_mat
        for end, sgn in ((0, -1.0), (len(path) - 1, 1.0)):
            kind = caps[0] if end == 0 else caps[1]
            if not kind:
                continue
            p, t = path[end], T[end] * sgn
            B = T[end].cross(N[end])
            prev = rings[end]
            if kind == 'dome':
                steps = 3
                for s in range(1, steps):
                    th = (math.pi / 2) * s / steps
                    cc = p + t * radius * math.sin(th)
                    rr = radius * math.cos(th)
                    ring = [bm.verts.new(cc + (N[end] * math.cos(k * TAU / seg) + B * math.sin(k * TAU / seg)) * rr) for k in range(seg)]
                    for k in range(seg):
                        k2 = (k + 1) % seg
                        faces.append(face(bm, [prev[k], ring[k], ring[k2], prev[k2]], cm, ref=p))
                    prev = ring
                tip = bm.verts.new(p + t * radius)
            else:
                tip = bm.verts.new(p)
            for k in range(seg):
                k2 = (k + 1) % seg
                faces.append(face(bm, [prev[k], tip, prev[k2]], cm, out=t))
    return [f for f in faces if f]


# ------------------------------------------------------------------ rounded box (cube-sphere layout, all quads)
def rbox(bm, size, r, seg=3, M=I4, mat=0, drop=(), flip=False):
    """box of `size` centred on the origin of M, every edge rounded with r. drop: faces left out, e.g. ('+y',)"""
    h = [s * 0.5 for s in size]
    r = max(0.0, min(r, *h))
    inner = [hh - r for hh in h]

    def ticks(ax):
        hh, inn = h[ax], inner[ax]
        if r < 1e-7 or seg == 0:
            return [-hh, hh]
        band = [inn + r * math.tan(math.pi / 4 * k / seg) for k in range(seg + 1)]
        t = sorted([-b for b in band] + band)
        out = []
        for x in t:
            if not out or abs(x - out[-1]) > 1e-7:
                out.append(x)
        return out

    TK = [ticks(0), ticks(1), ticks(2)]
    cache = {}
    R3 = M.to_3x3()

    def vert(p):
        key = tuple(round(c, 7) for c in p)
        if key in cache:
            return cache[key]
        q = Vector([max(-inner[k], min(inner[k], p[k])) for k in range(3)])
        d = Vector(p) - q
        pp = q + d.normalized() * r if d.length > 1e-9 else Vector(p)
        vv = bm.verts.new(M @ pp)
        cache[key] = vv
        return vv

    faces = []
    for ax in range(3):
        a1, a2 = [k for k in range(3) if k != ax]
        for sgn in (-1, 1):
            if (('-' if sgn < 0 else '+') + 'xyz'[ax]) in drop:
                continue
            nrm = Vector((0, 0, 0))
            nrm[ax] = sgn * (-1 if flip else 1)
            out = R3 @ nrm
            for i in range(len(TK[a1]) - 1):
                for j in range(len(TK[a2]) - 1):
                    cs = []
                    for (ii, jj) in ((i, j), (i + 1, j), (i + 1, j + 1), (i, j + 1)):
                        p = [0.0, 0.0, 0.0]
                        p[ax] = sgn * h[ax]
                        p[a1] = TK[a1][ii]
                        p[a2] = TK[a2][jj]
                        cs.append(vert(p))
                    f = face(bm, cs, mat, out=out)
                    if f:
                        faces.append(f)
    return faces


# ------------------------------------------------------------------ slab: a 2D outline pushed out with a rounded front edge
def _inset(loop, d):
    """offset a closed CCW 2D loop inward by d (miter)"""
    n = len(loop)
    out = []
    for i in range(n):
        p0, p1, p2 = loop[i - 1], loop[i], loop[(i + 1) % n]
        e1 = (p1 - p0).normalized()
        e2 = (p2 - p1).normalized()
        n1 = Vector((-e1.y, e1.x))
        n2 = Vector((-e2.y, e2.x))
        m = (n1 + n2)
        if m.length < 1e-9:
            m = n1
        m.normalize()
        c = max(0.3, m.dot(n1))
        out.append(p1 + m * (d / c))
    return out


def slab(bm, outline, depth, r, seg=3, origin=(0, 0, 0), basis=None, mat=0, back=False, side_mat=None, cap=True):
    """outline: closed 2D loop, counter-clockwise seen from the front. basis (u, v, w): w = front direction.
    Pushed from w = 0 to w = depth; the front edge rounded with r. back=True closes the back too."""
    u, v, w = basis or (V(1, 0, 0), V(0, 0, 1), V(0, -1, 0))
    O = Vector(origin)
    loop = [Vector(p) for p in outline]
    r = min(r, depth)
    levels = [(0.0, 0.0)] + ([(0.0, depth - r)] if depth - r > 1e-6 else [])
    for k in range(1, seg + 1):
        th = (math.pi / 2) * k / seg
        levels.append((r * (1 - math.cos(th)), depth - r + r * math.sin(th)))
    if r < 1e-6:
        levels = [(0.0, 0.0), (0.0, depth)]
    rings = []
    for off, ww in levels:
        L = _inset(loop, off) if off > 1e-9 else loop
        rings.append([bm.verts.new(O + u * p.x + v * p.y + w * ww) for p in L])
    n = len(loop)
    faces = []
    sm = mat if side_mat is None else side_mat
    for i in range(len(rings) - 1):
        for k in range(n):
            k2 = (k + 1) % n
            e = loop[k2] - loop[k]
            o2 = Vector((e.y, -e.x))             # outward of a CCW loop
            out = u * o2.x + v * o2.y + w * (1.0 if i >= 1 else 0.0) * 0.3
            f = face(bm, [rings[i][k], rings[i][k2], rings[i + 1][k2], rings[i + 1][k]], sm, out=out)
            if f:
                faces.append(f)
    if cap:
        f = face(bm, rings[-1], mat, out=w)
        if f:
            faces.append(f)
    if back:
        f = face(bm, rings[0], mat, out=-w)
        if f:
            faces.append(f)
    return faces


def rrect(w, h, r, seg=3, cx=0.0, cy=0.0):
    """rounded rectangle outline, CCW"""
    x0, x1, y0, y1 = cx - w / 2, cx + w / 2, cy - h / 2, cy + h / 2
    return fillet2([(x0, y0), (x1, y0), (x1, y1), (x0, y1)], r, seg, closed=True)


# ------------------------------------------------------------------ floor
def squash_to_floor(bm, cx, cz, zf=0.0, verts=None):
    """move every vertex below the floor along its ray from the door axis (cx, cz) in XZ onto the floor"""
    for v in (verts if verts is not None else bm.verts):
        x, y, z = v.co
        if z < zf - 1e-9:
            t = (zf - cz) / (z - cz)
            v.co = Vector((cx + (x - cx) * t, y, zf))


def drop_floor_faces(bm, zf=0.0, eps=1e-5):
    """delete faces that lie in the floor plane (hidden contact) and faces with no area"""
    kill = [f for f in bm.faces if all(abs(v.co.z - zf) < eps for v in f.verts) or f.calc_area() < 1e-9]
    bmesh.ops.delete(bm, geom=kill, context='FACES')
    loose = [v for v in bm.verts if not v.link_faces]
    bmesh.ops.delete(bm, geom=loose, context='VERTS')


# ------------------------------------------------------------------ objects
def mat_rot_y(phi):
    """rotation that takes +X to the direction at angle phi in the XZ plane (counter-clockwise seen from the front)"""
    return Matrix.Rotation(-phi, 4, 'Y')


def finish(bm, name, coll, mats, origin=(0, 0, 0), sharp_deg=38.0, merge=1e-5):
    bmesh.ops.remove_doubles(bm, verts=list(bm.verts), dist=merge)
    bmesh.ops.dissolve_degenerate(bm, edges=list(bm.edges), dist=1e-7)
    ngons = [f for f in bm.faces if len(f.verts) > 4]
    if ngons:
        bmesh.ops.triangulate(bm, faces=ngons, quad_method='BEAUTY', ngon_method='BEAUTY')
    me = bpy.data.meshes.get(name)
    if me is None:
        me = bpy.data.meshes.new(name)
    else:
        me.clear_geometry()
        me.materials.clear()
    bm.to_mesh(me)
    bm.free()
    me.transform(Matrix.Translation(-Vector(origin)))
    used = sorted({p.material_index for p in me.polygons})
    remap = {old: new for new, old in enumerate(used)}
    for p in me.polygons:
        p.material_index = remap[p.material_index]
    for old in used:
        me.materials.append(mats[old])
    me.shade_smooth()
    me.set_sharp_from_angle(angle=math.radians(sharp_deg))
    ob = bpy.data.objects.get(name)
    if ob is None:
        ob = bpy.data.objects.new(name, me)
    else:
        ob.data = me
    if ob.name not in coll.objects:
        for c in list(ob.users_collection):
            c.objects.unlink(ob)
        coll.objects.link(ob)
    ob.location = origin
    return ob


def audit(ob):
    me = ob.data
    q = sum(1 for p in me.polygons if len(p.vertices) == 4)
    t = sum(1 for p in me.polygons if len(p.vertices) == 3)
    ng = sum(1 for p in me.polygons if len(p.vertices) > 4)
    tris = sum(len(p.vertices) - 2 for p in me.polygons)
    return dict(name=ob.name, quads=q, tris=t, ngons=ng, game_tris=tris, verts=len(me.vertices), mats=[m.name for m in me.materials])


# ================================================================== added for the sliding door (sd_lib_extra.py)
def rr_loop(cx, cz, hw, hh, r, seg=8):
    """rounded rectangle, counter-clockwise in (x, z), 4 * (seg + 1) points; equal seg -> equal counts (loft-able).
    r may be 0 (the corner point then repeats; remove_doubles cleans it)."""
    r = max(0.0, min(r, hw, hh))
    pts = []
    for k, (sx, sz) in enumerate(((1, -1), (1, 1), (-1, 1), (-1, -1))):      # corner centres, CCW from the lower right
        ccx, ccz = cx + sx * (hw - r), cz + sz * (hh - r)
        a0 = math.radians(-90 + 90 * k)
        for s in range(seg + 1):
            a = a0 + (math.pi / 2) * s / seg
            pts.append(Vector((ccx + r * math.cos(a), ccz + r * math.sin(a))))
    return pts


def clip_poly(subject, clip):
    """Sutherland-Hodgman: subject polygon (2D Vectors) clipped to a CONVEX, CCW clip polygon"""
    out = [Vector(p) for p in subject]
    n = len(clip)
    for i in range(n):
        a, b = Vector(clip[i]), Vector(clip[(i + 1) % n])
        e = b - a

        def inside(p):
            return e.x * (p.y - a.y) - e.y * (p.x - a.x) >= -1e-9

        def inter(p, q):
            d = q - p
            den = e.x * d.y - e.y * d.x
            if abs(den) < 1e-12:
                return p
            t = (e.x * (a.y - p.y) - e.y * (a.x - p.x)) / den
            return p + d * t
        inp = out
        out = []
        if not inp:
            break
        s = inp[-1]
        for p in inp:
            if inside(p):
                if not inside(s):
                    out.append(inter(s, p))
                out.append(p)
            elif inside(s):
                out.append(inter(s, p))
            s = p
    clean = []
    for p in out:
        if not clean or (p - clean[-1]).length > 1e-7:
            clean.append(p)
    if len(clean) > 1 and (clean[0] - clean[-1]).length < 1e-7:
        clean.pop()
    return clean


def halfplane(side, xc, big=50.0):
    """a CCW rectangle covering x <= xc (side -1) or x >= xc (side +1), for clip_poly"""
    if side < 0:
        return [Vector((-big, -big)), Vector((xc, -big)), Vector((xc, big)), Vector((-big, big))]
    return [Vector((xc, -big)), Vector((big, -big)), Vector((big, big)), Vector((xc, big))]


def rect_poly(cx, cz, hw, hh):
    return [Vector((cx - hw, cz - hh)), Vector((cx + hw, cz - hh)), Vector((cx + hw, cz + hh)), Vector((cx - hw, cz + hh))]


def rot_poly(poly, ang, cx=0.0, cz=0.0):
    R = Matrix.Rotation(ang, 2)
    C = Vector((cx, cz))
    return [C + R @ (Vector(p) - C) for p in poly]


def notch_edge(poly, xc, cz, R, side, seg=12):
    """on a polygon with a straight vertical edge at x = xc, replace the middle of that edge by a half-circle notch of
    radius R about (0, cz) bulging into the polygon (side -1: the polygon lies at x <= xc; +1: at x >= xc)"""
    a = math.sqrt(max(R * R - xc * xc, 0.0))
    n = len(poly)
    for i in range(n):
        p, q = poly[i], poly[(i + 1) % n]
        if abs(p.x - xc) < 1e-7 and abs(q.x - xc) < 1e-7 and min(p.y, q.y) < cz - a and max(p.y, q.y) > cz + a:
            up = q.y > p.y
            if side < 0:
                angs = [math.radians(270 - 180 * k / seg) for k in range(seg + 1)]        # bottom, left, top
            else:
                angs = [math.radians(-90 + 180 * k / seg) for k in range(seg + 1)]        # bottom, right, top
            arc = [Vector((R * math.cos(t), cz + R * math.sin(t))) for t in angs]
            arc[0] = Vector((xc, cz - a))
            arc[-1] = Vector((xc, cz + a))
            if not up:
                arc.reverse()
            return poly[:i + 1] + arc + poly[i + 1:]
    raise ValueError("no straight edge at x = %g spanning the notch" % xc)


def rect_from_loop(loop, cx, cz, x0, x1, z0, z1):
    """project every point of a loop by a ray from (cx, cz) onto the rectangle's boundary; the four points nearest
    the corners snap onto the corners, so a loft from the loop to the rectangle keeps the rectangle's corners"""
    out = []
    for p in loop:
        d = Vector(p) - Vector((cx, cz))
        best = None
        for (ex, ez) in ((x0, None), (x1, None), (None, z0), (None, z1)):
            if ex is not None:
                if abs(d.x) < 1e-12:
                    continue
                t = (ex - cx) / d.x
            else:
                if abs(d.y) < 1e-12:
                    continue
                t = (ez - cz) / d.y
            if t <= 0:
                continue
            q = Vector((cx, cz)) + d * t
            if x0 - 1e-5 <= q.x <= x1 + 1e-5 and z0 - 1e-5 <= q.y <= z1 + 1e-5:
                if best is None or t < best[0]:
                    best = (t, q)
        out.append(best[1])
    for corner in (Vector((x0, z0)), Vector((x1, z0)), Vector((x1, z1)), Vector((x0, z1))):
        k = min(range(len(out)), key=lambda i: (out[i] - corner).length)
        out[k] = corner.copy()
    return out


def loft(bm, levels, mats=None, mat=0, center=(0.0, 0.0)):
    """levels: list of (loop2D, y, kind): equal point counts, consecutive levels bridged with quads. The band
    between level i and i+1 takes level i+1's kind: 'front' (-Y), 'back' (+Y), 'out' (away from center in XZ),
    'in' (toward center) or a Vector; and its material from mats[i] (or mat)."""
    rows = [[bm.verts.new(Vector((p.x, y, p.y))) for p in loop] for (loop, y, kind) in levels]
    C = Vector((center[0], center[1]))
    faces = []
    n = len(levels[0][0])
    for i in range(len(levels) - 1):
        kind = levels[i + 1][2]
        m = mats[i] if mats else mat
        for k in range(n):
            k2 = (k + 1) % n
            quad = [rows[i][k], rows[i][k2], rows[i + 1][k2], rows[i + 1][k]]
            c = sum((v.co for v in quad), Vector()) / 4.0
            if isinstance(kind, Vector):
                out = kind
            elif kind == 'front':
                out = Vector((0, -1, 0))
            elif kind == 'back':
                out = Vector((0, 1, 0))
            else:
                rad = Vector((c.x - C.x, 0.0, c.z - C.y))
                if rad.length < 1e-9:
                    rad = Vector((1, 0, 0))
                rad.normalize()
                out = rad if kind == 'out' else -rad
            f = face(bm, quad, m, out=out)
            if f:
                faces.append(f)
    return faces, rows


def round_levels(loop_fn, d1, y_in, y_out, r, seg=3, inward=False):
    """arc levels for a rounded edge: from the face at y_in (its loop at offset d1 - r) round to the wall at offset d1
    running toward y_out. loop_fn(d) gives the loop at offset d."""
    out = []
    sgn = 1.0 if y_out > y_in else -1.0
    for s in range(1, seg + 1):
        th = (math.pi / 2) * s / seg
        d = (d1 - r) + r * math.sin(th)
        y = y_in + sgn * r * (1 - math.cos(th))
        out.append((loop_fn(d), y, 'in' if inward else 'out'))
    return out


def poly_slab(bm, poly, y0, y1, r=0.004, seg=3, mat=0, side_mat=None, back=True):
    """a 2D polygon (CCW in x, z) pushed from y0 to y1 (the front is at y1), the front edge rounded with r"""
    depth = abs(y1 - y0)
    w = Vector((0, -1, 0)) if y1 < y0 else Vector((0, 1, 0))
    return slab(bm, [Vector(p) for p in poly], depth, r, seg, origin=(0, y0, 0), basis=(V(1, 0, 0), V(0, 0, 1), w),
                mat=mat, back=back, side_mat=side_mat)


def box_at(bm, center, size, r=0.003, seg=2, mat=0, rot=None):
    M = Matrix.Translation(Vector(center))
    if rot is not None:
        M = M @ rot
    return rbox(bm, size, r, seg, M=M, mat=mat)


def sag_path(p0, p1, sag, n=10, side=0.0):
    """a parabola between two points, hanging down by sag in the middle (and bulging `side` toward -Y)"""
    P0, P1 = Vector(p0), Vector(p1)
    pts = []
    for k in range(n + 1):
        t = k / n
        p = P0.lerp(P1, t)
        h = 4 * t * (1 - t)
        pts.append(Vector((p.x, p.y - side * h, p.z - sag * h)))
    return pts