# Mako v04: the hull, its openings, glass, lining (exec'd by build_mako.py after mk_shape.py).
#
# Topology (the Spearfish's clean method, models/tender-spearfish/build/sf_hull.py):
#   * one quad grid: 80 columns round, placed by arc length on each station's profile (mk_shape.py), pinned to its features
#     (the jaw's top corner or the body's waist Q, the sill S or the bezel's foot, the lips' edges, the hood's edge, the
#     bezel's shoulder, the shell's edge E, the rear glass bay, the hatch) and relaxed back to their even places within
#     0.35 m past a feature's ends where no face (a station pair 6 mm apart) takes the jump;
#   * openings are regions of whole faces: the canopy (the nose, then the bubble above the jaw's rim, less the lips and the
#     hood), the rear glass bay, the hatch; round each a raised rim and a sleeve down to the lining;
#   * the bezel's front and back faces and the lower body's front face are faces between station pairs 6 mm apart;
#     the shell's plating stands 10 mm proud of the dark body on a column pair 3 mm apart (a real step, not paint);
#   * normals: the grid's own (central differences, one-sided at the faces), stored per vertex (an_x, an_y, an_z) with a face
#     flag (on_skin) so the finish pass gives the smooth skin custom normals; steps, chamfers and the bezel stay crisp.
from mathutils.bvhtree import BVHTree
from mathutils.interpolate import poly_3d_calc

exec(open(os.path.join(HERE, "mk_topo.py")).read())
M = 80
HALF = M // 2                    # port columns j = 0 (keel) .. 40 (crest), starboard mirrored: k = 80 - j
J_WB, J_Q, J_WT = 10, 13, 16     # the rear glass bay's edges round the waist; frames on 9..10 and 16..17
J_S = 16                         # the sill (the bubble's lower edge) or the bezel's foot
J_L = (18, 19, 24, 25)           # the lips: foot of the lower slope, its top, top of the upper slope, its foot
J_SU, J_SH = 22, 27              # the bezel: the corner under the shoulder, the shoulder's tip
J_E = 28                         # the shell's edge; the step on 28..29
J_HF, J_H = 31, 32               # the hood's edge on the bubble, and its frame
J_HH = 35                        # the hatch's edge
J_GU = 30                        # v04 round 2: the line between the canopy's pale upper glass and its clear lower glass
J_CL0, J_CL1, J_CA0, J_CA1 = 5, 6, 9, 10   # the jaw's rear "C" (SIDE): its lower arm (z -0.73..-0.61), its upper arm (-0.26..-0.14)
X_CA, X_CL = 368.0, 342.0        # where the C's upper and lower arms end
FR_W = 0.040                     # frame widths (arc on the skin)
RIM_W = 0.010                    # how far a frame's inner edge stands proud
GLASS_D, LINER = 0.022, 0.060
RELAX = 0.35

Y_CL = Yi(CL_IMG)                # the clamp axis
HATCH_R = 0.25
Y_HT = Yi(590.0)                 # v04: the hatch, flush, one more plate behind the dorsal spine (which now ends at x 552)
HY_A, HY_B = Y_HT - HATCH_R, Y_HT + HATCH_R
Y_BULK = Yi(362.0)               # cockpit's rear bulkhead, just behind the bezel
Y_BAY = Yi(866.0)                # the rear bay's end (the lining's cap), behind the glass bay
BAY4 = dict(poly=P_BAY)


def side_span(poly, y):
    zs = []
    n = len(poly)
    for i in range(n):
        (y1, z1), (y2, z2) = poly[i], poly[(i + 1) % n]
        if (y1 - y) * (y2 - y) <= 0 and y1 != y2:
            zs.append(z1 + (z2 - z1) * (y - y1) / (y2 - y1))
    return (min(zs), max(zs)) if len(zs) >= 2 else None


def opening_range(poly, hmin=0.06):
    ys = [p[0] for p in poly]
    a, b = min(ys), max(ys)
    ok = []
    for i in range(801):
        y = a + (b - a) * i / 800
        s = side_span(poly, y)
        if s and s[1] - s[0] >= hmin:
            ok.append(y)
    return min(ok), max(ok)


BAY4["ya"], BAY4["yb"] = opening_range(P_BAY)

# ------------------------------------------------------------------ stations
def spacing(y):
    x = ximg(y)
    h = 0.11
    h = min(h, 0.010 + 0.035 * smoothstep(0.0, 0.25, y - Y_NOSE))          # the blunt nose: fine rings at the tip
    if x < 345:
        h = min(h, 0.045)
    elif x < 470:
        h = min(h, 0.060)
    if BAY4["ya"] - 0.1 <= y <= BAY4["yb"] + 0.1:
        h = min(h, 0.065)
    h = min(h, 0.05 + 0.08 * smoothstep(0.0, 0.4, Y_TAIL - y))
    return h


Y0 = Yi(X_TIP + 1.0)                              # first ring, 7.5 mm behind the point
ys = [Y0]
while ys[-1] < Y_TAIL - 1e-6:
    ys.append(min(Y_TAIL, ys[-1] + spacing(ys[-1] + 0.5 * spacing(ys[-1]))))
special = [Y_JAW, Y_B0 + CH_B, Y_B1 - CH_B, Yi(420.0), Yi(432.0), Yi(445.0), HY_A, HY_B, Y_BULK, Y_BAY,
           BAY4["ya"] - FR_W, BAY4["ya"], BAY4["yb"], BAY4["yb"] + FR_W] + [q for q in (Y_HD0, Y_LP0) if q is not None]
for s in special:
    ys = [y for y in ys if abs(y - s) > 0.45 * spacing(s)] + [s]
for s in (Y_B0, Y_LB, Y_B1 - STEP):               # the faces: a station and one 6 mm after it, nothing in between
    ys = [y for y in ys if not (s - 0.012 < y < s + STEP + 0.012)] + [s, s + STEP]
ys = sorted(set(round(y, 6) for y in ys))
if Y_TAIL - ys[-1] > 1e-4:
    ys.append(Y_TAIL)
YS = ys
NR = len(YS)
def st(y):
    return min(range(NR), key=lambda i: abs(YS[i] - y))
I_JAW = st(Y_JAW)
I_HD0 = st(Y_HD0) if Y_HD0 is not None else None
I_LP0 = st(Y_LP0) if Y_LP0 is not None else None
I_B0, I_B0S, I_LB, I_LBS, I_B1S, I_B1 = st(Y_B0), st(Y_B0 + STEP), st(Y_LB), st(Y_LB + STEP), st(Y_B1 - STEP), st(Y_B1)
I_HA, I_HB = st(HY_A), st(HY_B)
I_BULK, I_BAY = st(Y_BULK), st(Y_BAY)
BAY4["ia"], BAY4["ib"] = st(BAY4["ya"]), st(BAY4["yb"])
BAY4["fa"], BAY4["fb"] = st(BAY4["ya"] - FR_W), st(BAY4["yb"] + FR_W)
assert I_B0S == I_B0 + 1 and I_LBS == I_LB + 1 and I_B1 == I_B1S + 1, (I_B0, I_B0S, I_LB, I_LBS, I_B1S, I_B1)
assert I_JAW < I_B0 < I_LB < I_B1S < I_BULK < I_HA < I_HB < I_BAY, (I_JAW, I_B0, I_LB, I_B1S, I_BULK, I_HA, I_HB, I_BAY)
STEP_BANDS = {I_B0, I_LB, I_B1S}                  # bands that are faces (the jump between a station pair)

# ------------------------------------------------------------------ what each station is made of (its parts are forced: the warp
# below moves its vertices along the hull, and a vertex of the bezel stays a bezel vertex wherever its arc puts it)
def st_case(i):
    if YS[i] < Y_JAW:
        return "nose"                                            # the bubble alone, ahead of the jaw
    return "canopy" if i <= I_B0 else ("bezel" if i <= I_B1S else "body")
def st_lower(i):
    return "jaw" if i <= I_LB else "lowerbody"
def st_inset(i):
    if st_case(i) != "bezel":
        return None
    d = min(YS[i] - Y_B0, Y_B1 - YS[i])                          # the 2.5 cm chamfers at both faces
    return CH_B * max(0.0, min(1.0, 1.0 - (d - STEP) / (CH_B - STEP)))

# ------------------------------------------------------------------ the collar's arcs (v04, SIDE first): a warp along the hull
# SIDE draws the collar as a band round the bubble's rounded back: its front edge at x 266 near the top, 300 at mid height,
# 281 at its foot, its back edge 310 / 340 / 288. Each vertex of the bezel's stations moves along the hull by the arc's offset
# at its height; the stations ahead (0.45 m) and behind (0.60 m) take the move up gradually (no ring overtakes another), and
# nothing moves below the jaw's top (z -0.13), where SIDE draws straight lines.
_BF = pchip(sorted(BAND_F)); _BB = pchip(sorted(BAND_B))
D_AHEAD, D_BEHIND = 0.45, 0.60
def warp(i, z):
    f = smoothstep(-0.13, -0.02, z)
    if f <= 0.0:
        return 0.0
    y = YS[i]
    dF, dB = (_BF(z) - X_BZ0) * S, (_BB(z) - X_BZ1) * S
    if i <= I_B0:
        return dF * f * max(0.0, 1.0 - (Y_B0 - y) / D_AHEAD)
    if i >= I_B1:
        return dB * f * max(0.0, 1.0 - (y - Y_B1) / D_BEHIND)
    t = max(0.0, min(1.0, (y - YS[I_B0S]) / (YS[I_B1S] - YS[I_B0S])))
    return ((1 - t) * dF + t * dB) * f

# ------------------------------------------------------------------ profiles and the columns on them
def arc_at_z(P, arc, z, lo_name="Q"):
    """arc position of height z on the profile's flank (z rises along body profiles)"""
    for i in range(1, len(P.p)):
        z0, z1 = P.p[i - 1][1], P.p[i][1]
        if (z0 - z) * (z1 - z) <= 0 and z0 != z1:
            return arc[i - 1] + (arc[i] - arc[i - 1]) * (z - z0) / (z1 - z0)
    return arc[-1] if z > P.p[-1][1] else 0.0


def arc_at_lat_top(P, arc, lat):
    """arc position, counted from the crest backward, where the profile's top reaches lateral offset lat"""
    for i in range(len(P.p) - 1, 0, -1):
        l0, l1 = P.p[i][0], P.p[i - 1][0]
        if (l0 - lat) * (l1 - lat) <= 0 and l0 != l1:
            return arc[i] + (arc[i - 1] - arc[i]) * (lat - l0) / (l1 - l0)
    return arc[-1]


_PROF = {}
def prof_at(y, case, lower, inset):
    key = (round(y, 5), case, lower, None if inset is None else round(inset, 5))
    if key not in _PROF:
        P = profile(y, case, lower, inset)
        _PROF[key] = (P, arc_of(P))
    return _PROF[key]


def frac_at(name, y, case, lower):
    """a feature's arc position as a fraction of the profile's length (None if the profile has no such feature)"""
    P, arc = prof_at(y, case, lower, None)
    return arc[P.f[name]] / arc[-1] if name in P.f else None


# the pins each kind of station takes: (feature, column, relax ahead of this image x)
PINS = {"nose": [("Q", J_Q, X_JAW4), ("S", J_S, X_JAW4), ("GU", J_GU, None)],
        "canopy": [("Q", J_Q, X_JAW4), ("S", J_S, X_JAW4), ("GU", J_GU, None)],
        "bezel": [("Q", J_Q, None), ("S", J_S, None), ("SU", J_SU, None), ("SH", J_SH, None)],
        "body": [("Q", J_Q, None), ("E", J_E, None)]}
if X_LIP0 is not None:
    PINS["canopy"] += [(n, j, X_LIP0) for n, j in zip(("L0a", "L0b", "L1a", "L1b"), J_L)]
if X_HOOD0 is not None:
    PINS["canopy"] += [("H", J_H, X_HOOD0)]


def columns_at(y, case, lower, inset):
    """the HALF + 1 column points (lat, z) on the profile at y made of the given parts"""
    P, arc = prof_at(y, case, lower, inset)
    L = arc[-1]
    ctrl = {0: 0.0, HALF: L}
    dflt = lambda j: L * j / HALF
    def pin(j, tgt, w):
        if w > 0:
            ctrl[j] = dflt(j) + (tgt - dflt(j)) * w
    for name, j, x_relax in PINS[case]:
        if name in P.f:
            pin(j, arc[P.f[name]], 1.0)
        elif x_relax is not None and y < Yi(x_relax):              # ahead of where it begins: relax to the even places
            w = 1.0 - smoothstep(0.0, RELAX, Yi(x_relax) - y)
            fr = frac_at(name, Yi(x_relax), "canopy" if case == "nose" else case, lower)
            if fr is not None and w > 0:
                pin(j, fr * L, w)
    if "E" in P.f:                                           # the shell's step: the next column 3 mm on
        pin(J_E + 1, arc[P.f["E"]] + 0.003, 1.0)
    if "H" in P.f:                                           # the hood's frame on the bubble
        pin(J_HF, arc[P.f["H"]] - FR_W, 1.0)
    if case == "body" and BAY_ON:
        w = 1.0 if BAY4["ya"] <= y <= BAY4["yb"] else 1.0 - smoothstep(0.0, RELAX, min(abs(y - BAY4["ya"]), abs(y - BAY4["yb"])) - FR_W)
        if w > 0:                                            # the rear glass bay
            yc = min(max(y, BAY4["ya"]), BAY4["yb"])
            zlo, zhi = side_span(BAY4["poly"], yc)
            aB, aT = arc_at_z(P, arc, zlo), arc_at_z(P, arc, zhi)
            for j, tgt in ((J_WB - 1, aB - FR_W), (J_WB, aB), (J_WT, aT), (J_WT + 1, aT + FR_W)):
                pin(j, tgt, w)
    if case == "body":
        w = 1.0 if HY_A <= y <= HY_B else 1.0 - smoothstep(0.0, RELAX, min(abs(y - HY_A), abs(y - HY_B)) - 0.04)
        if w > 0:                                            # the hatch (kept on the shell, above its edge)
            aE = arc[P.f["E"]]
            pin(J_HH - 1, max(aE + 0.15, arc_at_lat_top(P, arc, HATCH_R + 0.04)), w)
            pin(J_HH, max(aE + 0.19, arc_at_lat_top(P, arc, HATCH_R)), w)
    if lower == "lowerbody":                                 # the jaw's rear "C": its two arms on the lower body
        w = 1.0 if y <= Yi(X_CA) else 1.0 - smoothstep(0.0, RELAX, y - Yi(X_CA))
        if w > 0:
            zca0 = -0.26 + 0.09 * smoothstep(Yi(335.0), Yi(X_CA), y)        # the upper arm tapers to a point at its end
            for j, zz in ((J_CL0, -0.73), (J_CL1, -0.61), (J_CA0, zca0), (J_CA1, -0.14)):
                pin(j, arc_at_z(P, arc, zz), w)
    keys = sorted(ctrl)
    vals = [ctrl[k] for k in keys]
    for q in range(1, len(vals)):                            # keep them in order
        vals[q] = max(vals[q], vals[q - 1] + 0.002)
    return [at_arc(P, arc, interp_lin(keys, vals, j)) for j in range(HALF + 1)]


def interp_lin(xs, ys_, x):
    if x <= xs[0]:
        return ys_[0]
    if x >= xs[-1]:
        return ys_[-1]
    for (a, va), (b, vb) in zip(zip(xs, ys_), zip(xs[1:], ys_[1:])):
        if a <= x <= b:
            return va + (vb - va) * (x - a) / (b - a)
    return ys_[-1]


def ring(i):
    """station i's port half: (lat, y, z) per column, each column moved along the hull by the collar's warp at its height
    and taken from the profile there (two passes: the move depends on the height, the height on where the profile is cut)"""
    case, lower, inset = st_case(i), st_lower(i), st_inset(i)
    base = columns_at(YS[i], case, lower, inset)
    out = []
    for j in range(HALF + 1):
        lat, z = base[j]
        y = YS[i]
        for _ in range(2):
            d = warp(i, z)
            if abs(d) < 1e-5:
                break
            y = YS[i] + d
            lat, z = columns_at(y, case, lower, inset)[j]
        out.append((lat, y, z))
    return out


RING3 = [ring(i) for i in range(NR)]
COL2D = [[(p[0], p[2]) for p in r] for r in RING3]
REPORT["warp_stations"] = sum(1 for r in RING3 if any(abs(p[1] - r[0][1]) > 1e-4 for p in r))
def vj(k):
    """the port column of ring vertex k"""
    k %= M
    return k if k <= HALF else M - k
def fj(k):
    """the port face index of ring face k (faces k = 40..79 mirror 39..0)"""
    k %= M
    return k if k < HALF else M - 1 - k
BASE = [[Vector(((1 if k <= HALF else -1) * RING3[i][vj(k)][0], RING3[i][vj(k)][1], RING3[i][vj(k)][2])) for k in range(M)] for i in range(NR)]

# ------------------------------------------------------------------ normals from the grid (one-sided at the faces' station pairs)
def _nrm(i, k):
    a = i - 1 if (i > 0 and (i - 1) not in STEP_BANDS) else i
    b = i + 1 if (i < NR - 1 and i not in STEP_BANDS) else i
    if a == b:
        a, b = max(0, i - 1), min(NR - 1, i + 1)
    tl = BASE[b][k] - BASE[a][k]
    tc = BASE[i][(k + 1) % M] - BASE[i][(k - 1) % M]
    n = tl.cross(tc)
    j = vj(k)
    p2 = COL2D[i][j]
    q2 = COL2D[i][min(HALF, j + 1)]; o2 = COL2D[i][max(0, j - 1)]
    out = Vector((q2[1] - o2[1], 0.0, -(q2[0] - o2[0])))       # the profile's own outward normal (lat, z)
    if j in (0, HALF):
        out = Vector((0, 0, 1 if j == HALF else -1))
    if k > HALF:
        out.x = -out.x
    if n.length < 1e-12:
        return out.normalized()
    if n.dot(out) < 0:
        n = -n
    return n.normalized()
NRM = [[_nrm(i, k) for k in range(M)] for i in range(NR)]

# ------------------------------------------------------------------ roles of the grid faces (i, k): stations i..i+1, columns k..k+1
ROLE = {}
for i in range(NR - 1):
    ym = 0.5 * (YS[i] + YS[i + 1])
    for k in range(M):
        f = fj(k)
        if i < I_JAW:
            r = "glassup" if J_GU <= f < M - J_GU else "glass"
        elif i < I_B0:                                       # the canopy
            if f < J_Q:
                r = "jaw"
            elif f < J_S:
                r = "frame"                                  # the jaw's dark top rim: the glass's sill
            elif Y_LP0 is not None and ym >= Y_LP0 and J_L[0] <= f < J_L[3]:
                r = "lip"
            elif Y_HD0 is not None and ym >= Y_HD0 and f >= J_H:
                r = "hood"
            elif Y_HD0 is not None and ym >= Y_HD0 and f == J_HF:
                r = "frame"
            else:
                r = "glassup" if f >= J_GU else "glass"
        elif i < I_B1:                                       # the bezel, its two faces included
            if f >= J_S:
                r = "collar"
            elif f >= J_Q:
                r = "sill"
            else:
                r = "jaw" if ym < Y_LB else "lower"
        else:                                                # the body
            r = "shell" if f >= J_E else "lower"
        if r == "lower" and ((ym <= Yi(X_CA) and J_CA0 <= f < J_CA1) or (ym <= Yi(X_CL) and J_CL0 <= f < J_CL1)):
            r = "jaw"                                        # the C's arms in the jaw's paint
        if r == "lower" and i >= I_B1 and ym <= Yi(452.0) and J_Q <= f < J_E:
            r = "recess"                                     # v04 round 2: the lit cavity behind the collar, over the side body
        if i >= I_B1:
            if I_HA <= i < I_HB and f >= J_HH:
                r = "hatch"
            if BAY_ON and BAY4["ia"] <= i < BAY4["ib"] and J_WB <= f < J_WT:
                r = "window"
            elif BAY_ON and BAY4["fa"] <= i < BAY4["fb"] and ((J_WB <= f < J_WT) or f in (J_WB - 1, J_WT)):
                r = "frame"
        ROLE[(i, k)] = r
HOLES = ("glass", "glassup", "hatch", "window")
PAINT_ROLES = ("collar", "shell", "hood", "lip")


def base_raise(i, j):
    if i >= I_B1 and j >= J_E + 1:
        return H_SHELL
    return 0.0


def boundary_loops(kind):
    """closed loops of grid vertices round the faces of one hole kind (the canopy's two glasses are one opening)"""
    edges = {}
    kinds = ("glass", "glassup") if kind == "glass" else (kind,)
    for (i, k), r in ROLE.items():
        if r not in kinds:
            continue
        q = [(i, k % M), (i, (k + 1) % M), (i + 1, (k + 1) % M), (i + 1, k % M)]
        for a, b in zip(q, q[1:] + q[:1]):
            key = (a, b) if a < b else (b, a)
            edges[key] = edges.get(key, 0) + 1
    border = [e for e, n in edges.items() if n == 1]
    border = [e for e in border if not (e[0][0] == 0 and e[1][0] == 0)]   # the open front ring (the nose's glass fan) is no border
    nxt = {}
    for a, b in border:
        nxt.setdefault(a, []).append(b)
        nxt.setdefault(b, []).append(a)
    loops, seen = [], set()
    for start in nxt:
        if start in seen:
            continue
        loop, prev, cur = [start], None, start
        seen.add(start)
        while True:
            cand = [v for v in nxt[cur] if v != prev and (v not in seen or (v == start and len(loop) > 2))]
            if not cand:
                break
            v = cand[0]
            if v == start:
                break
            loop.append(v); seen.add(v); prev, cur = cur, v
        loops.append(loop)
    return loops

LOOPS = {kind: boundary_loops(kind) for kind in ("glass", "hatch", "window")}
RIM = {}
for kind, rim in (("glass", RIM_W), ("window", RIM_W)):      # the hatch is flush: no rim
    for loop in LOOPS[kind]:
        for ik in loop:
            RIM[ik] = rim

def P(ik):
    i, k = ik[0], ik[1] % M
    return BASE[i][k] + NRM[i][k] * (base_raise(i, vj(k)) + RIM.get((i, k), 0.0))

# ------------------------------------------------------------------ the hull mesh: skin, steps, sleeves, stern
hull = MB("SM_Mako_Hull", [PAINT, DARK, GRAPH, JAW, FIN])
MI_P, MI_D, MI_G, MI_J, MI_F = hull.mi(PAINT), hull.mi(DARK), hull.mi(GRAPH), hull.mi(JAW), hull.mi(FIN)
bm = bmesh.new()
lay_n = [bm.verts.layers.float.new(n) for n in ("an_x", "an_y", "an_z")]
lay_s = bm.faces.layers.int.new("on_skin")
V = {}
def hv(ik):
    i, k = ik[0], ik[1] % M
    if (i, k) not in V:
        v = bm.verts.new(P((i, k)))
        for c, l in zip(NRM[i][k], lay_n):
            v[l] = c
        V[(i, k)] = v
    return V[(i, k)]
MI_R = hull.mi(RECESS)
MATI = {"frame": MI_G, "jaw": MI_J, "sill": MI_F, "lower": MI_D, "recess": MI_R}
for (i, k), role in ROLE.items():
    if role in HOLES:
        continue
    f = bm.faces.new([hv(a) for a in ((i, k), (i, (k + 1) % M), (i + 1, (k + 1) % M), (i + 1, k))])
    f.material_index = MI_P if role in PAINT_ROLES else MATI[role]
    fjk = fj(k)
    crisp = (i in STEP_BANDS or role in ("frame", "collar", "sill") or (i >= I_B1 and fjk == J_E)
             or (role == "lip" and fjk in (J_L[0], J_L[2])) or (I_B0 <= i < I_B1)
             or any(a in RIM for a in ((i, k % M), (i, (k + 1) % M), (i + 1, (k + 1) % M), (i + 1, k % M))))
    f[lay_s] = 0 if crisp else 1
def sleeve_dir(kind, i, k):
    """the way a hole's sleeve goes in: the skin's normal; the hatch's without its along-the-hull part, so its fore and aft
    walls stand upright and stay parallel to the lid's as it swings about its side hinge"""
    n = NRM[i][k % M]
    if kind == "hatch":
        n = Vector((n.x, 0.0, n.z)).normalized()
    return n
for kind in ("glass", "hatch", "window"):                   # sleeves: from each hole's raised edge down to the lining
    for loop in LOOPS[kind]:
        rows = [[hv(ik) for ik in loop]]
        for d in (GLASS_D, LINER):
            rows.append([bm.verts.new(BASE[i][k % M] - sleeve_dir(kind, i, k) * d) for i, k in loop])
        ring_rows(bm, rows, closed=True, mat_index=MI_G)
iE = NR - 1                                                  # the stern: the last ring steps into a recessed panel
cenE = Vector((0, YS[iE], 0.5 * (COL2D[iE][0][1] + COL2D[iE][HALF][1])))
rE = [hv((iE, k)) for k in range(M)]
rows = [rE]
for sc_, dy in ((0.94, 0.0), (0.90, -0.03), (0.80, -0.03)):
    rows.append([bm.verts.new(cenE + (v.co - cenE) * sc_ + Vector((0, dy, 0))) for v in rE])
ring_rows(bm, rows[:3], closed=True, mat_index=MI_D)
ring_rows(bm, rows[2:], closed=True, mat_index=MI_G)
coons_cap(bm, rows[-1][::-1], MI_G)
bm.normal_update()
bmesh.ops.recalc_face_normals(bm, faces=bm.faces[:])
top_f = max(bm.faces, key=lambda f: f.calc_center_median().z)
if top_f.normal.z < 0:
    bmesh.ops.reverse_faces(bm, faces=bm.faces[:])
bm.normal_update()
hull.merge(bm, None)
REPORT["hull_stations"] = NR
REPORT["hull_columns"] = M
REPORT["boundary_loops"] = {k: [len(l) for l in v] for k, v in LOOPS.items()}

# ------------------------------------------------------------------ the skin as a surface for everything that sits on it
_SV, _SN, _SF = [], [], []
_SI = {}
for i in range(NR):
    for k in range(M):
        _SI[(i, k)] = len(_SV)
        _SV.append(P((i, k)))
        _SN.append(NRM[i][k])
for i in range(NR - 1):
    for k in range(M):
        _SF.append((_SI[(i, k)], _SI[(i, (k + 1) % M)], _SI[(i + 1, (k + 1) % M)], _SI[(i + 1, k)]))
_tip = len(_SV)
_SV.append(Vector((0, Y_NOSE, DZf(Y_NOSE)))); _SN.append(Vector((0, -1, 0)))
for k in range(M):
    _SF.append((_SI[(0, (k + 1) % M)], _SI[(0, k)], _tip))
_cap = len(_SV)
_SV.append(cenE + Vector((0, 0.0, 0))); _SN.append(Vector((0, 1, 0)))
for k in range(M):
    _SF.append((_SI[(iE, k)], _SI[(iE, (k + 1) % M)], _cap))
def _out(f):                                     # wind a face outward (the inside test reads the faces' own normals)
    a, b, c = (_SV[q] for q in f[:3])
    n = (b - a).cross(c - a)
    m = sum((_SN[q] for q in f), Vector())
    return f if n.dot(m) >= 0 else tuple(reversed(f))
_ng = (NR - 1) * M                               # the grid's quads share one winding: set it from a face on the mid-body's flank
_probe = _out(_SF[(NR // 2) * M + HALF // 2])
_flip = _probe != _SF[(NR // 2) * M + HALF // 2]
_SF = [tuple(reversed(f)) if _flip else f for f in _SF[:_ng]] + [_out(f) for f in _SF[_ng:]]
SKIN = BVHTree.FromPolygons(_SV, _SF)


def _snorm(p, fi):
    ids = _SF[fi]
    w = poly_3d_calc([_SV[q] for q in ids], p)
    n = Vector()
    for wi, q in zip(w, ids):
        n += _SN[q] * wi
    return n.normalized() if n.length > 1e-9 else Vector((0, 0, 1))


def normal_at(p):
    hit, nrm, fi, dist = SKIN.find_nearest(Vector(p))
    return _snorm(hit, fi) if hit is not None else Vector((0, 0, 1))


_DIRS = [Vector(d).normalized() for d in ((0, 0, 1), (0, 0, -1), (1, 0, 0.13), (-1, 0, 0.13), (0.3, 1, 0.2))]
def is_inside(p):
    """inside the skin: rays in five directions vote (the first face each meets faces away from the ray from inside); a
    nearest-point sign test fails beside the bezel's steps"""
    votes = 0
    for d in _DIRS:
        hit, nrm, fi, dist = SKIN.ray_cast(p, d)
        votes += -1 if hit is None else (1 if nrm.dot(d) > 0 else -1)
    return votes > 0


def inside_depth(p):
    """depth of p under the skin (negative: outside)"""
    p = Vector(p)
    if p.y < Y_NOSE or p.y > Y_TAIL + 0.005:
        return -1.0
    hit, nrm, fi, dist = SKIN.find_nearest(p)
    if hit is None:
        return -1.0
    return dist if is_inside(p) else -dist


def hull_z_at(x, y, top=True, inset=0.0):
    o = Vector((x, y, 4.0 if top else -4.0))
    hit, nrm, fi, dist = SKIN.ray_cast(o, Vector((0, 0, -1 if top else 1)))
    if hit is None:
        prof = profile(y).p
        return max(z for l, z in prof) if top else min(z for l, z in prof)
    n = _snorm(hit, fi)
    return hit.z + (-1 if top else 1) * inset / max(0.3, abs(n.z))


def hull_x_at(z, y, inset=0.0):
    hit, nrm, fi, dist = SKIN.ray_cast(Vector((4.0, y, z)), Vector((-1, 0, 0)))
    if hit is None:
        return 0.0
    n = _snorm(hit, fi)
    return hit.x - inset / max(0.3, abs(n.x))


def on_top(lat, y, lift=0.0):
    p = Vector((lat, y, hull_z_at(lat, y, True)))
    return p + normal_at(p) * lift if lift else p


def on_side(sg, y, z, lift=0.0):
    p = Vector((sg * hull_x_at(z, y), y, z))
    return p + normal_at(p) * lift if lift else p

# ------------------------------------------------------------------ glass: the canopy's and the bay's own grid, sunk to the sleeve
def glass_sink(i):
    """the nose is all glass: no sink at its point, the full 22 mm from 0.3 m back"""
    return GLASS_D * smoothstep(0.0, 0.30, YS[i] - Y_NOSE)
glass = MB("SM_Mako_Glass", [GLASS, GLASSUP, BAYGLASS])
bmg = bmesh.new()
G = {}
def gv(i, k):
    k %= M
    if (i, k) not in G:
        G[(i, k)] = bmg.verts.new(BASE[i][k] - NRM[i][k] * glass_sink(i))
    return G[(i, k)]
MG_C, MG_U, MG_B = glass.mi(GLASS), glass.mi(GLASSUP), glass.mi(BAYGLASS)
for (i, k), role in ROLE.items():
    if role in ("glass", "glassup", "window"):
        f_ = bmg.faces.new((gv(i, k), gv(i, k + 1), gv(i + 1, k + 1), gv(i + 1, k)))
        ym_ = 0.5 * (YS[i] + YS[i + 1])
        f_.material_index = MG_U if role == "glassup" else (MG_B if role == "window" and ym_ > Yi(BAY_FRONT_X) else MG_C)
tipg = bmg.verts.new((0, Y_NOSE, DZf(Y_NOSE)))
for k in range(M):                                 # the nose point: a fan (the bubble's tip)
    bmg.faces.new((gv(0, k + 1), gv(0, k), tipg)).material_index = MG_C
bmg.normal_update()
for f in bmg.faces:
    c = f.calc_center_median()
    n_out = normal_at(c)
    if f.normal.dot(n_out) < 0:
        f.normal_flip()
glass.merge(bmg, None)

# ------------------------------------------------------------------ cabin lining: the grid sunk 60 mm, from the jaw's front to the bay's end
cab = MB("SM_Mako_Cabin", [CAB, FLOOR, SEAT, SCREEN, CABGLOW, GRAPH, STEEL])
bml = bmesh.new()
lay_lin = bml.verts.layers.int.new("is_lining")      # build-time tag: the cabin check skips the lining itself
LV = {}
def lv(i, k):
    k %= M
    if (i, k) not in LV:
        LV[(i, k)] = bml.verts.new(BASE[i][k] - NRM[i][k] * LINER)
        LV[(i, k)][lay_lin] = 1
    return LV[(i, k)]
for i in range(I_JAW, I_BAY):
    for k in range(M):
        if ROLE[(i, k)] in HOLES:
            continue
        bml.faces.new((lv(i, k), lv(i + 1, k), lv(i + 1, (k + 1) % M), lv(i, (k + 1) % M)))
coons_cap(bml, [lv(I_BAY, k) for k in range(M)])
for v in bml.verts:
    v[lay_lin] = 1
bml.normal_update()
for f in bml.faces:                               # seen from inside: against the skin's normal
    c = f.calc_center_median()
    if f.normal.dot(normal_at(c)) > 0:
        f.normal_flip()
cab.merge(bml, CAB)

# the cockpit's rear bulkhead: the lining's section filled, both faces (the cockpit and the rear bay see it)
for d, sgn in ((0.0, -1), (0.03, 1)):
    bmb = bmesh.new()
    ring = [bmb.verts.new(BASE[I_BULK][k] - NRM[I_BULK][k] * LINER + Vector((0, d, 0))) for k in range(M)]
    coons_cap(bmb, ring)
    lb = bmb.verts.layers.int.new("is_lining")
    for v in bmb.verts:
        v[lb] = 1
    bmb.normal_update()
    if sum(f.normal.y for f in bmb.faces) * sgn < 0:
        bmesh.ops.reverse_faces(bmb, faces=bmb.faces[:])
    cab.merge(bmb, CAB)
REPORT["stage"] = "hull+glass+lining"