# 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: # Critique round 1 (M2): 72 columns (was 144) and stations by need (about 110, was 161): the Spearfish's and the Mako's
# density (64-80 columns). The pinned columns are the same features, fewer free columns between them.
M = 72
HALF = M // 2
J_LR = (10, 11, 12, 13)          # lower groove: shell edge, floor, floor, band edge
J_FBB = (16, 17)                 # the frame's lower edge: shell, frame
J_D0, J_SGB, J_D1, J_SGT = 18, 19, 22, 23
J_UR = (24, 25, 26, 27)          # upper groove: band (or frame) edge, floor, floor, shell edge
J_FG = 29                        # the front glass's edge
J_HC = 32                        # the collar's square (lat 0.43)
WALL = 0.008                    # the grooves' walls: 8 mm chamfers (round 1: 3 mm walls twisted where the grooves run diagonally)
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.012            # the groove's front wall over the frame: 45 mm deep, a 12 mm chamfer
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.030                 # the frame band's front end: a 30 mm chamfer (round 1: a 6 mm step
                                                          # folded its corner quads with the lower edge's)
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.20
    h = min(h, 0.035 + 0.165 * smoothstep(0.0, 0.50, Y_TAIL - y))      # the stern shoulder
    if x > 884.0:
        h = min(h, 0.030)                                                 # the raked nose face (the front glass's rows)
    elif x > 858.0:
        h = min(h, 0.055)
    elif x > 834.0:
        h = min(h, 0.045)
    if 806.0 <= x <= 834.0:
        h = min(h, 0.065)                                                 # the lower groove's descent onto the pod
    if 740.0 <= x <= 770.0:
        h = min(h, 0.070)                                                 # the pod's tail (inside its tail cone)
    for e in (DOOR_YA, DOOR_YB, SG_YA, SG_YB):                          # the openings' rounded corners: the grid hole's
        if abs(y - e) < 0.25:                                            # chords must stay within the 8 mm door gap
            h = min(h, 0.016 + 0.10 * smoothstep(0.0, 0.25, abs(y - e)))
    if abs(y - Y_CL) < A_SQ + 0.25:
        h = min(h, 0.098)
    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, DOOR_YA + 0.012, DOOR_YB - 0.012, SG_YA + 0.012, SG_YB - 0.012, 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:
        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"]
        # ride the NACELLE's own circle (its axis and drawn radius), not the grid's lobe: inside the nose piece the lobe sinks
        # into the hull and a point at a fixed angle round its sunk centre lands anywhere on the profile (round 1: the
        # columns jumped up the flank at x 845 and folded a whole station)
        rr_ = pod_r_true(min(ximg(y), 836.0))
        a = _nearest_arc(P, arc, (POD_L + rr_ * math.cos(th), POD_Z + rr_ * 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 + 0.030}            # the frame band's lower edge: a 30 mm chamfer (was a 3 mm 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))
            # at the nose's front (x 903 -> 912.8) the relaxing features all end at its lower corner and packed the free
            # columns 2 mm apart (slivers in the chin's cap): their pull fades out there, the free columns spread evenly
            w *= 1.0 - smoothstep(903.0, 911.5, ximg(y))
            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, set(ctrl)
    return solve(ctrl), P, arc, set(ctrl)


# Critique round 1 (M2, the wavy cheek): the free columns (those no feature pins at a station) are smoothed ALONG the hull,
# station to station, as fractions of their profile's length (a Gaussian over +-3 stations, sigma 6 cm); the pinned columns
# stay exact; then each station is made monotone again between its pins. The lobe's fillets and the chin's rising keel change
# a profile's length unevenly from one station to the next, which the per-station placement passed on as zig-zags.
SMOOTH_SIG = 0.06
_SM = {}


def smoothed_arcs(pod):
    if pod in _SM:
        return _SM[pod]
    raw = [column_arcs(i, pod) for i in range(NR)]
    Ls = [r[2][-1] for r in raw]
    F = [[a / L for a in r[0]] for r, L in zip(raw, Ls)]
    SF = []
    for i in range(NR):
        pins = raw[i][3]
        f = list(F[i])
        for j in range(1, HALF):
            if j in pins:
                continue
            acc = ws = 0.0
            for k in range(max(0, i - 3), min(NR, i + 4)):
                w = math.exp(-0.5 * ((YS[k] - YS[i]) / SMOOTH_SIG) ** 2)
                acc += w * F[k][j]
                ws += w
            f[j] = acc / ws
        SF.append(f)
    # the station pairs (steps 3-20 mm apart) keep the same free columns on both stations (a pinned one wins), so the step's
    # wall quads are not sheared
    for a_, b_ in ((I_BF1, I_BF0), (I_TG0b, I_TG0), (I_TG1, I_TG1a), (I_FE1, I_FE0)):
        pa, pb = raw[a_][3], raw[b_][3]
        for j in range(1, HALF):
            if j in pa and j not in pb:
                SF[b_][j] = SF[a_][j]
            elif j in pb and j not in pa:
                SF[a_][j] = SF[b_][j]
            elif j not in pa and j not in pb:
                SF[a_][j] = SF[b_][j] = 0.5 * (SF[a_][j] + SF[b_][j])
    out = []
    for i in range(NR):
        pins = raw[i][3]
        L = Ls[i]
        f = SF[i]
        gap = 0.002 / L
        for j in range(1, HALF + 1):                      # monotone, forward then backward, pins kept
            if j not in pins:
                f[j] = max(f[j], f[j - 1] + gap)
        for j in range(HALF - 1, 0, -1):
            if j not in pins:
                f[j] = min(f[j], f[j + 1] - gap)
        out.append([fj * L for fj in f])
    _SM[pod] = out
    return out


def ring2d(i):
    """station i's port columns as (lat, z)"""
    _, P, arc, _p = column_arcs(i)
    arcs = smoothed_arcs(True)[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")
# critique round 1: the light strips glow only on a thin line in the middle of each groove floor (M5): strip_v runs 0 -> 1 across
# a floor, strip_w is the floor's width (m); the skin's own UV islands (m11): UV_Grid in arc-length metres, uv_grid = 1 on its faces
lay_stv = bm.verts.layers.float.new("strip_v")
lay_stw = bm.verts.layers.float.new("strip_w")
lay_uvg = bm.faces.layers.int.new("uv_grid")
lay_uv = bm.loops.layers.uv.new("UV_Grid")
V = {}


def _strip_attr(i, j):
    """(v, w) of grid vertex (i, port column j): on a groove floor's two columns, else (-1, 0)"""
    for J in (J_LR, J_UR):
        if j in (J[1], J[2]):
            if J is J_LR and lower_groove_depth(i) < 0.5:
                return -1.0, 0.0
            if J is J_UR or band_zone(i) or YS[i] >= Y_TG1a - 1e-7:
                w = (P(i, J[2]) - P(i, J[1])).length
                return (0.0 if j == J[1] else 1.0), w
    if j >= J_UR[3] and i in (I_TG1a, I_TG0b):                # the groove across the nose top: its floor between two stations
        return (0.0 if i == I_TG1a else 1.0), abs(YS[I_TG0b] - YS[I_TG1a])
    return -1.0, 0.0


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[lay_stv], v[lay_stw] = _strip_attr(i, vj(k))
        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"


# the skin's UV chart: u = distance round each station's ring (levelled points, keel -> crest), v = one value per station, the
# running sum of the median step between consecutive stations (v per column folded the chart where columns run diagonally);
# u monotone along a ring and v monotone along the hull: every quad maps to a trapezoid, no folds. One island, port at +u and
# starboard at -u joined along the keel, cut at the crest.
UG = [[0.0] * (HALF + 1) for _ in range(NR)]
for i in range(NR):
    for j in range(1, HALF + 1):
        UG[i][j] = UG[i][j - 1] + (P(i, j) - P(i, j - 1)).length
VS = [0.0]
for i in range(1, NR):
    steps = sorted((P(i, j) - P(i - 1, j)).length for j in range(HALF + 1))
    VS.append(VS[-1] + steps[len(steps) // 2])
VG = [[VS[i]] * (HALF + 1) for i in range(NR)]
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])
        sgu = 1.0 if k < HALF else -1.0
        for l_, (a, b) in zip(f.loops, q):
            l_[lay_uv].uv = (sgu * UG[a][vj(b)], VG[a][vj(b)])
        f[lay_uvg] = 1
        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))
# the shroud is a STEPPED cylinder (the concept; critique round 1 m6): forward two thirds 14 mm under the drawn radius, a 45 deg
# step at two thirds of its length aft (x 525.4), then the lip ring at the drawn r 0.697 (SIDE's silhouette); the slot at x 530
# stays on the forward part
R_FWD = R_OUT - 0.014
Y_STEP = Y_TAIL + 0.66 * (Y_LIP - Y_TAIL)
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_FWD, Y_TAIL + 0.038, "shell"),
    (R_FWD, Y_STEP - 0.010, "shell"),
    (R_FWD + 0.003, Y_STEP - 0.006, "shell"),
    (R_OUT - 0.003, Y_STEP + 0.006, "frame"),
    (R_OUT, Y_STEP + 0.010, "frame"),
    (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([R_LAST[k] - 0.004 for k in range(M)])]      # hugs the stepped last ring (no fold)
    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)]
    if q_ in (0, 1):                                # the duct ring strip's floor (its line glows, M5)
        for v_ in cur:
            v_[lay_stv], v_[lay_stw] = float(q_), 0.030
    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
# (round 1) the flat front closed by a Coons grid over the ring (its four sides keel -> port mid -> crest -> starboard mid):
# the old ladder of strips between mirror pairs made 1-2 mm slivers where the section's top is nearly flat
Q4 = M // 8                                          # the patch's corners at the section's four shoulders (not at its flat
for f in coons_cap(bm, cur[Q4:] + cur[:Q4], MI["shell"]):   # keel and crest, where a corner of 180 deg folds the cells)
    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])


GLASS_T = 0.004


def glass_slab(bmg):
    """(critique round 2, M3) each pane a closed slab GLASS_T thick (inward), not one surface: a single refracting surface bent
    every view ray into the cabin toward the glass's normal (the concept camera saw the dashboard, not the port wall's screens),
    and from inside the cabin a one-sided pane is culled in Unreal"""
    bmg.normal_update()
    out = sum((f.normal * f.calc_area() for f in bmg.faces), Vector()).normalized()
    c0 = sum((v.co for v in bmg.verts), Vector()) / len(bmg.verts)
    res = bmesh.ops.solidify(bmg, geom=bmg.faces[:], thickness=GLASS_T)
    new = [v for v in res["geom"] if isinstance(v, bmesh.types.BMVert)]
    if new and (sum((v.co for v in new), Vector()) / len(new) - c0).dot(out) > 0:      # it went outward: the other way
        for v in new:
            v.co -= out * (2 * GLASS_T)
    bmg.normal_update()
    bmesh.ops.recalc_face_normals(bmg, faces=bmg.faces[:])
    bmg.normal_update()
    if sum((f.normal * f.calc_area()).dot(out) for f in bmg.faces if (f.calc_center_median() - c0).dot(out) > -GLASS_T * 0.5) < 0:
        bmesh.ops.reverse_faces(bmg, faces=bmg.faces[:])


def glass_side(poly, sg, step=0.05, rows=6):
    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_slab(bmg)
    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=10):
    n = 18
    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_slab(bmg)
    glass.merge(bmg, ACRYL)


glass_front()

# ------------------------------------------------------------------ swept trims of one constant section
def adaptive(poly, max_step=0.07, max_turn_deg=11.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.07, flange=0.056, max_turn=11.0):
    """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, max_turn)
    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)]
        # (round 2) shortened only where the outline turns OUTWARD (a concave stretch, where an outward flange would fold);
        # round a convex corner the outward flange fans out and keeps its length (a shortened one left gaps at the windows'
        # corners)
        def _concave(j):
            t0 = (P3[j] - P3[j - 1]).normalized()
            t1 = (P3[(j + 1) % n] - P3[j]).normalized()
            return (t1 - t0).dot(O3[j]) > 1e-4
        cc = [_concave(j) for j in range(n)]
        cc = [cc[j - 1] or cc[j] or cc[(j + 1) % n] for j in range(n)]
        ext = [sleeve(P3[j], N3[j], O3[j] * (min(1.0, 0.8 * rcs[j] / flange) if cc[j] else 1.0)) for j in range(n)]
        REPORT.setdefault("sweep_concave_points", []).append(sum(cc))
        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.016, 0.010), (-0.008, 0.011), (-0.016, 0.002), (0.0, -0.010)]   # 5 points (round 1: 7)


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)))


# critique round 2 (m2): the windows' inner flange reaches 0.16 m out under the lining (was 0.05): the lining's opening is up
# to two grid faces wider than the window, and through the 5 cm flange's gap the skin's underside was seen from outside
# through the glass (from inside the cabin, a see-through gap in Unreal). The lining's opening is fitted to the trims WITHOUT
# their flanges (FRAMES_TREE from a probe sweep), so the lining runs on over the flange and the reveal is closed.
WIN_FL = 0.16
# the cabin's front end wall (the lining's first ring, I_LF = I_FG1, sunk along the raked nose's normals): its aft-most point
Y_ENDWALL = max(YS[I_FG1] - NRM[I_FG1][k].y * liner_at(YS[I_FG1]) for k in range(M))


def _skin_along(q, n):
    """the base skin's point and normal under q along the opening's own normal n (a vertical or sideways ray picked the wrong
    face on the nose's steep sides)"""
    hit = SKIN.ray_cast(q + n * 0.25, -n, 0.6)
    if hit[0] is None:
        return None, None
    return hit[0], _snorm(hit[0], hit[2])


def win_sleeve(p, n, o, fl=WIN_FL, proj=None):
    """down to 22 mm past the lining, then the inner flange out over the lining's opening (which is up to two grid faces wider
    than the outline), its edge turned back toward the lining. (round 2, m2) The flange follows the lining's curve 16 mm on the
    cabin side of it (proj: the skin's point and normal under a point; the lining lies liner_at(y) under the skin), so a long
    flange neither leaves the lining nor runs through it"""
    d = LV_FRAME + liner_at(p.y) + 0.022
    S0, S1 = p - n * (GLASS_D + 0.004), p - n * d
    if fl <= 0.0:
        return [S0, S1]
    q_ = S1 + o * fl                                    # no further forward than 12 mm aft of the cabin's front end wall (the
    y_lim = Y_ENDWALL + 0.012                           # front glass's lower edge and corners: the flange would cross the wall)
    if q_.y < y_lim and S1.y > y_lim:
        fl = max(0.05, fl * (S1.y - y_lim) / max(1e-6, S1.y - q_.y))
    elif q_.y < y_lim:
        fl = 0.05
    k = max(1e-6, o.length)
    od = o / k
    for _ in range(8):
        out = [S0, S1]
        for t, extra in ((0.5 * fl, 0.016), (fl - 0.006, 0.016), (fl, 0.008)):
            Q, N = _skin_along(p + od * (t * k), n) if (proj is not None and fl > 0.06) else (None, None)   # short: straight
            if Q is None:
                out.append(S1 + o * t + (n * 0.008 if extra < 0.01 else Vector()))
            else:
                out.append(Q - N * (liner_at(Q.y) + extra))
        if fl <= 0.05 or min(q.y for q in out[2:]) >= y_lim or S1.y < y_lim:
            break
        fl = max(0.05, fl * 0.7)                         # the flange as it lies would cross the end wall: shorter
    return out


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,
                lambda p, n, o, sg=sg: win_sleeve(p, n, o, fl=WIN_FL - 0.06 * smoothstep(0.3, 0.8, -o.z / max(1e-6, o.length)),
                                                  proj=lambda q: on_side(sg, q.y, q.z)), flange=WIN_FL,
                max_turn=6.0)                       # (round 2, m3: the corners in 6 deg steps; 11 deg twisted the trim)
frame_sweep(frames, [(u, v) for u, v in FG_POLY], lambda u, v: on_top(u, v, LV_FRAME), WIN_PROF, FRAME,
            lambda p, n, o: win_sleeve(p, n, o, fl=0.05 + (WIN_FL - 0.05) * smoothstep(0.05, 0.25, p.z),
                                       proj=lambda q: on_top(q.x, q.y)), flange=WIN_FL, max_step=0.06, max_turn=6.0)
# (the front glass's flange is long only along its upper half, under the ceiling; at its lower corners the lining turns from
# the nose's end into the side walls within a few cm and a long flange ran through it)

# 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.06)

# ------------------------------------------------------------------ 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):
    _, P_, arc_, _p = column_arcs(i, False)
    arcs = smoothed_arcs(False)[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_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)
# the lining's opening keeps clear of the trims, their sleeves and the first 5 cm of their flanges (round 1's), not of the
# windows' long flanges beyond that (they run under the lining)
_probe = MB("_frames_probe", [FRAME])
for sg in (1, -1):
    frame_sweep(_probe, P_SG, lambda u, v, sg=sg: on_side(sg, u, v, LV_FRAME), WIN_PROF, FRAME, lambda p, n, o: win_sleeve(p, n, o, 0.05))
frame_sweep(_probe, [(u, v) for u, v in FG_POLY], lambda u, v: on_top(u, v, LV_FRAME), WIN_PROF, FRAME,
            lambda p, n, o: win_sleeve(p, n, o, 0.05), max_step=0.06)
for sg in (1, -1):
    frame_sweep(_probe, P_DOOR, lambda u, v, sg=sg: on_side(sg, u, v, LV_BAND), DOOR_PROF, FRAME, door_sleeve(sg), max_step=0.06)
FRAMES_TREE = BVHTree.FromBMesh(_probe.bm)
_probe.bm.free()


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 or min(math.hypot(p.x, p.y - Y_CL) for p in q) < 0.535):
        return True                                             # the collar's bore (r 0.5067) and its flange under the ceiling
    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)]
    q4_ = len(ring_) // 8
    coons_cap(bmc_, ring_[q4_:] + ring_[:q4_], 0)       # (round 1) a Coons grid, not strips between mirror pairs (slivers)
    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"