"""
Parts made for this hull (not kit): ray-cast fitting against the ship's shell, bulkheads with door openings,
decks, slabs and the carousel pylon. Quad grids, welded, UV-mapped; triangles only where a grid meets the curved
shell. No booleans.
"""
import math

import bmesh
import bpy
from mathutils import Vector
from mathutils.bvhtree import BVHTree


class Fit(object):
    """Ray casts against named objects of the open file (their evaluated meshes, world space)."""

    def __init__(self, names):
        dg = bpy.context.evaluated_depsgraph_get()
        self.trees = {}
        for n in names:
            ob = bpy.data.objects.get(n)
            if ob is None:
                continue
            ev = ob.evaluated_get(dg)
            me = ev.to_mesh()
            mw = ob.matrix_world
            verts = [tuple(mw @ v.co) for v in me.vertices]
            polys = [tuple(p.vertices) for p in me.polygons]
            self.trees[n] = BVHTree.FromPolygons(verts, polys, epsilon=0.0)
            ev.to_mesh_clear()

    def ray(self, o, d, only=None, dist=200.0):
        best = None
        for n, t in self.trees.items():
            if only is not None and n not in only:
                continue
            loc, nrm, idx, dd = t.ray_cast(Vector(o), Vector(d), dist)
            if loc is not None and (best is None or dd < best[1]):
                best = (loc, dd)
        return best[0] if best else None

    def side_y(self, x, z, sgn, only=None, y0=0.0):
        p = self.ray((x, y0, z), (0, sgn, 0), only)
        return None if p is None else p.y

    def ceil_z(self, x, y, z0, only=None):
        p = self.ray((x, y, z0), (0, 0, 1), only)
        return None if p is None else p.z

    def floor_z(self, x, y, z0, only=None):
        p = self.ray((x, y, z0), (0, 0, -1), only)
        return None if p is None else p.z


# ---------------------------------------------------------------------------------------------------------------
def _emit(bm, vm, co):
    key = (round(co[0], 4), round(co[1], 4), round(co[2], 4))
    v = vm.get(key)
    if v is None:
        v = bm.verts.new(co)
        vm[key] = v
    return v


def fix_tjunctions(bm, tol=1e-5):
    """Where a grid steps, a long boundary edge can meet two short ones: insert the missing vertices and re-split
    that face, so the sheet is properly welded before it gets its thickness."""
    for _ in range(4):
        todo = []
        for e in [e for e in bm.edges if len(e.link_faces) == 1]:
            a = e.verts[0].co
            ab = e.verts[1].co - a
            l2 = ab.length_squared
            if l2 < 1e-10:
                continue
            inside = set()
            for v in bm.verts:
                if v is e.verts[0] or v is e.verts[1]:
                    continue
                ap = v.co - a
                t = ap.dot(ab) / l2
                if tol < t < 1 - tol and (ap - ab * t).length < tol:
                    inside.add(round(t, 7))
            if inside:
                todo.append((e, sorted(inside)))
        if not todo:
            break
        for e, ts in todo:
            if not e.is_valid:
                continue
            v0, v1 = e.verts[0], e.verts[1]
            a = v0.co.copy()
            b = v1.co.copy()
            cur, start, prev_t = e, v0, 0.0
            for t in ts:
                fac = (t - prev_t) / (1.0 - prev_t)
                ne, nv = bmesh.utils.edge_split(cur, start, fac)
                nv.co = a + (b - a) * t
                cur = next(x for x in nv.link_edges if v1 in x.verts)
                start, prev_t = nv, t
        bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-5)
    ng = [f for f in bm.faces if len(f.verts) > 4]
    if ng:
        r = bmesh.ops.triangulate(bm, faces=ng)
        try:
            bmesh.ops.join_triangles(bm, faces=r["faces"], angle_face_threshold=0.02, angle_shape_threshold=3.14)
        except Exception:  # noqa: BLE001
            pass


def _lines(a, b, step, extra=()):
    vals = {round(a, 5), round(b, 5)}
    k = math.ceil(a / step - 1e-9)
    while k * step < b - 1e-9:
        if k * step > a + 1e-9:
            vals.add(round(k * step, 5))
        k += 1
    for e in extra:
        if a - 1e-9 <= e <= b + 1e-9:
            vals.add(round(e, 5))
    return sorted(vals)


def _finish(name, bm, coll, parent, material, thick, axis, uv_scale=0.25):
    """Give the sheet its thickness, box UVs and one material; link the object."""
    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-5)
    fix_tjunctions(bm)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    if thick > 0:
        res = bmesh.ops.solidify(bm, geom=bm.faces[:], thickness=thick)
        bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    uv = bm.loops.layers.uv.verify()
    for f in bm.faces:
        n = f.normal
        ax = max(range(3), key=lambda i: abs(n[i]))
        for l in f.loops:
            c = l.vert.co
            l[uv].uv = ((c.y, c.z) if ax == 0 else (c.x, c.z) if ax == 1 else (c.x, c.y))
            l[uv].uv = (l[uv].uv[0] * uv_scale, l[uv].uv[1] * uv_scale)
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me)
    bm.free()
    if material is not None:
        me.materials.append(material)
    ob = bpy.data.objects.new(name, me)
    coll.objects.link(ob)
    if parent is not None:
        ob.parent = parent
    return ob


def bulkhead_x(name, coll, parent, fit, only, x_ray, x_a, x_b, zbot, ztop, openings=(), material=None, cell=1.2, inset=0.0):
    """Wall across the ship between x_a and x_b (its two faces), fitted to the shell found by rays cast at x_ray.
    openings: (y0, y1, z0, z1). The top closes onto the shell's ridge."""
    zs = _lines(zbot, ztop, cell, [v for o in openings for v in (o[2], o[3])])
    ext = {}
    for z in list(zs):
        lo, hi = fit.side_y(x_ray, z, -1, only), fit.side_y(x_ray, z, 1, only)
        if lo is None or hi is None or hi - lo < 0.3:
            zs.remove(z)
            continue
        ext[z] = (lo + inset, hi - inset)
    okeys = [v for o in openings for v in (o[0], o[1])]
    ymin = min(e[0] for e in ext.values())
    ymax = max(e[1] for e in ext.values())
    ylines = _lines(ymin, ymax, cell, okeys)
    bm = bmesh.new()
    vm = {}
    for z0, z1 in zip(zs[:-1], zs[1:]):
        (l0, r0), (l1, r1) = ext[z0], ext[z1]
        lo, hi = max(l0, l1) + 0.12, min(r0, r1) - 0.12
        mids = [y for y in ylines if lo < y < hi]
        c0, c1 = [l0] + mids + [r0], [l1] + mids + [r1]
        for j in range(len(c0) - 1):
            cz = (z0 + z1) / 2
            cy = (c0[j] + c0[j + 1] + c1[j] + c1[j + 1]) / 4
            if any(a0 - 1e-6 <= cy <= a1 + 1e-6 and b0 - 1e-6 <= cz <= b1 + 1e-6 for (a0, a1, b0, b1) in openings):
                continue
            q = [_emit(bm, vm, (x_a, c0[j], z0)), _emit(bm, vm, (x_a, c1[j], z1)),
                 _emit(bm, vm, (x_a, c1[j + 1], z1)), _emit(bm, vm, (x_a, c0[j + 1], z0))]
            try:
                bm.faces.new(q)
            except ValueError:
                pass
    # close onto the ridge above the last line
    zt = zs[-1]
    apex = fit.ceil_z(x_ray, 0.5 * (ext[zt][0] + ext[zt][1]), zt - 0.05, only)
    if apex is not None and apex - zt > 0.01:
        row = sorted([v for k, v in vm.items() if abs(k[2] - round(zt, 4)) < 1e-4], key=lambda v: v.co.y)
        top = _emit(bm, vm, (x_a, 0.5 * (ext[zt][0] + ext[zt][1]), apex))
        for a, b in zip(row[:-1], row[1:]):
            try:
                bm.faces.new((a, top, b))
            except ValueError:
                pass
    ob = _finish(name, bm, coll, parent, material, 0.0, 0)
    # thickness toward x_b (kept as a real second sheet, welded round the edges)
    me = ob.data
    bm = bmesh.new()
    bm.from_mesh(me)
    for f in bm.faces:
        f.normal_update()
    sign = 1.0 if x_b > x_a else -1.0
    flip = [f for f in bm.faces if f.normal.x * sign > 0]
    if flip:
        bmesh.ops.reverse_faces(bm, faces=flip)
    res = bmesh.ops.extrude_face_region(bm, geom=bm.faces[:], use_keep_orig=True)
    new_verts = [g for g in res["geom"] if isinstance(g, bmesh.types.BMVert)]
    bmesh.ops.translate(bm, verts=new_verts, vec=(x_b - x_a, 0, 0))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    uv = bm.loops.layers.uv.verify()
    for f in bm.faces:
        n = f.normal
        ax = max(range(3), key=lambda i: abs(n[i]))
        for l in f.loops:
            c = l.vert.co
            p = (c.y, c.z) if ax == 0 else (c.x, c.z) if ax == 1 else (c.x, c.y)
            l[uv].uv = (p[0] * 0.25, p[1] * 0.25)
    bm.to_mesh(me)
    bm.free()
    return ob


def deck_fitted(name, coll, parent, fit, only, x0, x1, ztop, thick, material=None, cell=1.2, inset=0.05, holes=(), half_max=None):
    """Deck between x0 and x1 fitted to the shell (narrowest of its top, middle and underside)."""
    xs = _lines(x0, x1, cell, [v for h in holes for v in (h[0], h[1])])
    ext = {}
    for x in xs:
        lo = hi = None
        for z in (ztop - 0.01, ztop - thick / 2, ztop - thick + 0.01):
            a, b = fit.side_y(x, z, -1, only), fit.side_y(x, z, 1, only)
            if a is None or b is None:
                raise RuntimeError("%s: no shell at x=%.2f z=%.2f" % (name, x, z))
            lo = a if lo is None else max(lo, a)
            hi = b if hi is None else min(hi, b)
        lo, hi = lo + inset, hi - inset
        if half_max is not None:
            lo, hi = max(lo, -half_max), min(hi, half_max)
        ext[x] = (lo, hi)
    hkeys = [v for h in holes for v in (h[2], h[3])]
    ylines = _lines(min(e[0] for e in ext.values()), max(e[1] for e in ext.values()), cell, hkeys)
    bm = bmesh.new()
    vm = {}
    for xa, xb in zip(xs[:-1], xs[1:]):
        (l0, r0), (l1, r1) = ext[xa], ext[xb]
        lo, hi = max(l0, l1) + 0.12, min(r0, r1) - 0.12
        mids = [y for y in ylines if lo < y < hi]
        c0, c1 = [l0] + mids + [r0], [l1] + mids + [r1]
        for j in range(len(c0) - 1):
            cx = (xa + xb) / 2
            cy = (c0[j] + c0[j + 1] + c1[j] + c1[j + 1]) / 4
            if any(h[0] - 1e-6 <= cx <= h[1] + 1e-6 and h[2] - 1e-6 <= cy <= h[3] + 1e-6 for h in holes):
                continue
            q = [_emit(bm, vm, (xa, c0[j], ztop)), _emit(bm, vm, (xb, c1[j], ztop)),
                 _emit(bm, vm, (xb, c1[j + 1], ztop)), _emit(bm, vm, (xa, c0[j + 1], ztop))]
            try:
                bm.faces.new(q)
            except ValueError:
                pass
    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-5)
    fix_tjunctions(bm)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    up = [f for f in bm.faces if f.normal.z < 0]
    if up:
        bmesh.ops.reverse_faces(bm, faces=up)
    res = bmesh.ops.extrude_face_region(bm, geom=bm.faces[:], use_keep_orig=True)
    bmesh.ops.translate(bm, verts=[g for g in res["geom"] if isinstance(g, bmesh.types.BMVert)], vec=(0, 0, -thick))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return _finish(name, bm, coll, parent, material, 0.0, 2)


def slab(name, coll, parent, x0, x1, y0, y1, ztop, thick, material=None, cell=1.2):
    """Plain rectangular slab as a quad grid."""
    xs, ys = _lines(x0, x1, cell), _lines(y0, y1, cell)
    bm = bmesh.new()
    vm = {}
    for xa, xb in zip(xs[:-1], xs[1:]):
        for ya, yb in zip(ys[:-1], ys[1:]):
            bm.faces.new([_emit(bm, vm, (xa, ya, ztop)), _emit(bm, vm, (xb, ya, ztop)),
                          _emit(bm, vm, (xb, yb, ztop)), _emit(bm, vm, (xa, yb, ztop))])
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    up = [f for f in bm.faces if f.normal.z < 0]
    if up:
        bmesh.ops.reverse_faces(bm, faces=up)
    res = bmesh.ops.extrude_face_region(bm, geom=bm.faces[:], use_keep_orig=True)
    bmesh.ops.translate(bm, verts=[g for g in res["geom"] if isinstance(g, bmesh.types.BMVert)], vec=(0, 0, -thick))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return _finish(name, bm, coll, parent, material, 0.0, 2)


def box(name, coll, parent, x0, x1, y0, y1, z0, z1, material=None, cell=1.2):
    """Plain box (a pier, a lintel): gridded top and bottom, one quad per grid step on each side."""
    return slab(name, coll, parent, min(x0, x1), max(x0, x1), min(y0, y1), max(y0, y1), max(z0, z1), abs(z1 - z0), material, cell)


def ledge_fitted(name, coll, parent, fit, only, x0, x1, y_in, sgn, ztop, thick, material=None, cell=1.2, inset=0.02):
    """Shelf along the ship from the line y = y_in (inner edge, straight) out to the shell on side sgn (+1 port)."""
    xs = _lines(x0, x1, cell)
    outer = {}
    for x in xs:
        v = None
        for z in (ztop - 0.005, ztop - thick / 2, ztop - thick + 0.005):
            y = fit.side_y(x, z, sgn, only)
            if y is None:
                raise RuntimeError("%s: no shell at x=%.2f z=%.2f" % (name, x, z))
            v = y if v is None else (min(v, y) if sgn > 0 else max(v, y))
        outer[x] = v - sgn * inset
    bm = bmesh.new()
    vm = {}
    for xa, xb in zip(xs[:-1], xs[1:]):
        q = [_emit(bm, vm, (xa, y_in, ztop)), _emit(bm, vm, (xb, y_in, ztop)),
             _emit(bm, vm, (xb, outer[xb], ztop)), _emit(bm, vm, (xa, outer[xa], ztop))]
        bm.faces.new(q)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    up = [f for f in bm.faces if f.normal.z < 0]
    if up:
        bmesh.ops.reverse_faces(bm, faces=up)
    res = bmesh.ops.extrude_face_region(bm, geom=bm.faces[:], use_keep_orig=True)
    bmesh.ops.translate(bm, verts=[g for g in res["geom"] if isinstance(g, bmesh.types.BMVert)], vec=(0, 0, -thick))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    return _finish(name, bm, coll, parent, material, 0.0, 2)


class FitFirst(object):
    """Fit against several shells in a fixed order of preference: the first one of `order` that a ray meets counts
    (the hull where it exists, the roof pod where the hull is open into it)."""

    def __init__(self, fit, order):
        self.fit, self.order = fit, order

    def side_y(self, x, z, sgn, only=None, y0=0.0):
        for n in self.order:
            v = self.fit.side_y(x, z, sgn, (n,), y0)
            if v is not None:
                return v
        return None

    def ceil_z(self, x, y, z0, only=None):
        for n in self.order:
            v = self.fit.ceil_z(x, y, z0, (n,))
            if v is not None:
                return v
        return None


def drawn_tree(ob, evaluated=False):
    """BVH of an object's surface AS BLENDER DRAWS IT: its own loop triangles, in world space.
    (Not BVHTree.FromPolygons on the polygons: that splits a quad along its own choice of diagonal, and on this
    hull's warped quads the two surfaces lie up to 11 mm apart. Measured against the wrong one, the roof pod's
    edge seemed 5 mm under the hull's skin all round and was in fact up to 3 mm above it in twenty places.)"""
    if evaluated:
        ev = ob.evaluated_get(bpy.context.evaluated_depsgraph_get())
        me = ev.to_mesh()
    else:
        me = ob.data
    me.calc_loop_triangles()
    mw = ob.matrix_world
    tree = BVHTree.FromPolygons([tuple(mw @ v.co) for v in me.vertices], [tuple(t.vertices) for t in me.loop_triangles], epsilon=0.0)
    if evaluated:
        ev.to_mesh_clear()
    return tree


class Envelope(object):
    """The ship's outer skin where two shells count as one: the hull, and the roof pod that sits on it. Their own
    meshes are used, without modifiers (the Mask on the hull only hides faces that lie under the pod), triangulated
    the way Blender draws them (drawn_tree). Hull skin that lies under the pod is not outer skin: the two spaces
    are one room there."""

    def __init__(self, hull="M_SHIP_Hull", pod="M_SHIP_Dorsal_Roof"):
        self.trees = {}
        for key, n in (("hull", hull), ("pod", pod)):
            ob = bpy.data.objects.get(n) if n else None
            if ob is None:
                continue
            self.trees[key] = drawn_tree(ob)

    def first(self, o, d, dist=400.0):
        """First skin met from o along d: (point, distance, 'hull' | 'pod'), or None."""
        best = None
        o, d = Vector(o), Vector(d)
        for k, t in self.trees.items():
            loc, nrm, idx, dd = t.ray_cast(o, d, dist)
            if loc is not None and (best is None or dd < best[1]):
                best = (loc, dd, k)
        return best

    def side_y(self, x, z, sgn):
        """The outermost skin on side sgn (+1 port) at height z, met coming from outside."""
        b = self.first((x, sgn * 100.0, z), (0.0, -sgn, 0.0))
        return None if b is None else b[0].y

    def top_z(self, x, y):
        """The outermost skin over (x, y), met coming from above."""
        b = self.first((x, y, 100.0), (0.0, 0.0, -1.0))
        return None if b is None else b[0].z

    def clearance(self, p):
        """Distance from p to the outer skin."""
        p = Vector(p)
        d = 1e9
        pod = self.trees.get("pod")
        if pod is not None:
            loc, nrm, idx, dd = pod.find_nearest(p)
            if loc is not None:
                d = dd
        hull = self.trees.get("hull")
        if hull is not None:
            loc, nrm, idx, dd = hull.find_nearest(p)
            if loc is not None and dd < d:
                under_pod = pod is not None and pod.ray_cast(loc + nrm * 0.001, nrm, 30.0)[0] is not None
                if not under_pod:
                    d = dd
        return d

    def inside(self, p):
        """Is p inside the ship: under the hull's skin, or above it with the pod overhead?
        (Not by counting rays that get away: the hull has real openings, the well and the canopy glazing.)"""
        p = Vector(p)
        hull, pod = self.trees.get("hull"), self.trees.get("pod")
        if hull is not None:
            loc, nrm, idx, dd = hull.find_nearest(p)
            if loc is not None and (p - loc).dot(nrm) < 0.0:
                return True
        if pod is not None:
            # (the second and third ray: one cast exactly along a mesh edge, the pod's centre seam, can slip through)
            for d in ((0.0, 0.0, 1.0), (0.0, 0.003, 1.0), (0.003, 0.0, 1.0)):
                if pod.ray_cast(p, Vector(d).normalized(), 60.0)[0] is not None:
                    return True
        return False


def _half_width(env, xs, z, inset, step=0.002):
    """Half the width of a wall at height z that keeps `inset` clear of the outer skin at every station of xs
    (both sides, the narrower one: the wall is symmetric)."""
    w = None
    for sgn in (1.0, -1.0):
        ys = [env.side_y(x, z, sgn) for x in xs]
        if any(v is None for v in ys):
            return None
        y = min(abs(v) for v in ys) - inset
        n = 0
        while y > 0.2 and min(env.clearance((x, sgn * y, z)) for x in xs) < inset:
            y -= step
            n += 1
            if n > 1500:
                return None
        if y <= 0.2:
            return None
        w = y if w is None else min(w, y)
    return w


def _top_height(env, xs, y, inset, step=0.002):
    """Height of a wall's top at +-y that keeps `inset` clear of the outer skin at every station of xs."""
    sides = (1.0, -1.0) if y > 1e-6 else (1.0,)
    zs = [env.top_z(x, s * y) for x in xs for s in sides]
    if any(v is None for v in zs):
        return None
    z = min(zs) - inset
    n = 0
    while min(env.clearance((x, s * y, z)) for x in xs for s in sides) < inset:
        z -= step
        n += 1
        if n > 1500:
            return None
    return z


def section_wall_mesh(name, env, x_a, x_b, zbot, openings=(), material=None, cell=1.2, inset=0.02,
                      tol_out=0.005, tol_gap=0.015, min_h=0.06, crown_drop=0.12, info=None):
    """Mesh of a wall across the ship, between its two faces x_a and x_b, that closes the whole section there: the
    hull and, above it, the roof pod (env: an Envelope). Its edge stays `inset` inside the outer skin, measured as
    a true distance at both faces and the middle, so nothing of it can show outside.

    The sheet is a grid of quads in rows. Row lines are put where the shape needs them: a straight edge between
    two lines may pass the true edge by no more than tol_out (that is what keeps it inside where the skin is
    hollow, at the corner between hull and pod and up the pod's foot) and may fall short of it by no more than
    tol_gap. The top row's cells reach up to the roof at the grid's column lines, where the roof is nearly level
    (within crown_drop of its highest point). openings: (y0, y1, z0, z1)."""
    xs = (min(x_a, x_b), 0.5 * (x_a + x_b), max(x_a, x_b))
    cache = {}

    def half(z):
        k = round(z, 4)
        if k not in cache:
            cache[k] = _half_width(env, xs, k, inset)
        return cache[k]

    h0 = half(zbot)
    if h0 is None:
        raise RuntimeError("%s: no shell at z=%.2f" % (name, zbot))
    okeys = [v for o in openings for v in (o[0], o[1])]
    cols = [y for y in _lines(-h0, h0, cell, okeys) if abs(y) < h0 - 0.01]

    # ---- the crown: the column lines over which the roof is nearly level ---------------------------------------
    apex = _top_height(env, xs, 0.0, inset)
    if apex is None:
        raise RuntimeError("%s: no roof over the centre line" % name)
    tops = {0.0: apex}
    yc = 0.0
    for y in sorted(v for v in cols if v > 1e-6):
        t = _top_height(env, xs, y, inset)
        if t is None or apex - t > crown_drop:
            break
        tops[round(y, 5)] = tops[round(-y, 5)] = t
        yc = y
    z_r = tops[round(yc, 5)]                       # the height at which the wall is 2 * yc wide

    # ---- row lines ---------------------------------------------------------------------------------------------
    # (only the lines that have to be there to begin with: the foot, the openings, the crown's base. The rest is
    # put where the shape asks for it, and no row is taller than a cell.)
    zs = sorted({round(zbot, 5), round(z_r, 5)} | {round(v, 5) for o in openings for v in (o[2], o[3]) if zbot < v < z_r})
    val = {}
    for z in zs:
        val[z] = half(z)
        if val[z] is None:
            raise RuntimeError("%s: no shell at z=%.2f" % (name, z))
    val[round(z_r, 5)] = max(yc, 1e-4)
    i = 0
    while i < len(zs) - 1:
        z0, z1 = zs[i], zs[i + 1]
        h = z1 - z0
        at = None
        if h >= 2.0 * min_h:
            n = max(2, int(math.ceil(h / 0.02)))
            worst = 1.0
            for k in range(1, n):
                zm = z0 + h * k / n
                ym = half(zm)
                if ym is None:
                    continue
                chord = val[z0] + (val[z1] - val[z0]) * k / n
                score = max((chord - ym) / tol_out, (ym - chord) / tol_gap)
                if score > worst:
                    worst, at = score, zm
            if at is not None:
                at = min(max(at, z0 + min_h), z1 - min_h)
            elif h > cell + 1e-6:
                at = z0 + h / math.ceil(h / cell - 1e-9)
        if at is not None:
            zm = round(at, 4)
            val[zm] = half(zm)
            zs.insert(i + 1, zm)
            continue
        i += 1
    # the row under the crown has to be wider than the crown by a usable cell
    while len(zs) > 2 and val[zs[-2]] < yc + 0.15:
        del zs[-2]

    # ---- the sheet ---------------------------------------------------------------------------------------------
    bm = bmesh.new()
    vm = {}

    def face(pts):
        try:
            bm.faces.new([_emit(bm, vm, (x_a, y, z)) for y, z in pts])
        except ValueError:
            pass

    def is_open(cy, cz):
        return any(a0 - 1e-6 <= cy <= a1 + 1e-6 and b0 - 1e-6 <= cz <= b1 + 1e-6 for (a0, a1, b0, b1) in openings)

    for z0, z1 in zip(zs[:-2], zs[1:-1]):
        w0, w1 = val[z0], val[z1]
        lim = min(w0, w1) - 0.12
        mids = [y for y in cols if -lim < y < lim]
        c0, c1 = [-w0] + mids + [w0], [-w1] + mids + [w1]
        for j in range(len(c0) - 1):
            if is_open((c0[j] + c0[j + 1] + c1[j] + c1[j + 1]) / 4.0, (z0 + z1) / 2.0):
                continue
            face([(c0[j], z0), (c1[j], z1), (c1[j + 1], z1), (c0[j + 1], z0)])
    # the top row: from the last line up to the roof
    z_p = zs[-2]
    w_p = val[z_p]
    lim = w_p - 0.12
    bottom = [-w_p] + [y for y in cols if -lim < y < lim] + [w_p]
    eps = 1e-6
    for a, b in zip(bottom[:-1], bottom[1:]):
        if is_open((a + b) / 2.0, (z_p + z_r) / 2.0):
            continue
        if a >= -yc - eps and b <= yc + eps and b - a > eps:
            face([(a, z_p), (a, tops[round(a, 5)]), (b, tops[round(b, 5)]), (b, z_p)])
        elif b <= -yc + eps:
            face([(a, z_p), (-yc, z_r), (b, z_p)])
        else:
            face([(a, z_p), (yc, z_r), (b, z_p)])

    bmesh.ops.remove_doubles(bm, verts=bm.verts[:], dist=1e-5)
    fix_tjunctions(bm)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    # thickness toward x_b: a real second sheet, welded round the edge and round the openings
    for f in bm.faces:
        f.normal_update()
    sign = 1.0 if x_b > x_a else -1.0
    flip = [f for f in bm.faces if f.normal.x * sign > 0]
    if flip:
        bmesh.ops.reverse_faces(bm, faces=flip)
    res = bmesh.ops.extrude_face_region(bm, geom=bm.faces[:], use_keep_orig=True)
    bmesh.ops.translate(bm, verts=[g for g in res["geom"] if isinstance(g, bmesh.types.BMVert)], vec=(x_b - x_a, 0, 0))
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    uv = bm.loops.layers.uv.verify()
    for f in bm.faces:
        n = f.normal
        ax = max(range(3), key=lambda k: abs(n[k]))
        for l in f.loops:
            c = l.vert.co
            p = (c.y, c.z) if ax == 0 else (c.x, c.z) if ax == 1 else (c.x, c.y)
            l[uv].uv = (p[0] * 0.25, p[1] * 0.25)
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me)
    bm.free()
    if material is not None:
        me.materials.append(material)
    me["fit"] = ("closes the section of hull and roof pod between x %.2f and %.2f; edge %.0f mm inside the outer skin "
                 "(shipfit.section_wall_mesh)" % (xs[0], xs[2], inset * 1000.0))
    if info is not None:
        info.update({"row_lines": [round(z, 3) for z in zs], "half_widths": [round(val[z], 3) for z in zs], "columns": [round(y, 3) for y in cols],
                     "crown_half": yc, "crown_base_z": round(z_p, 3), "roof_edge_at_centre_z": round(apex, 3)})
    return me


def section_wall(name, coll, parent, env, x_a, x_b, zbot, **kw):
    """section_wall_mesh as a new object."""
    me = section_wall_mesh(name, env, x_a, x_b, zbot, **kw)
    ob = bpy.data.objects.new(name, me)
    coll.objects.link(ob)
    if parent is not None:
        ob.parent = parent
    return ob


def check_inside(ob, env, spacing=0.05, clear_min=0.0):
    """Every edge of a mesh object sampled every `spacing`: is each point inside the ship, and how near is the
    outer skin? (A vertex-only test misses a straight edge that cuts a hollow corner of the skin.)"""
    mw = ob.matrix_world
    me = ob.data
    seen, n, outside, near, cmin, worst = set(), 0, [], 0, 1e9, None
    for e in me.edges:
        a, b = mw @ me.vertices[e.vertices[0]].co, mw @ me.vertices[e.vertices[1]].co
        k = max(1, int(math.ceil((b - a).length / spacing)))
        for j in range(k + 1):
            p = a.lerp(b, j / k)
            key = (round(p.x, 4), round(p.y, 4), round(p.z, 4))
            if key in seen:
                continue
            seen.add(key)
            n += 1
            c = env.clearance(p)
            if not env.inside(p):
                outside.append((key, round(c, 4)))
                continue
            if c < cmin:
                cmin, worst = c, key
            if c < clear_min:
                near += 1
    return {"samples": n, "outside": len(outside), "furthest_outside_m": round(max([c for k, c in outside] or [0.0]), 4),
            "outside_sample": [k for k, c in sorted(outside, key=lambda r: -r[1])[:6]],
            "least_clearance_m": round(cmin, 4), "least_clearance_at": worst, "closer_than_%.0f_mm" % (clear_min * 1000): near}


def turned(name, coll, parent, profile, segments=24, materials=(), mat_of=None):
    """Lathe a (radius, z) profile about Z: quads on the wall, triangle fans on the two ends."""
    bm = bmesh.new()
    rings = []
    for r, z in profile:
        rings.append([bm.verts.new((r * math.cos(2 * math.pi * k / segments), r * math.sin(2 * math.pi * k / segments), z))
                      for k in range(segments)])
    uv = bm.loops.layers.uv.verify()
    lens = [0.0]
    for (r0, z0), (r1, z1) in zip(profile[:-1], profile[1:]):
        lens.append(lens[-1] + math.hypot(r1 - r0, z1 - z0))
    total = max(lens[-1], 1e-6)
    for i in range(len(rings) - 1):
        for k in range(segments):
            k2 = (k + 1) % segments
            f = bm.faces.new((rings[i][k], rings[i][k2], rings[i + 1][k2], rings[i + 1][k]))
            if mat_of:
                f.material_index = mat_of(i)
            us = (k / segments, (k + 1) / segments, (k + 1) / segments, k / segments)
            vs = (lens[i] / total, lens[i] / total, lens[i + 1] / total, lens[i + 1] / total)
            for l, u, v in zip(f.loops, us, vs):
                l[uv].uv = (u, v)
    for ring, z, flip in ((rings[0], profile[0][1], True), (rings[-1], profile[-1][1], False)):
        c = bm.verts.new((0, 0, z))
        for k in range(segments):
            k2 = (k + 1) % segments
            f = bm.faces.new((c, ring[k2], ring[k]) if flip else (c, ring[k], ring[k2]))
            for l in f.loops:
                l[uv].uv = (0.5 + l.vert.co.x * 0.2, 0.5 + l.vert.co.y * 0.2)
    bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
    me = bpy.data.meshes.new(name)
    bm.to_mesh(me)
    bm.free()
    for m in materials:
        me.materials.append(m)
    ob = bpy.data.objects.new(name, me)
    coll.objects.link(ob)
    if parent is not None:
        ob.parent = parent
    return ob


def audit(objs):
    """Counts for a list of mesh objects: faces, triangles, n-gons, open and non-manifold edges, zero-area faces."""
    tot = dict(meshes=0, faces=0, tris=0, ngons=0, open=0, nonmanifold=0, zero=0, no_uv=0)
    for ob in objs:
        bm = bmesh.new()
        bm.from_mesh(ob.data)
        tot["meshes"] += 1
        tot["faces"] += len(bm.faces)
        tot["tris"] += sum(1 for f in bm.faces if len(f.verts) == 3)
        tot["ngons"] += sum(1 for f in bm.faces if len(f.verts) > 4)
        tot["open"] += sum(1 for e in bm.edges if len(e.link_faces) == 1)
        tot["nonmanifold"] += sum(1 for e in bm.edges if len(e.link_faces) > 2)
        tot["zero"] += sum(1 for f in bm.faces if f.calc_area() < 1e-8)
        tot["no_uv"] += 0 if ob.data.uv_layers else 1
        bm.free()
    return tot