# Grouper: the hull grid, its openings, glass, frames and the cabin lining (exec'd by build_grouper.py after gr_shape.py).
#
# Topology (the Spearfish's and the Mako's clean method):
#   * ONE quad grid: M = 144 columns round (port j = 0 keel .. 72 crest, starboard mirrored k = 144 - j), placed by arc length
#     on each station's profile (gr_shape.py) and pinned to the measured lines: the waist band's two grooves (four columns
#     each: edge, 3 mm wall, floor, 3 mm wall, edge), the window frame's lower edge (a 3 mm pair: the step), the door's top and
#     bottom, the side glass's top and bottom, the front glass's edge, the collar's square; pins relax back to their even
#     places within 0.30 m past a feature's ends;
#   * stations spaced by need (8.5 cm amidships, 1.6 cm on the nose's raked face), with station PAIRS 3-6 mm apart where a
#     transverse step is: the band's front end / the frame's aft edge (x 839.6), the strip's groove across the nose top
#     (two walls), the frame's front end (x 911);
#   * LEVELS: every vertex is raised along the surface normal by its region's level: shell 0, the waist band -30 mm, the
#     grooves' floors -50 mm (beside the band) or -20 mm (beside the frame), the window frame band +25 mm. Faces between
#     vertices of different levels are the steps' walls (real faces, not paint);
#   * OPENINGS are regions of whole faces: the side windows and the front window (in the frame band), the side doors (in the
#     band), the collar's bore (a square of the grid mapped onto the circle r 0.512 in plan; the collar hides the warp);
#     their edges follow the smooth outlines station by station, and swept trims of one constant section (the Spearfish's
#     frame_sweep) cover them; the glass is its own quad ladder over the smooth outline, sunk 22 mm;
#   * the stern closes into the duct (the duct's rings continue the grid: the strip's groove round the hull, the shroud,
#     its rounded lip, the glowing inner wall, the stern plate filled by a Coons patch); the nose's chin closes with three
#     rings drawn in to a Coons patch;
#   * normals: the grid's own (central differences, one-sided at station pairs), 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 stay crisp.
from mathutils.bvhtree import BVHTree
from mathutils.interpolate import poly_3d_calc

M = 144
HALF = M // 2
J_LR = (18, 19, 20, 21)          # lower groove: shell edge, floor, floor, band edge
J_FBB = (32, 33)                 # the frame's lower edge: shell, frame
J_D0, J_SGB, J_D1, J_SGT = 35, 37, 43, 45
J_UR = (46, 47, 48, 49)          # upper groove: band (or frame) edge, floor, floor, shell edge
J_FG = 51                        # the front glass's edge
J_HC = 64                        # the collar's square (lat 0.43)
WALL = 0.003
RELAX = 0.30
LV_BAND, LV_FLOOR_B, LV_FLOOR_F, LV_FRAME = -0.030, -0.050, -0.020, 0.025
GLASS_D = 0.022
LINER = 0.113
HMIN = 0.03

# ------------------------------------------------------------------ the measured outlines (SIDE px unless said)
X_BF0 = 839.6                                   # the band's front wall = the window frame's aft edge
Y_BF0 = Yi(X_BF0); Y_BF1 = Y_BF0 - 0.006
X_TG0, X_TG1 = 887.2, 888.4                     # the strip's groove across the nose top (0.04 m aft of the glass)
Y_TG0 = Yi(X_TG0); Y_TG0b = Y_TG0 - 0.003
Y_TG1a = Yi(X_TG1); Y_TG1 = Y_TG1a - 0.003
X_FG0, X_FG1 = 890.2, 909.4                     # the front glass's top and bottom stations (FRONT_WINDOW on SIDE's nose)
Y_FG0, Y_FG1 = Yi(X_FG0), Yi(X_FG1)
X_FE0 = 911.0                                   # the frame band's front end at the chin
Y_FE0 = Yi(X_FE0); Y_FE1 = Y_FE0 - 0.006
X_BULK = 715.0                                  # the cabin's rear bulkhead (aft of the doors and of the hatch's bore)
Y_BULK = Yi(X_BULK)

# the frame band's outer edge (SIDE_WINDOW_FRAME_OUT, its aft corners squared to meet the upper groove and the band's front)
FRAME_TOP = [(839.5, 345.0), (845.0, 341.0), (870.0, 339.0), (887.0, 339.0), (893.0, 341.0)]
FRAME_BOT = [(839.5, 375.0), (841.0, 377.0), (848.0, 386.0), (860.0, 396.0), (880.0, 391.0), (905.0, 391.5), (911.5, 389.5)]


def lin_px(tab, x, k=1):
    """a table [(x, v0, v1 ...)] at image x, linear, clamped"""
    if x <= tab[0][0]:
        return tab[0][k]
    if x >= tab[-1][0]:
        return tab[-1][k]
    for a, b in zip(tab, tab[1:]):
        if a[0] <= x <= b[0]:
            t = (x - a[0]) / (b[0] - a[0])
            return a[k] + (b[k] - a[k]) * t
    return tab[-1][k]


# the measured lines as smooth curves (a lightly smoothed monotone cubic of image x): linear tables kink at every reading,
# and the columns pinned to them would zig-zag
def _sm(tab, k=1, sig=1.6):
    return pchip(smooth_tab([(r[0], r[k]) for r in tab], lambda x: sig))


_LR_U, _LR_L = _sm(WAIST_LOWER_RIM, 1), _sm(WAIST_LOWER_RIM, 2)
_UR_U, _UR_L = _sm(WAIST_UPPER_RIM, 1), _sm(WAIST_UPPER_RIM, 2)
_FT = _sm(FRAME_TOP, 1, 2.0)
_FB = _sm(FRAME_BOT, 1, 2.0)


def rim_lower(x):
    """(upper edge, lower edge) image y of the lower groove at image x"""
    return _LR_U(x), max(_LR_L(x), _LR_U(x) + 1.6)


def rim_upper(x):
    """(upper edge, lower edge) image y of the upper groove: the waist's upper rim, forward of the band along the frame's top"""
    if x < X_BF0:
        return _UR_U(x), max(_UR_L(x), _UR_U(x) + 1.6)
    le = _FT(x)
    return le - 3.0, le


def rounded_rect(x0, y0, x1, y1, r, n=6):
    pts = []
    for (cx, cy, a0) in ((x1 - r, y0 + r, -90), (x1 - r, y1 - r, 0), (x0 + r, y1 - r, 90), (x0 + r, y0 + r, 180)):
        for k in range(n + 1):
            a = math.radians(a0 + 90.0 * k / n)
            pts.append((cx + r * math.cos(a), cy + r * math.sin(a)))
    return pts


DOOR_PX = rounded_rect(DOOR["x"][0], DOOR["y"][0], DOOR["x"][1], DOOR["y"][1], DOOR["corner_r"])
P_DOOR = side_poly(DOOR_PX)                               # (Y, Z), the doorway's outline
P_SG = side_poly(smooth_closed(SIDE_WINDOW_GLASS, 2.5))     # (Y, Z), the side glass (its corners eased 2.5 px)


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


DOOR_YA, DOOR_YB = opening_range(P_DOOR)
SG_YA, SG_YB = opening_range(P_SG)

# the front glass: FRONT_WINDOW on SIDE's nose profile; half width 0.443 at its top (x 890.2) to 0.539 at its bottom
# (x 909.4), corners rounded 7 cm along the face. Chart: (lat, Y).
FG_R = 0.07


def _crest_s(y0, y1, n=60):
    """distance along the crest line between two stations"""
    s = 0.0
    for k in range(n):
        ya, yb = y0 + (y1 - y0) * k / n, y0 + (y1 - y0) * (k + 1) / n
        s += math.hypot(yb - ya, TOPf(yb) - TOPf(ya))
    return s


FG_LEN = _crest_s(Y_FG1, Y_FG0)


def fg_hw(y):
    """the front glass's half width at station y (its corners rounded)"""
    if y < Y_FG1 - 1e-6 or y > Y_FG0 + 1e-6:
        return None
    zc = TOPf(y)
    yimg = Z0 - zc / S
    t = max(0.0, min(1.0, (yimg - FRONT_WINDOW["y_top"]) / (FRONT_WINDOW["y_bot"] - FRONT_WINDOW["y_top"])))
    hw = Lp(FRONT_WINDOW["half_top_px"]) + (Lp(FRONT_WINDOW["half_bot_px"]) - Lp(FRONT_WINDOW["half_top_px"])) * t
    d = min(_crest_s(y, Y_FG0, 12), _crest_s(Y_FG1, y, 12))
    if d < FG_R:
        hw -= FG_R - math.sqrt(max(0.0, FG_R * FG_R - (FG_R - d) ** 2))
    return max(0.05, hw)


# the collar's square (plan): mapped onto the bore's circle
A_SQ = 0.43
R_BORE_HOLE = 0.512                       # the hull's and the lining's hole round the bore (the collar's wall is r 0.5067)

# ------------------------------------------------------------------ stations
def spacing(y):
    x = ximg(y)
    h = 0.10
    h = min(h, 0.026 + 0.074 * smoothstep(0.0, 0.30, Y_TAIL - y))       # the stern shoulder
    if x > 884.0:
        h = min(h, 0.019)                                                 # the raked nose face
    elif x > 858.0:
        h = min(h, 0.032)
    elif x > 834.0:
        h = min(h, 0.026)
    if 808.0 <= x <= 834.0:
        h = min(h, 0.030)                                                 # the lower groove's descent, the pod's shoulder
    if 836.0 <= x <= 850.0:
        h = min(h, 0.016)                                                 # the pod's rounded nose
    if 718.0 <= x <= 765.0:
        h = min(h, 0.034)                                                 # the pod's tail
    for e in (DOOR_YA, DOOR_YB, SG_YA, SG_YB):
        if abs(y - e) < 0.22:
            h = min(h, 0.020 + 0.05 * smoothstep(0.02, 0.22, abs(y - e)))
    if abs(y - Y_CL) < A_SQ + 0.25:
        h = min(h, 0.062)
    return h


ys = [Y_NOSE]
while ys[-1] < Y_TAIL - 1e-6:
    ys.append(min(Y_TAIL, ys[-1] + spacing(ys[-1] + 0.5 * spacing(ys[-1]))))
SPECIAL = [DOOR_YA, DOOR_YB, SG_YA, SG_YB, Y_FG0, Y_FG1, Y_CL - A_SQ, Y_CL + A_SQ, Y_BULK, Yi(830.0), Yi(826.0)]
for s in SPECIAL:
    ys = [y for y in ys if abs(y - s) > 0.45 * spacing(s)] + [s]
PAIRS = [(Y_BF1, Y_BF0), (Y_TG0b, Y_TG0), (Y_TG1, Y_TG1a), (Y_FE1, Y_FE0)]
for a, b in PAIRS:
    ys = [y for y in ys if not (a - 0.008 < y < b + 0.008)]
for a, b in PAIRS:
    ys += [a, b]
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_BF0, I_BF1 = st(Y_BF0), st(Y_BF1)
I_TG0, I_TG0b, I_TG1a, I_TG1 = st(Y_TG0), st(Y_TG0b), st(Y_TG1a), st(Y_TG1)
I_FE0, I_FE1 = st(Y_FE0), st(Y_FE1)
I_FG0, I_FG1 = st(Y_FG0), st(Y_FG1)
I_BULK = st(Y_BULK)
I_DA, I_DB = st(DOOR_YA), st(DOOR_YB)
I_SA, I_SB = st(SG_YA), st(SG_YB)
I_CA, I_CB = st(Y_CL - A_SQ), st(Y_CL + A_SQ)
for a, b in ((I_BF1, I_BF0), (I_TG0b, I_TG0), (I_TG1, I_TG1a), (I_FE1, I_FE0)):
    assert b == a + 1, ("pair", a, b)
REPORT["hull_stations"] = NR

# ------------------------------------------------------------------ features: target arcs per station
_PROF = {}


def prof_at(y, pod=True):
    key = (round(y, 6), pod)
    if key not in _PROF:
        P = profile(y, pod)
        _PROF[key] = (P, arc_of(P))
    return _PROF[key]


def has_pod(i):
    return True


def _clamp_lower_z(P, z):
    """the lower groove stays on the pod's lower outer side (where SIDE draws it plunging to the silhouette's bottom, x 826-830,
    the pod's bottom has no room for the columns of the pod's whole outer side above it)"""
    zmin = min(p[1] for p in P.p if p[0] > 0.45)
    zlim = zmin + 0.004
    if P.c is not None and P.r > 0.05:
        zlim = max(zlim, P.c[1] - 0.80 * P.r)
    return max(z, zlim)


_LR_REF = {}
X_LR_RIDE = 826.0


def _nearest_arc(P, arc, q):
    best = (9.0, 0.0)
    for i in range(1, len(P.p)):
        (a0, b0), (a1, b1) = P.p[i - 1], P.p[i]
        dx, dz = a1 - a0, b1 - b0
        L2 = dx * dx + dz * dz or 1e-12
        t = max(0.0, min(1.0, ((q[0] - a0) * dx + (q[1] - b0) * dz) / L2))
        d = math.hypot(a0 + dx * t - q[0], b0 + dz * t - q[1])
        if d < best[0]:
            best = (d, arc[i - 1] + (arc[i] - arc[i - 1]) * t)
    return best[1]


def f_lr(P, arc, y):
    """the lower groove: on the measured rim to x 826; from there, where SIDE draws it plunging to the pod's bottom, its
    columns ride on the pod's lobe at a fixed angle round the lobe's centre until the lobe ends (x 849), so they move with the
    receding nose instead of sliding along it (the groove itself runs out by x 830)"""
    x = ximg(y)
    if x > X_LR_RIDE and P.c is not None:
        if not _LR_REF:
            y0 = Yi(X_LR_RIDE)
            P0, arc0 = prof_at(y0)
            d0 = f_lr_rim(P0, arc0, y0)
            am = 0.5 * (d0[J_LR[0]] + d0[J_LR[3]])
            pm = at_arc(P0, arc0, am)
            _LR_REF["th"] = math.atan2(pm[1] - P0.c[1], pm[0] - P0.c[0])
            _LR_REF["w"] = d0[J_LR[3]] - d0[J_LR[0]]
        th, w = _LR_REF["th"], _LR_REF["w"]
        a = _nearest_arc(P, arc, (P.c[0] + P.r * math.cos(th), P.c[1] + P.r * math.sin(th)))
        return {J_LR[0]: a - 0.5 * w, J_LR[1]: a - 0.5 * w + WALL, J_LR[2]: a + 0.5 * w - WALL, J_LR[3]: a + 0.5 * w}
    return f_lr_rim(P, arc, y)


def f_lr_rim(P, arc, y):
    x = ximg(y)
    yu, yl = rim_lower(min(x, 830.0))
    z_lo, z_hi = _clamp_lower_z(P, Zi(yl)), _clamp_lower_z(P, Zi(yu))
    a0, a3 = arc_at_z_outer(P, arc, z_lo), arc_at_z_outer(P, arc, z_hi)
    if a3 - a0 < 0.02:
        m = 0.5 * (a0 + a3)
        a0, a3 = m - 0.01, m + 0.01
    return {J_LR[0]: a0, J_LR[1]: a0 + WALL, J_LR[2]: a3 - WALL, J_LR[3]: a3}


def f_ur(P, arc, y):
    yu, yl = rim_upper(ximg(y))
    a0, a3 = arc_at_z_outer(P, arc, Zi(yl)), arc_at_z_outer(P, arc, Zi(yu))
    if a3 - a0 < 0.02:
        m = 0.5 * (a0 + a3)
        a0, a3 = m - 0.01, m + 0.01
    return {J_UR[0]: a0, J_UR[1]: a0 + WALL, J_UR[2]: a3 - WALL, J_UR[3]: a3}


def f_fbb(P, arc, y):
    a = arc_at_z_outer(P, arc, Zi(_FB(ximg(y))))
    return {J_FBB[0]: a, J_FBB[1]: a + WALL}


def f_door(P, arc, y):
    sp = span_at(P_DOOR, min(max(y, DOOR_YA), DOOR_YB))
    return {J_D0: arc_at_z_outer(P, arc, sp[0]), J_D1: arc_at_z_outer(P, arc, sp[1])}


def f_sg(P, arc, y):
    sp = span_at(P_SG, min(max(y, SG_YA), SG_YB))
    return {J_SGB: arc_at_z_outer(P, arc, sp[0]), J_SGT: arc_at_z_outer(P, arc, sp[1])}


def f_fg(P, arc, y):
    hw = fg_hw(min(max(y, Y_FG1), Y_FG0))
    return {J_FG: arc_at_lat_top(P, arc, hw)}


def f_hc(P, arc, y):
    return {J_HC: arc_at_lat_top(P, arc, A_SQ)}


# (feature, valid y range [lo, hi], relax beyond). The frame's lower edge is also pinned for 0.27 m aft of the band's front
# wall (a band column there), so the free columns do not jump across the wall's station pair.
FEATURES = [(f_lr, (Yi(X_POD_FACE), Y_TAIL)), (f_ur, (Y_TG1a, Y_TAIL)), (f_fbb, (Y_FE0, Yi(824.0))), (f_door, (DOOR_YA, DOOR_YB)),
            (f_sg, (SG_YA, SG_YB)), (f_fg, (Y_FG1, Y_FG0)), (f_hc, (Y_CL - A_SQ, Y_CL + A_SQ))]
_FRAC = {}


def feature_fracs(fi, yb, pod):
    """a feature's columns at the end of its range, as distances along the profile: from the keel for the lower groove, from
    the crest for everything above it (the pod's lobe changes the profile's length below them)"""
    key = (fi, round(yb, 6))
    if key not in _FRAC:
        P, arc = prof_at(yb, pod)
        L = arc[-1]
        from_keel = FEATURES[fi][0] is f_lr
        _FRAC[key] = {j: (a if from_keel else a - L, from_keel) for j, a in FEATURES[fi][0](P, arc, yb).items()}
    return _FRAC[key]


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 column_arcs(i, pod=None):
    """the HALF + 1 column arc positions on station i's profile"""
    y = YS[i]
    if pod is None:
        pod = has_pod(i)
    P, arc = prof_at(y, pod)
    L = arc[-1]
    ctrl = {0: 0.0, HALF: L}
    soft = []                                         # features past their ends: (weight, {j: arc})
    for fi, (f, (lo, hi)) in enumerate(FEATURES):
        if f is f_lr and not pod:
            continue                                  # the lining (sections without the pods) needs no groove columns
        if lo - 1e-7 <= y <= hi + 1e-7:
            for j, a in f(P, arc, y).items():
                ctrl[j] = a
        else:
            yb = lo if y < lo else hi
            w = 1.0 - smoothstep(0.0, RELAX, abs(y - yb))
            if w > 0.0:
                soft.append((w, {j: min(max((d if fk else L + d), 0.0), L) for j, (d, fk) in feature_fracs(fi, yb, True).items()}))

    def solve(c):
        keys = sorted(c)
        vals = [c[k] for k in keys]
        for q in range(1, len(vals)):
            vals[q] = max(vals[q], vals[q - 1] + 0.002)
        return [interp_lin(keys, vals, j) for j in range(HALF + 1)]
    if soft:
        # a relaxing feature blends from its own place to where the other pins alone would put its columns (continuous: the
        # columns do not jump when its weight reaches zero)
        free = solve(ctrl)
        c2 = dict(ctrl)
        for w, tg in soft:
            for j, a in tg.items():
                if j not in ctrl:
                    c2[j] = free[j] + (a - free[j]) * w
        return solve(c2), P, arc
    return solve(ctrl), P, arc


def ring2d(i):
    """station i's port columns as (lat, z)"""
    arcs, P, arc = column_arcs(i)
    pts = [at_arc(P, arc, a) for a in arcs]
    pts[0] = (0.0, pts[0][1])
    pts[HALF] = (0.0, pts[HALF][1])
    return pts


COL2D = [ring2d(i) for i in range(NR)]


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 column of ring face k (faces k = 72..143 mirror 71..0)"""
    k %= M
    return k if k < HALF else M - 1 - k


BASE = [[Vector(((1 if k <= HALF else -1) * COL2D[i][vj(k)][0], YS[i], COL2D[i][vj(k)][1])) for k in range(M)] for i in range(NR)]

# ---- the collar's square: the grid's square round the clamp axis (columns J_HC .. M - J_HC, stations I_CA .. I_CB) mapped
# onto the circle r 0.512 in plan; the grid round it follows within 1.6 times the square (hidden under the collar's skirt,
# r 0.665, except a gentle pull on the diagonals); z back on the hull's own upper curve.
def z_top(x, y):
    zt, zb, zw, W, nu, nl = hull_params(y)
    q = min(1.0, abs(x) / W)
    return zw + (zt - zw) * (1.0 - q ** nu) ** (1.0 / nu)


A_SQ_Y = 0.5 * (YS[I_CB] - YS[I_CA])
Y_SQC = 0.5 * (YS[I_CB] + YS[I_CA])


def warp_collar(B):
    """map the collar's square of a grid onto the bore's circle in plan, the grid round it following within 1.6 times the
    square; z back on the hull's own upper curve"""
    for i in range(NR):
        if abs(YS[i] - Y_SQC) > 1.7 * A_SQ_Y + 0.05:
            continue
        for k in range(M):
            p = B[i][k]
            if p.z < 0.5:
                continue
            u, v = p.x / A_SQ, (p.y - Y_SQC) / A_SQ_Y
            rho = max(abs(u), abs(v))
            if rho < 1e-6 or rho > 1.6:
                continue
            dx, dy = p.x, p.y - Y_SQC
            rr = math.hypot(dx, dy)
            w = smoothstep(1.0, 1.6, rho)
            s_ = (1 - w) * R_BORE_HOLE * rho / rr + w
            nx, ny = dx * s_, Y_SQC + dy * s_
            B[i][k] = Vector((nx, ny, z_top(nx, ny)))


warp_collar(BASE)
REPORT["collar_square"] = {"stations": [I_CA, I_CB], "half_y": round(A_SQ_Y, 4), "centre_y": round(Y_SQC, 4)}

# ------------------------------------------------------------------ normals from the grid (one-sided across station pairs)
PAIR_BANDS = {I_BF1, I_TG0b, I_TG1, I_FE1}        # the band between i and i + 1 is a step


def _nrm(i, k, B=None, C2=None):
    B = BASE if B is None else B
    C2 = COL2D if C2 is None else C2
    a = i - 1 if (i > 0 and (i - 1) not in PAIR_BANDS) else i
    b = i + 1 if (i < NR - 1 and i not in PAIR_BANDS) else i
    if a == b:
        a, b = max(0, i - 1), min(NR - 1, i + 1)
    tl = B[b][k] - B[a][k]
    tc = B[i][(k + 1) % M] - B[i][(k - 1) % M]
    n = tl.cross(tc)
    j = vj(k)
    q2 = C2[i][min(HALF, j + 1)]; o2 = C2[i][max(0, j - 1)]
    out = Vector((q2[1] - o2[1], 0.0, -(q2[0] - o2[0])))
    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)]

# the skin's normals from the analytic surface (they become the custom normals): at each vertex the cross product of its own
# section's tangent and the tangent along the hull (the nearest points on the sections 8 mm ahead and behind). The grid's own
# differences go wrong where columns sweep diagonally across stations (the grooves' descent onto the pods, the pods' ends).
# Vertices of the collar's warped patch take the deck's height function's gradient.
import numpy as _np
_PA = {}


def _prof_np(y, pod=True):
    key = (round(y, 6), pod)
    if key not in _PA:
        _PA[key] = _np.array(profile(y, pod).p)
    return _PA[key]


def _nearest2(P2, q):
    a, b = P2[:-1], P2[1:]
    d = b - a
    L2 = (d * d).sum(1)
    L2[L2 < 1e-14] = 1e-14
    t = _np.clip(((q - a) * d).sum(1) / L2, 0.0, 1.0)
    c = a + d * t[:, None]
    i = int(_np.argmin(((c - q) ** 2).sum(1)))
    return c[i], d[i] / _np.sqrt(L2[i])


def analytic_normals(B, C2, pod=True, rows=None):
    out = [list(r) for r in NRM] if B is BASE else [list(r) for r in NRM_L]
    dy = 0.008
    for i in (range(NR) if rows is None else rows):
        y = YS[i]
        P0, Pp, Pm = _prof_np(y, pod), _prof_np(y + dy, pod), _prof_np(y - dy, pod)
        for k in range(M):
            p = B[i][k]
            if abs(p.y - y) > 1e-5:                              # the collar's warped patch: z = z_top(x, y)
                e = 0.004
                gx = (z_top(p.x + e, p.y) - z_top(p.x - e, p.y)) / (2 * e)
                gy = (z_top(p.x, p.y + e) - z_top(p.x, p.y - e)) / (2 * e)
                n = Vector((-gx, -gy, 1.0)).normalized()
            else:
                sg_ = 1.0 if (k <= HALF and p.x >= 0) or p.x > 1e-9 else -1.0
                q = _np.array((abs(p.x), p.z))
                c0, t0 = _nearest2(P0, q)
                cp, _ = _nearest2(Pp, q)
                cm, _ = _nearest2(Pm, q)
                tc = Vector((sg_ * t0[0], 0.0, t0[1]))
                ty = Vector((sg_ * (cp[0] - cm[0]), 2 * dy, cp[1] - cm[1]))
                n = tc.cross(ty)
                if n.length < 1e-12:
                    continue
                n.normalize()
            if n.dot(out[i][k]) < 0:
                n = -n
            out[i][k] = n
    return out


NRM = analytic_normals(BASE, COL2D)

# ------------------------------------------------------------------ levels and roles
def band_zone(i):
    return YS[i] >= Y_BF0 - 1e-7


def frame_zone(i):
    return Y_FE0 - 1e-7 <= YS[i] <= Y_BF1 + 1e-7


def lower_groove_depth(i):
    """the lower groove runs out on the pod's bottom (x 824 -> 830)"""
    return 1.0 - smoothstep(824.0, 830.0, ximg(YS[i]))


def level(i, j):
    if band_zone(i):
        if j <= J_LR[0]:
            return 0.0
        if j in (J_LR[1], J_LR[2]):
            d = lower_groove_depth(i)
            return LV_BAND + (LV_FLOOR_B - LV_BAND) * d
        if j <= J_UR[0]:
            return LV_BAND
        if j in (J_UR[1], J_UR[2]):
            return LV_FLOOR_B
        return 0.0
    if frame_zone(i):
        y = YS[i]
        if j <= J_FBB[0]:
            return 0.0
        if j <= J_UR[0]:
            return LV_FRAME
        if j in (J_UR[1], J_UR[2]):
            return LV_FLOOR_F if y >= Y_TG1a - 1e-7 else LV_FRAME
        # j >= J_UR[3]: the shell above the frame, the groove across the nose top, the frame over the front glass
        if y >= Y_TG0 - 1e-7:
            return 0.0
        if y >= Y_TG1a - 1e-7:
            return LV_FLOOR_F
        return LV_FRAME
    return 0.0                                                  # the chin, ahead of the frame's front end


LV = [[level(i, vj(k)) for k in range(M)] for i in range(NR)]


def hole_of(i, k):
    """the opening a face (i, k) belongs to, or None"""
    f = fj(k)
    ya, yb = YS[i], YS[i + 1]
    if I_DA <= i < I_DB and J_D0 <= f < J_D1:
        return "door_" + ("L" if k < HALF else "R")
    if I_SA <= i < I_SB and J_SGB <= f < J_SGT:
        return "sg_" + ("L" if k < HALF else "R")
    if I_FG1 <= i < I_FG0 and f >= J_FG:
        return "fg"
    if I_CA <= i < I_CB and f >= J_HC:
        return "collar"
    return None


HOLE = {(i, k): hole_of(i, k) for i in range(NR - 1) for k in range(M)}


def P(i, k):
    k %= M
    return BASE[i][k] + NRM[i][k] * LV[i][k]


# ------------------------------------------------------------------ seam fields for the surfaces pass (5 + signed distance, 0 = none)
def seam_chine(p):
    x = ximg(p.y)
    if not (545.0 <= x <= 632.0):
        return 0.0
    zc = Zi(lin_px(CHINE_SIDE, x))
    d = p.z - zc
    return 5.0 + d if abs(d) < 0.3 and abs(p.x) > 0.2 else 0.0


def seam_belly(p):
    d = p.y - Yi(BELLY_SEAM_X)
    return 5.0 + d if abs(d) < 0.3 and p.z < -0.9 else 0.0


# ------------------------------------------------------------------ the hull mesh
hull = MB("SM_Grouper_Hull", [GUN, BAND, FRAME, STRIP, DUCTG, RECESS])
MI = {"shell": hull.mi(GUN), "band": hull.mi(BAND), "frame": hull.mi(FRAME), "floor": hull.mi(STRIP), "duct": hull.mi(DUCTG),
      "recess": hull.mi(RECESS)}
bm = bmesh.new()
lay_n = [bm.verts.layers.float.new(n) for n in ("an_x", "an_y", "an_z")]
lay_sc = bm.verts.layers.float.new("seam_chine")
lay_sb = bm.verts.layers.float.new("seam_belly")
lay_s = bm.faces.layers.int.new("on_skin")
V = {}


def hv(i, k):
    k %= M
    if (i, k) not in V:
        p = P(i, k)
        v = bm.verts.new(p)
        for c, l in zip(NRM[i][k], lay_n):
            v[l] = c
        v[lay_sc] = seam_chine(p)
        v[lay_sb] = seam_belly(p)
        V[(i, k)] = v
    return V[(i, k)]


def lv_mat(l):
    if l >= LV_FRAME - 1e-6:
        return "frame"
    if l <= LV_FLOOR_B + 1e-6 or abs(l - LV_FLOOR_F) < 1e-6:
        return "floor"
    if abs(l - LV_BAND) < 1e-6:
        return "band"
    return "shell"


ROLE = {}
for i in range(NR - 1):
    for k in range(M):
        if HOLE[(i, k)]:
            continue
        q = [(i, k), (i, (k + 1) % M), (i + 1, (k + 1) % M), (i + 1, k)]
        lvs = [LV[a][b] for a, b in q]
        f = bm.faces.new([hv(a, b) for a, b in q])
        if max(lvs) - min(lvs) < 1e-6:
            role = lv_mat(lvs[0])
            f[lay_s] = 1
        else:                                                    # a step's wall: the higher side's material
            hi = max(lvs)
            role = "frame" if hi >= LV_FRAME - 1e-6 else ("band" if abs(hi - LV_BAND) < 1e-6 else "shell")
            if abs(hi - LV_BAND) < 1e-6 and min(lvs) < LV_BAND - 1e-6:
                role = "band"
            f[lay_s] = 0
        f.material_index = MI[role]
        ROLE[(i, k)] = role

# ---- the stern: the duct's rings continue the grid (the strip's groove round the hull, the shroud, its rounded lip, the
# glowing inner wall, the stern plate). Angles of the last ring's columns about the duct's axis.
iE = NR - 1
DUCT_C = Vector((0.0, 0.0, Zi(DUCT["y_c"])))
R_OUT = Lp(DUCT["r_out"])                                       # 0.697
Y_LIP = Yi(DUCT["x_aft"])                                       # 3.961
R_CORNER = Lp(DUCT["aft_corner_r"])
R_IN = Lp(DUCT["r_inner"]) * DUCT["r_out"] / 41.7               # the inner wall, BACK's ring scaled to SIDE's lip
R_THROAT = Lp(DUCT["r_throat"]) * DUCT["r_out"] / 41.7
Y_PLATE = Yi(541.0) + 0.06                                      # the stern plate inside the duct (behind the rotor)
ANG = []
for k in range(M):
    p = BASE[iE][k]
    ANG.append(math.atan2(p.z - DUCT_C.z, p.x))
DUCT_RINGS = [  # (radius, y, material of the band from the previous ring)
    (R_OUT - 0.025, Y_TAIL + 0.004, "shell"),
    (R_OUT - 0.025, Y_TAIL + 0.034, "floor"),                   # the duct ring strip (gr_measure LIGHT_STRIPS duct_ring)
    (R_OUT, Y_TAIL + 0.038, "shell"),
    (R_OUT, Y_LIP - R_CORNER, "shell")]
for q in range(1, 5):
    a = math.pi / 2 * q / 4
    DUCT_RINGS.append((R_OUT - R_CORNER + R_CORNER * math.cos(a), Y_LIP - R_CORNER + R_CORNER * math.sin(a), "shell"))
DUCT_RINGS += [(R_IN + 0.035, Y_LIP, "frame"), (R_IN + 0.008, Y_LIP - 0.010, "frame"), (R_IN, Y_LIP - 0.035, "frame"),
               (R_THROAT, Yi(527.0), "duct"), (R_THROAT + 0.004, Yi(536.0), "duct"), (R_IN - 0.02, Y_PLATE + 0.02, "duct"),
               (R_IN - 0.045, Y_PLATE, "duct"), (Lp(DUCT["r_hub"]) + 0.03, Y_PLATE, "frame")]
prev = [hv(iE, k) for k in range(M)]
R_LAST = [math.hypot(P(iE, k).x, P(iE, k).z - DUCT_C.z) for k in range(M)]
for q_, (r_, y_, role) in enumerate(DUCT_RINGS):
    if q_ == 0:                                     # the groove's inner wall stays inside every vertex of the hull's last ring
        cur = [bm.verts.new(Vector((rr * math.cos(ANG[k]), y_, DUCT_C.z + rr * math.sin(ANG[k]))))
               for k, rr in enumerate([min(r_, R_LAST[k] - 0.004) for k in range(M)])]
    else:
        cur = [bm.verts.new(Vector((r_ * math.cos(ANG[k]), y_, DUCT_C.z + r_ * math.sin(ANG[k])))) for k in range(M)]
    for k in range(M):
        f = bm.faces.new((prev[k], prev[(k + 1) % M], cur[(k + 1) % M], cur[k]))
        f.material_index = MI[role]
        f[lay_s] = 0
    prev = cur
for f in coons_cap(bm, prev, MI["frame"]):
    f[lay_s] = 0

# ---- the nose's chin: a chamfer ring drawn in 15 mm, then the flat front closed by a ladder of horizontal strips between the
# mirror pairs (port column j and starboard M - j), two triangles where the ladder meets the keel and the crest
i0 = 0
cen = sum((BASE[0][k] for k in range(M)), Vector()) / M
prev = [hv(0, k) for k in range(M)]
cur = []
for k in range(M):
    p = BASE[0][k]
    d = (cen - p)
    d.y = 0.0
    cur.append(bm.verts.new(p + d.normalized() * min(0.015, 0.4 * d.length) + Vector((0, -0.010, 0))))
for k in range(M):
    f = bm.faces.new((prev[k], prev[(k + 1) % M], cur[(k + 1) % M], cur[k]))
    f.material_index = MI["shell"]
    f[lay_s] = 0
for j in range(HALF):
    a, b = cur[j], cur[j + 1]
    c_, d_ = cur[(M - j - 1) % M], cur[(M - j) % M]
    vs_ = [a, b] + ([c_] if (M - j - 1) % M != j + 1 else []) + ([d_] if (M - j) % M != j else [])
    if len(set(vs_)) >= 3:
        f = bm.faces.new(list(dict.fromkeys(vs_)))
        f.material_index = MI["shell"]
        f[lay_s] = 0
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_columns"] = M

# ------------------------------------------------------------------ the base 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(BASE[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(cen + Vector((0, -0.02, 0))); _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(Vector((0.0, Y_TAIL + 0.05, DUCT_C.z))); _SN.append(Vector((0, 1, 0)))
for k in range(M):
    _SF.append((_SI[(iE, k)], _SI[(iE, (k + 1) % M)], _cap))


def _wind(f):
    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))


_SF = [_wind(f) for f in _SF]
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):
    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 base skin (negative: outside)"""
    p = Vector(p)
    if p.y < Y_NOSE - 0.03 or p.y > Y_TAIL + 0.04:
        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):
    hit, nrm, fi, dist = SKIN.ray_cast(Vector((x, y, 5.0 if top else -5.0)), Vector((0, 0, -1 if top else 1)))
    if hit is None:
        return None
    return hit.z


def on_top(lat, y, lift=0.0):
    hit, nrm, fi, dist = SKIN.ray_cast(Vector((lat, y, 5.0)), Vector((0, 0, -1)))
    if hit is None:
        return None, None
    n = _snorm(hit, fi)
    return hit + n * lift, n


def on_side(sg, y, z, lift=0.0):
    hit, nrm, fi, dist = SKIN.ray_cast(Vector((sg * 5.0, y, z)), Vector((-sg, 0, 0)))
    if hit is None:
        return None, None
    n = _snorm(hit, fi)
    return hit + n * lift, n


def hull_x_at(z, y):
    hit, nrm, fi, dist = SKIN.ray_cast(Vector((5.0, y, z)), Vector((-1, 0, 0)))
    return hit.x if hit is not None else 0.0


# ------------------------------------------------------------------ glass: each window's own ladder over its smooth outline, sunk 22 mm
glass = MB("SM_Grouper_Glass", [ACRYL])


def glass_side(poly, sg, step=0.03, rows=10):
    ys_ = [p[0] for p in poly]
    y0, y1 = min(ys_) + 0.002, max(ys_) - 0.002
    n = max(4, int(math.ceil((y1 - y0) / step)))
    cols = []
    for j in range(n + 1):
        t = 0.5 - 0.5 * math.cos(math.pi * j / n)
        y = y0 + (y1 - y0) * t
        sp = span_at(poly, y)
        if sp is None:
            continue
        lo, hi = sp
        col = []
        for r in range(rows + 1):
            p, nn = on_side(sg, y, lo + (hi - lo) * r / rows)
            col.append(p - nn * GLASS_D)
        cols.append(col)
    bmg = bmesh.new()
    VV = [[bmg.verts.new(p) for p in col] for col in cols]
    for j in range(len(VV) - 1):
        for r in range(rows):
            bmg.faces.new((VV[j][r], VV[j + 1][r], VV[j + 1][r + 1], VV[j][r + 1]))
    bmg.normal_update()
    if sum(f.normal.x * sg for f in bmg.faces) < 0:
        bmesh.ops.reverse_faces(bmg, faces=bmg.faces[:])
    glass.merge(bmg, ACRYL)


for sg in (1, -1):
    glass_side(P_SG, sg)


def fg_chart():
    """the front glass's outline in the chart (lat, Y), a closed loop"""
    n = 28
    yy = [Y_FG1 + (Y_FG0 - Y_FG1) * (0.5 - 0.5 * math.cos(math.pi * q / n)) for q in range(n + 1)]
    right = [(fg_hw(y), y) for y in yy]
    return right + [(-l, y) for l, y in reversed(right)]


FG_POLY = fg_chart()


def glass_front(nl=16):
    n = 28
    bmg = bmesh.new()
    rows = []
    for q in range(n + 1):
        y = Y_FG1 + 0.001 + (Y_FG0 - Y_FG1 - 0.002) * (0.5 - 0.5 * math.cos(math.pi * q / n))
        hw = fg_hw(y)
        row = []
        for c in range(nl + 1):
            p, nn = on_top(-hw + 2 * hw * c / nl, y)
            row.append(bmg.verts.new(p - nn * GLASS_D))
        rows.append(row)
    for q in range(n):
        for c in range(nl):
            bmg.faces.new((rows[q][c], rows[q][c + 1], rows[q + 1][c + 1], rows[q + 1][c]))
    bmg.normal_update()
    if sum(f.normal.y for f in bmg.faces) > 0:
        bmesh.ops.reverse_faces(bmg, faces=bmg.faces[:])
    glass.merge(bmg, ACRYL)


glass_front()

# ------------------------------------------------------------------ swept trims of one constant section
def adaptive(poly, max_step=0.05, max_turn_deg=6.0):
    pts = resample_closed(poly, 0.004)
    n = len(pts)
    keep = [0]
    acc_d = acc_t = 0.0
    for i in range(1, n):
        (y0, z0), (y1, z1), (y2, z2) = pts[i - 1], pts[i], pts[(i + 1) % n]
        a1 = math.atan2(z1 - z0, y1 - y0); a2 = math.atan2(z2 - z1, y2 - y1)
        acc_t += abs((a2 - a1 + math.pi) % (2 * math.pi) - math.pi)
        acc_d += math.hypot(y1 - y0, z1 - z0)
        if acc_d >= max_step or acc_t >= math.radians(max_turn_deg):
            keep.append(i); acc_d = acc_t = 0.0
    return [pts[i] for i in keep]


def frame_sweep(mb, poly, to3d, prof, mat, sleeve=None, max_step=0.05):
    """a trim along a closed chart outline: to3d(u, v) -> (point, normal) on its reference surface; prof [(a, h)] (a out of
    the opening along the surface, h out of it); sleeve(p, n, o) -> extra points continuing the section inward. The outward
    direction is worked out in 3D (perpendicular to the outline on the surface, away from its middle: a chart's axes stretch
    on the raked nose), and smoothed over the neighbours; a sleeve's flange is shortened where the outline turns tighter than
    its length (it would fold over itself)"""
    pts2 = adaptive(poly, max_step)
    n = len(pts2)
    P3, N3 = [], []
    for (u, v) in pts2:
        p, nrm = to3d(u, v)
        P3.append(p); N3.append(nrm)
    cen = sum(P3, Vector()) / n
    O3, RC = [], []
    for j in range(n):
        t = P3[(j + 1) % n] - P3[j - 1]
        o = t.cross(N3[j])
        o = (o - o.dot(N3[j]) * N3[j]).normalized()
        RC.append(t)
        O3.append(o)
    if sum(O3[j].dot(P3[j] - cen) for j in range(n)) < 0:
        O3 = [-o for o in O3]
    O3 = [(O3[j - 1] + 2 * O3[j] + O3[(j + 1) % n]).normalized() for j in range(n)]
    rc = []
    for j in range(n):
        t0 = (P3[j] - P3[j - 1]); t1 = (P3[(j + 1) % n] - P3[j])
        ang = t0.angle(t1) if t0.length > 1e-9 and t1.length > 1e-9 else 0.0
        rc.append(0.5 * (t0.length + t1.length) / ang if ang > 1e-6 else 9.0)
    loops = [[P3[j] + O3[j] * a + N3[j] * h for j in range(n)] for a, h in prof]
    if sleeve is not None:
        rcs = [min(rc[j - 1], rc[j], rc[(j + 1) % n]) for j in range(n)]
        ext = [sleeve(P3[j], N3[j], O3[j] * min(1.0, 0.8 * rcs[j] / 0.056)) for j in range(n)]
        for q in range(len(ext[0])):
            loops.append([ext[j][q] for j in range(n)])
    grid_rings(mb, loops, mat, expect=(1, n // 4, tuple(N3[n // 4])))
    return P3, N3, O3


frames = MB("_frames", [FRAME])
# the windows' inner lip: 16 mm over the glass, 12 mm proud of the frame band, its outer edge sunk 3 mm into the band,
# then a sleeve straight down to the lining
WIN_PROF = [(0.022, -0.003), (0.018, 0.008), (0.006, 0.012), (-0.010, 0.010), (-0.016, 0.004), (-0.016, -0.004),
            (0.0, -0.010)]


def liner_at(y):
    """the lining's depth under the skin: 113 mm (SECTION A-A's wall), thinning to 56 mm on the nose, where the hull curves
    faster than that (a parallel surface 113 mm in would fold there)"""
    return LINER * (1.0 - 0.5 * smoothstep(858.0, 885.0, ximg(y)))


def win_sleeve(p, n, o):
    """down to 10 mm inside the lining, then an inner flange 50 mm out over the lining's opening (which is one grid face
    wider than the outline), its edge turned back toward the lining"""
    d = LV_FRAME + liner_at(p.y) + 0.022
    return [p - n * (GLASS_D + 0.004), p - n * d, p - n * d + o * 0.050, p - n * (d - 0.008) + o * 0.056]


for sg in (1, -1):
    frame_sweep(frames, P_SG, lambda u, v, sg=sg: on_side(sg, u, v, LV_FRAME), WIN_PROF, FRAME, win_sleeve)
frame_sweep(frames, [(u, v) for u, v in FG_POLY], lambda u, v: on_top(u, v, LV_FRAME), WIN_PROF, FRAME, win_sleeve,
            max_step=0.04)

# the doorway's trim on the band: a flange and a rounded top round the opening, then a sleeve going in HORIZONTALLY (across the
# door's vertical hinge axis), so the leaf, which is the same horizontal extrusion, swings out past its jambs
DOOR_PROF = [(0.034, -0.003), (0.030, 0.006), (0.014, 0.010), (0.003, 0.008), (0.0, 0.002)]


def door_sleeve(sg):
    """horizontal, to 10 mm past the lining (the band is 30 mm in, the lining 113 mm under the skin), then an inner flange
    50 mm out over the lining's opening, parallel to the lining"""
    def f(p, n, o):
        sx = Vector((sg, 0, 0))
        d = (liner_at(p.y) - (-LV_BAND) + 0.022) / max(0.35, abs(n.x))
        q = p - sx * d
        o2 = (o - o.dot(n) * n).normalized()
        return [p - sx * 0.03, q, q + o2 * 0.050, q + o2 * 0.056 + n * 0.008]
    return f


DOOR_TRIM = {}
for sg in (1, -1):
    DOOR_TRIM[sg] = frame_sweep(frames, P_DOOR, lambda u, v, sg=sg: on_side(sg, u, v, LV_BAND), DOOR_PROF, FRAME, door_sleeve(sg),
                                max_step=0.04)

# ------------------------------------------------------------------ the cabin lining: the grid sunk 113 mm, from the rear bulkhead to the
# front glass's bottom (the openings left open: the trims' sleeves close them)
cab = MB("SM_Grouper_Cabin", [CABL, CABF, SEAT, SCREEN, CABD, FRAME])


def sink_of(i):
    zs = [COL2D[i][j][1] for j in range(HALF + 1)]
    return min(liner_at(YS[i]), 0.32 * (max(zs) - min(zs)))


I_LF = I_FG1
# the lining follows the hull's main body only: the same pinned columns on each station's section WITHOUT the pods (sunk into a
# pod's lobe the lining would fold through its thin nose and tail); the openings' columns land at the same heights
def ring2d_nopod(i):
    arcs, P_, arc_ = column_arcs(i, False)
    pts = [at_arc(P_, arc_, a) for a in arcs]
    pts[0] = (0.0, pts[0][1]); pts[HALF] = (0.0, pts[HALF][1])
    return pts


COL2D_L = [ring2d_nopod(i) if I_LF - 1 <= i <= I_BULK + 1 else COL2D[i] for i in range(NR)]
BASE_L = [[Vector(((1 if k <= HALF else -1) * COL2D_L[i][vj(k)][0], YS[i], COL2D_L[i][vj(k)][1])) for k in range(M)] for i in range(NR)]
warp_collar(BASE_L)
NRM_L = [[_nrm(i, k, BASE_L, COL2D_L) if I_LF - 1 <= i <= I_BULK + 1 else NRM[i][k] for k in range(M)] for i in range(NR)]

bml = bmesh.new()
lay_lin = bml.verts.layers.int.new("is_lining")
LVX = {}


def lv(i, k):
    k %= M
    if (i, k) not in LVX:
        v = bml.verts.new(BASE_L[i][k] - NRM_L[i][k] * sink_of(i))
        v[lay_lin] = 1
        LVX[(i, k)] = v
    return LVX[(i, k)]


P_DOOR_OUT = side_poly(rounded_rect(DOOR["x"][0] - 1.2, DOOR["y"][0] - 1.2, DOOR["x"][1] + 1.2, DOOR["y"][1] + 1.2, DOOR["corner_r"] + 1.2))
GLASS_TREE = BVHTree.FromBMesh(glass.bm)
FRAMES_TREE = BVHTree.FromBMesh(frames.bm)


def lining_hole(i, k):
    """the lining's openings: the windows' one grid face wider than the hull's (two on the nose, where the sleeves drift
    sideways along the fast-curving surface; they go in along the normal); the
    doorway's is the door's outline (+2 cm) in (y, z): its sleeve goes in horizontally, so that is its exact footprint; and
    no lining face within 2 cm of the glass"""
    if HOLE.get((i, k)) and not HOLE[(i, k)].startswith("door"):
        return True
    for di in (-2, -1, 0, 1, 2):
        for dk in (-2, -1, 0, 1, 2):
            h = HOLE.get((i + di, (k + dk) % M))
            if h and h != "collar" and not h.startswith("door") and (max(abs(di), abs(dk)) < 2 or ximg(YS[i]) > 868.0):
                return True
    q = [BASE_L[a][b % M] - NRM_L[a][b % M] * sink_of(a) for a, b in ((i, k), (i + 1, k), (i + 1, k + 1), (i, k + 1))]
    c = sum(q, Vector()) / 4
    if abs(c.x) > 0.6 and pip((c.y, c.z), P_DOOR_OUT):
        return True
    hit = GLASS_TREE.find_nearest(c)
    if hit[0] is not None and hit[3] < 0.02:
        return True
    if abs(c.x) > 0.25 or c.y < Y_FG0 + 0.1:                     # no lining face within 25 mm of a window or door trim
        hit = FRAMES_TREE.find_nearest(c)
        if hit[0] is not None and hit[3] < 0.025:
            return True
    if c.z > 0.8 and math.hypot(c.x, c.y - Y_CL) < 0.56:        # the collar's bore (r 0.5067) and its flange under the ceiling
        return True
    return False


# every second column and station (it is only seen through the windows and the hatch): a coarse face is left out where any
# of the four fine faces it covers is an opening
LST = list(range(I_LF, I_BULK + 1, 2))
if LST[-1] != I_BULK:
    LST.append(I_BULK)
for a_, b_ in zip(LST, LST[1:]):
    for k in range(0, M, 2):
        if any(lining_hole(i, kk) for i in range(a_, b_) for kk in (k, k + 1)):
            continue
        bml.faces.new((lv(a_, k), lv(b_, k), lv(b_, (k + 2) % M), lv(a_, (k + 2) % M)))
for v in bml.verts:
    v[lay_lin] = 1
bml.normal_update()
for f in bml.faces:                                  # seen from inside
    c = f.calc_center_median()
    if f.normal.dot(normal_at(c)) > 0:
        f.normal_flip()
cab.merge(bml, CABL)
for ii, sgn in ((I_BULK, -1), (I_LF, 1)):            # the rear bulkhead faces forward, the nose's end wall aft
    bmc_ = bmesh.new()
    lc_ = bmc_.verts.layers.int.new("is_lining")
    ring_ = [bmc_.verts.new(BASE_L[ii][k] - NRM_L[ii][k] * sink_of(ii)) for k in range(0, M, 2)]
    M2, H2 = M // 2, HALF // 2
    for j in range(H2):                              # horizontal strips between the mirror pairs (triangles at keel and crest)
        vs_ = [ring_[j], ring_[j + 1], ring_[(M2 - j - 1) % M2], ring_[(M2 - j) % M2]]
        vs_ = list(dict.fromkeys(vs_))
        if len(vs_) >= 3:
            bmc_.faces.new(vs_)
    near_trim = [f for f in bmc_.faces if (lambda h: h[0] is not None and h[3] < 0.025)(FRAMES_TREE.find_nearest(f.calc_center_median()))]
    if near_trim:                                    # the window trims' flanges close the end wall there
        bmesh.ops.delete(bmc_, geom=near_trim, context='FACES')
    for v in bmc_.verts:
        v[lc_] = 1
    bmc_.normal_update()
    if sum(f.normal.y for f in bmc_.faces) * sgn < 0:
        bmesh.ops.reverse_faces(bmc_, faces=bmc_.faces[:])
    cab.merge(bmc_, CABL)
REPORT["stage"] = "hull+glass+frames+lining"