# Grouper, round 3 (2026-10-22): the hull's SKIN, rebuilt after the user's "fix the topology and lines needs to be smooth"
# (exec'd by gr_hull.py after gr_shape.py; shares the build's namespace).
#
# What round 2 had (and why it went): one station grid (planar rings x 72 columns) whose columns were pinned to coarse, noisy
# readings of the sheet's lines (the grooves bent in straight pieces between stations 6-20 cm apart); the transverse steps
# (the band's front wall, the groove across the nose top, the frame's front end, the openings' corners) were station PAIRS
# 3-12 mm apart that ran as thin rings round the whole hull; the columns met the nose's face and bunched; the groove turned
# across the nose top at a square corner.
#
# Round 3:
#   * every design line is a FAIR curve (gr_fair.py: tracked on the sheet to 0.05 px, a smoothing spline with a penalty on
#     the curvature's variation; checks/lines_fit.json);
#   * the skin is a STRUCTURED quad layout of a rounded box:
#       CYL  the cylinder: rings round the hull from the duct (x 543) forward; planar stations placed by need (where a line
#            or a silhouette curves, graded <= 1.25 x from one interval to the next) up to the band's front wall, then
#            transition rings that lean forward into C;
#       F    the nose cap, a Coons patch forward of the ring C: the flank's columns (the groove, the window rows, FRAME_BOT)
#            wrap round the nose as its u-lines (so the groove crosses the nose top as the SAME four loops, through a round
#            corner), the deck's middle columns come down the face as its v-lines and run aft under the chin as the belly's
#            columns. Its four corners are the skin's only poles (valence 3): two on the deck at the middle of the groove's
#            corners (with 8 rows of the deck wrapped round each, so the groove turns over many vertices), two on the
#            underside behind the chin;
#   * openings are HOLE RINGS: a rectangle of the grid is taken out and a ring of quads joins its edge to the opening's
#     outline (offset under the trim's lip), so no extra station or column is needed at an opening's rounded corners;
#   * no station pairs: a step (the band's front wall, the frame's chamfers) is one ordinary interval of the grid;
#   * the free vertices are RELAXED on the analytic surface (tangential Laplacian, projected back each pass; the design
#     lines' vertices slide along their curves);
#   * the analytic surface is sampled densely (DY x ND) for the projection and for every query of the later stages
#     (on_side, on_top, hull_z_at, inside_depth ...), independent of the skin's own topology.
import numpy as np
import json as _json
from mathutils.bvhtree import BVHTree
from mathutils.interpolate import poly_3d_calc

SKIN_T0 = time.time()
WALL = 0.008                    # the grooves' walls: 8 mm chamfers
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
X_BF0 = 839.6                   # the band's front wall (the window frame's aft edge)
Y_BF0 = Yi(X_BF0)
X_BULK = 715.0                  # the cabin's rear bulkhead
Y_BULK = Yi(X_BULK)

# ================================================================== 1. the analytic surface, densely sampled
ND = 300                        # points per port half-section (keel -> crest)
DY = 0.006                      # station step of the dense sampling (m)


def _resample(pts, n):
    P = np.asarray(pts, float)
    seg = np.hypot(*np.diff(P, axis=0).T)
    s = np.concatenate([[0.0], np.cumsum(seg)])
    t = np.linspace(0.0, s[-1], n + 1)
    return np.stack([np.interp(t, s, P[:, 0]), np.interp(t, s, P[:, 1])], 1)


_YD = np.arange(Y_NOSE, Y_TAIL + 1e-9, DY)
if _YD[-1] < Y_TAIL - 1e-6:
    _YD = np.append(_YD, Y_TAIL)
_YD[-1] = Y_TAIL
DENSE = np.stack([_resample(profile(float(y)).p, ND) for y in _YD])          # (n, ND+1, 2): lat, z
DENSE[:, 0, 0] = 0.0
DENSE[:, -1, 0] = 0.0
_nd = len(_YD)
_ring = 2 * ND                                                           # port 0..ND, starboard ND+1..2ND-1 (mirror)
_DV = np.zeros((_nd, _ring, 3))
_DV[:, :ND + 1, 0] = DENSE[:, :, 0]
_DV[:, :ND + 1, 2] = DENSE[:, :, 1]
_DV[:, ND + 1:, 0] = -DENSE[:, ND - 1:0:-1, 0]
_DV[:, ND + 1:, 2] = DENSE[:, ND - 1:0:-1, 1]
_DV[:, :, 1] = _YD[:, None]
_dverts = [Vector(p) for p in _DV.reshape(-1, 3)]
_dfaces = []
for i in range(_nd - 1):
    a, b = i * _ring, (i + 1) * _ring
    for k in range(_ring):
        k2 = (k + 1) % _ring
        _dfaces.append((a + k, a + k2, b + k2, b + k))
# the nose's flat front (x 913): a planar Coons grid over the first ring, its corners at the section's shoulders
_nq = _ring // 4
_q0 = ND // 4
_cring = [((k + _q0) % _ring) for k in range(_ring)]
_cP = [_DV[0, k] for k in _cring]
_cG = {}
for ii in range(1, _nq):
    for jj in range(1, _nq):
        u, v = ii / _nq, jj / _nq
        B_, T_ = _cP[ii], _cP[(3 * _nq - ii) % _ring]
        L_, R_ = _cP[(4 * _nq - jj) % _ring], _cP[_nq + jj]
        p = ((1 - v) * B_ + v * T_ + (1 - u) * L_ + u * R_
             - ((1 - u) * (1 - v) * _cP[0] + u * (1 - v) * _cP[_nq] + (1 - u) * v * _cP[3 * _nq] + u * v * _cP[2 * _nq]))
        _cG[(ii, jj)] = len(_dverts)
        _dverts.append(Vector(p))


def _cidx(ii, jj):
    if jj == 0:
        return _cring[ii]
    if ii == _nq:
        return _cring[_nq + jj]
    if jj == _nq:
        return _cring[(3 * _nq - ii) % _ring]
    if ii == 0:
        return _cring[(4 * _nq - jj) % _ring]
    return _cG[(ii, jj)]


_n_side_faces = len(_dfaces)
for ii in range(_nq):
    for jj in range(_nq):
        _dfaces.append((_cidx(ii, jj), _cidx(ii + 1, jj), _cidx(ii + 1, jj + 1), _cidx(ii, jj + 1)))
# orientation: the side faces share one winding (turned so the top-most faces up), the cap's faces forward (-Y)
_DVa = np.array([tuple(v) for v in _dverts])
_F4 = np.array(_dfaces)
_fn = np.cross(_DVa[_F4[:, 2]] - _DVa[_F4[:, 0]], _DVa[_F4[:, 3]] - _DVa[_F4[:, 1]])
_side = np.arange(len(_F4)) < _n_side_faces
_top = int(np.argmax(np.where(_side, _DVa[_F4].mean(1)[:, 2], -9.0)))
if _fn[_top, 2] < 0:
    _fn[_side] *= -1
    _F4[_side] = _F4[_side][:, ::-1]
_capf = ~_side
if (_fn[_capf, 1] > 0).sum() > (_capf.sum() // 2):
    _fn[_capf] *= -1
    _F4[_capf] = _F4[_capf][:, ::-1]
_DN = np.zeros_like(_DVa)
for c in range(4):
    np.add.at(_DN, _F4[:, c], _fn)
_DN /= np.maximum(1e-12, np.linalg.norm(_DN, axis=1))[:, None]
# (round 3) the side's normals smoothed over the dense grid: the sections are polylines (the fillets' arcs in 8 mm chords, 9
# degrees apart on a 5 cm ball) whose vertices shift from section to section as the fillets grow, so the face normals carried
# a few degrees of noise that changed with y: the skin, sampling it at its own stations, showed zebra stripes in steps. Two
# passes of a [1 2 1] x [1 4 6 4 1] kernel (about 6 mm along the hull, 2 cm round the section) take it out; the smallest
# designed radius (the pod's fillets, 4 cm and up) keeps its normals.
_ns = _DN[:_nd * _ring].reshape(_nd, _ring, 3).copy()
for _pass in range(2):
    _k = (np.roll(_ns, 2, 1) + 4 * np.roll(_ns, 1, 1) + 6 * _ns + 4 * np.roll(_ns, -1, 1) + np.roll(_ns, -2, 1)) / 16.0
    _y = _k.copy()
    _y[1:-1] = (_k[:-2] + 2 * _k[1:-1] + _k[2:]) / 4.0
    _y /= np.maximum(1e-12, np.linalg.norm(_y, axis=2))[:, :, None]
    _ns = _y
_ns[:, 0, 0] = 0.0                                    # the keel's and the crest's normals stay in the centre plane
_ns[:, ND, 0] = 0.0
_ns /= np.maximum(1e-12, np.linalg.norm(_ns, axis=2))[:, :, None]
_DN[:_nd * _ring] = _ns.reshape(-1, 3)
DSKIN = BVHTree.FromPolygons(_dverts, [tuple(f) for f in _F4])
_DNv = [Vector(n) for n in _DN]
_F4l = [tuple(f) for f in _F4]


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


def proj(p):
    """the nearest point of the analytic surface and its normal"""
    hit, nrm, fi, dist = DSKIN.find_nearest(Vector(p))
    return hit, snorm_dense(hit, fi)


def normal_at(p):
    return proj(p)[1]


# the base-surface queries the later stages use (round 2's names): the dense surface instead of the skin's own grid
SKIN = DSKIN
_SV = _dverts


def _snorm(p, fi):
    return snorm_dense(p, fi)


_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)
        if hit is not None and hit.y > Y_TAIL - 1e-4 and d.y > 0:
            hit = None                                   # the open stern (the duct is not part of this surface)
        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): along the normal near the surface, by ray votes far from it"""
    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
    n = snorm_dense(hit, fi)
    d = (p - hit).dot(n)
    if dist > 0.25:                                    # far from the surface: vote with rays (deep inside or far out)
        return dist if is_inside(p) else -dist
    return -d


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


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)


REPORT["skin_dense"] = {"stations": int(_nd), "points_per_half": ND, "seconds": round(time.time() - SKIN_T0, 2)}

# ================================================================== 2. fair curves (gr_fair.py) and helpers on sections
FAIR = _json.load(open(os.path.join(HERE, "fair_curves.json")))


def fair_xy(name):
    c = FAIR[name]
    return np.stack([np.array(c["x"]), np.array(c["y"])], 1)


def fair_fn(name):
    """a fair curve that is a function of the sheet's x: y(x) (clamped)"""
    P = fair_xy(name)
    o = np.argsort(P[:, 0])
    xs, ys = P[o, 0], P[o, 1]
    return lambda x: float(np.interp(x, xs, ys))


UR_U, UR_L = fair_fn("ur_u"), fair_fn("ur_l")       # the upper groove: its upper (shell) and lower (band / frame) edges
LR_U, LR_L = fair_fn("lr_u"), fair_fn("lr_l")       # the lower groove: its upper (band) and lower (shell) edges
X_UR_END = float(fair_xy("ur_u")[:, 0].max())
X_LR_END = float(fair_xy("lr_u")[:, 0].max())


def frame_bottom_curve():
    """FRAME_BOT: the lower part of the window frame's fair outline, from its aft edge's foot to the nose's front (x, y px)"""
    P = fair_xy("frame_out")
    n = len(P)
    i0 = int(np.argmin(np.hypot(P[:, 0] - 839.6, P[:, 1] - 366.0)))
    i1 = int(np.argmin(np.hypot(P[:, 0] - 911.5, P[:, 1] - 389.5)))
    # the way round that passes the lowest point
    path_a = [P[(i0 + k) % n] for k in range((i1 - i0) % n + 1)]
    path_b = [P[(i0 - k) % n] for k in range((i0 - i1) % n + 1)]
    path = path_a if max(p[1] for p in path_a) > max(p[1] for p in path_b) else path_b
    return np.array(path)


FB_PX = frame_bottom_curve()


def prof_np(y):
    return np.asarray(profile(float(y)).p, float)


class Sec:
    """a planar section at y: the port half's polyline (lat, z) keel -> crest, its arc lengths, helpers"""

    def __init__(self, y, pod=True):
        self.y = y
        P = profile(float(y), pod)
        self.prof = P
        self.p = np.asarray(P.p, float)
        seg = np.hypot(*np.diff(self.p, axis=0).T)
        self.s = np.concatenate([[0.0], np.cumsum(seg)])
        self.L = float(self.s[-1])

    def at(self, a):
        a = min(max(a, 0.0), self.L)
        return np.array([np.interp(a, self.s, self.p[:, 0]), np.interp(a, self.s, self.p[:, 1])])

    def pt3(self, a, sg=1.0):
        q = self.at(a)
        return Vector((sg * q[0], self.y, q[1]))

    def arc_z_outer(self, z):
        """arc position of the outermost point at height z (what SIDE sees)"""
        p, s = self.p, self.s
        best = None
        z0, z1 = p[:-1, 1], p[1:, 1]
        m = ((z0 - z) * (z1 - z) <= 0) & (z0 != z1)
        for i in np.nonzero(m)[0]:
            t = (z - z0[i]) / (z1[i] - z0[i])
            lat = p[i, 0] + (p[i + 1, 0] - p[i, 0]) * t
            if best is None or lat > best[0]:
                best = (lat, s[i] + (s[i + 1] - s[i]) * t)
        if best is None:
            return self.L if z > p[:, 1].max() else 0.0
        return best[1]

    def arc_lat_top(self, lat):
        p, s = self.p, self.s
        for i in range(len(p) - 1, 0, -1):
            l0, l1 = p[i, 0], p[i - 1, 0]
            if (l0 - lat) * (l1 - lat) <= 0 and l0 != l1:
                return s[i] + (s[i - 1] - s[i]) * (lat - l0) / (l1 - l0)
        return self.L

    def nearest_arc(self, q):
        a, b = self.p[:-1], self.p[1:]
        d = b - a
        L2 = np.maximum(1e-14, (d * d).sum(1))
        t = np.clip(((q - a) * d).sum(1) / L2, 0, 1)
        c = a + d * t[:, None]
        i = int(np.argmin(((c - q) ** 2).sum(1)))
        return self.s[i] + (self.s[i + 1] - self.s[i]) * t[i]


_SECS = {}


def sec_at(y, pod=True):
    key = (round(float(y), 6), pod)
    if key not in _SECS:
        _SECS[key] = Sec(y, pod)
    return _SECS[key]


# ================================================================== 3. the design lines in 3D (port half; starboard = mirror)
def _polyline_len(P):
    return np.concatenate([[0.0], np.cumsum(np.linalg.norm(np.diff(P, axis=0), axis=1))])


class Curve3:
    """a dense 3D polyline with nearest-point and arc-length helpers (a design line the relaxation keeps vertices on)"""

    def __init__(self, P, name=""):
        self.P = np.asarray(P, float)
        self.s = _polyline_len(self.P)
        self.L = float(self.s[-1])
        self.name = name

    def at(self, t):
        t = min(max(t, 0.0), self.L)
        return Vector(tuple(np.array([np.interp(t, self.s, self.P[:, k]) for k in range(3)])))

    def nearest(self, q, t_hint=None, window=0.4):
        q = np.asarray(tuple(q), float)
        if t_hint is None:
            lo, hi = 0, len(self.P) - 1
        else:
            lo = max(0, int(np.searchsorted(self.s, t_hint - window)) - 1)
            hi = min(len(self.P) - 1, int(np.searchsorted(self.s, t_hint + window)) + 1)
        a, b = self.P[lo:hi], self.P[lo + 1:hi + 1]
        d = b - a
        L2 = np.maximum(1e-14, (d * d).sum(1))
        t = np.clip(((q - a) * d).sum(1) / L2, 0, 1)
        c = a + d * t[:, None]
        i = int(np.argmin(((c - q) ** 2).sum(1)))
        return Vector(tuple(c[i])), float(self.s[lo + i] + (self.s[lo + i + 1] - self.s[lo + i]) * t[i])

    def mirrored(self):
        P = self.P.copy()
        P[:, 0] *= -1
        return Curve3(P, self.name + "_R")


def smooth3(C3, sigma, iters=2, keep_ends=True):
    """a design line made geodesically fair: Gaussian-smoothed along its length (sigma m; odd reflection at the ends keeps
    them and their tangents) and projected back onto the surface, twice"""
    P0 = C3.P
    for _ in range(iters):
        s_ = _polyline_len(P0)
        L_ = s_[-1]
        ds = sigma / 6.0
        n = max(8, int(L_ / ds))
        su = np.linspace(0.0, L_, n + 1)
        Pu = np.stack([np.interp(su, s_, P0[:, k]) for k in range(3)], 1)
        r = int(3 * sigma / (L_ / n)) + 1
        ext = np.vstack([2 * Pu[0] - Pu[r:0:-1], Pu, 2 * Pu[-1] - Pu[-2:-r - 2:-1]])
        w = np.exp(-0.5 * ((np.arange(-r, r + 1) * (L_ / n)) / sigma) ** 2)
        w /= w.sum()
        sm = np.stack([np.convolve(ext[:, k], w, mode="valid") for k in range(3)], 1)
        sm = np.array([tuple(proj(Vector(tuple(p)))[0]) for p in sm])
        if keep_ends:
            sm[0], sm[-1] = Pu[0], Pu[-1]
        P0 = sm
    return Curve3(P0, C3.name)


def _smooth_poly(P, iters=3):
    P = np.asarray(P, float).copy()
    for _ in range(iters):
        P[1:-1] = 0.25 * P[:-2] + 0.5 * P[1:-1] + 0.25 * P[2:]
    return P


# ---- the upper groove: four lines (UR0 band/frame edge, UR1 floor, UR2 floor, UR3 shell edge), from the duct round the nose's
# upper corner and across the nose top to the centre line. On the flank (to X_TURN0) each line is where SIDE draws it (the
# outermost point of the section at its fair height); then its centre line continues in PLAN (TOP: it runs on in along the
# nose's taper and turns, in a round G2 corner, onto the crossing at x X_CROSS, where TOP draws the groove's two lines at
# x 884.5 and 888.5), each point on the surface under it (top-down), and the four lines are offsets of it in the surface.
X_TURN0 = 868.0                 # where the plan turn starts (the side-view lines are kept aft of it)
X_CROSS = 886.6                 # the crossing's centre line (TOP: 884.5 / 888.5)
LAT_TURN_END = 0.42             # the turn has become the straight crossing here (m from the centre line)


def _flank_arcs(x):
    S_ = sec_at(Yi(x))
    a0 = S_.arc_z_outer(Zi(UR_L(x)))
    a3 = S_.arc_z_outer(Zi(UR_U(x)))
    if a3 - a0 < 0.03:
        m = 0.5 * (a0 + a3); a0, a3 = m - 0.015, m + 0.015
    return S_, (a0, a0 + WALL, a3 - WALL, a3)


def _quintic(p0, d0, dd0, p1, d1, dd1, n):
    """a quintic Hermite (positions, first and second derivatives at both ends), n+1 samples"""
    t = np.linspace(0.0, 1.0, n + 1)[:, None]
    h0 = 1 - 10 * t ** 3 + 15 * t ** 4 - 6 * t ** 5
    h1 = t - 6 * t ** 3 + 8 * t ** 4 - 3 * t ** 5
    h2 = 0.5 * t ** 2 - 1.5 * t ** 3 + 1.5 * t ** 4 - 0.5 * t ** 5
    h3 = 10 * t ** 3 - 15 * t ** 4 + 6 * t ** 5
    h4 = -4 * t ** 3 + 7 * t ** 4 - 3 * t ** 5
    h5 = 0.5 * t ** 3 - t ** 4 + 0.5 * t ** 5
    return h0 * p0 + h1 * d0 + h2 * dd0 + h3 * p1 + h4 * d1 + h5 * dd1


def build_upper_groove(turn_len=1.0):
    # flank: x from 545 to X_TURN0, every 0.5 px
    xs = np.arange(543.0, X_TURN0 + 1e-6, 0.5)               # (from the duct's front: the fair curve held 2 px)
    lines = [[], [], [], []]
    cen = []
    for x in xs:
        S_, arcs = _flank_arcs(x)
        for k in range(4):
            lines[k].append(tuple(S_.pt3(arcs[k])))
        cen.append(tuple(S_.pt3(0.5 * (arcs[0] + arcs[3]))))
    cen = np.array(cen)
    w_flank = float(np.linalg.norm(np.array(lines[3][-1]) - np.array(lines[0][-1])))
    # the plan path's start: position, tangent and second derivative of the centre line in plan (x_m, lat) at X_TURN0
    pl = np.stack([cen[:, 1], cen[:, 0]], 1)                       # (y, lat)
    p0 = pl[-1]
    # its tangent and second derivative, by arc length in plan (a quadratic through the last 9 samples)
    seg = np.linalg.norm(np.diff(pl[-9:], axis=0), axis=1)
    t_ = np.concatenate([[0.0], np.cumsum(seg)])
    cx = np.polyfit(t_, pl[-9:, 0], 2); cl = np.polyfit(t_, pl[-9:, 1], 2)
    T0 = np.array([np.polyval(np.polyder(cx), t_[-1]), np.polyval(np.polyder(cl), t_[-1])])
    A0 = np.array([np.polyval(np.polyder(cx, 2), t_[-1]), np.polyval(np.polyder(cl, 2), t_[-1])])
    nT = np.linalg.norm(T0)
    T0, A0 = T0 / nT, A0 / nT ** 2
    p1 = np.array([Yi(X_CROSS), LAT_TURN_END])
    chord = np.linalg.norm(p1 - p0)
    Lh = turn_len * chord
    d0v, dd0v = T0 * Lh, A0 * Lh * Lh
    d1v = np.array([0.0, -1.0]) * Lh
    turn = _quintic(p0, d0v, dd0v, p1, d1v, np.zeros(2), 160)
    cross = np.stack([np.full(60, Yi(X_CROSS)), np.linspace(LAT_TURN_END, 0.0, 61)[1:]], 1)
    plan = np.vstack([turn[1:], cross])                            # (y, lat)
    cen_t = []
    for y_, lat_ in plan:
        p_, n_ = on_top(float(lat_), float(y_))
        cen_t.append(tuple(p_))
    cen_t = np.array(cen_t)
    # offsets in the surface: the in-surface normal of the centre line, the width easing from the flank's to w_cross
    full_c = np.vstack([cen, cen_t])
    out = [list(l) for l in lines]
    n_f = len(cen)
    b_prev = (Vector(lines[3][-1]) - Vector(lines[0][-1])).normalized()     # toward the shell side (UR3) at the flank's end
    for q in range(len(cen_t)):
        i = n_f + q
        c = Vector(tuple(full_c[i]))
        t = Vector(tuple(full_c[min(i + 1, len(full_c) - 1)] - full_c[max(i - 1, 0)])).normalized()
        n_ = proj(c)[1]
        b = n_.cross(t).normalized()
        if b.dot(b_prev) < 0:                                       # continuous, toward the shell side (UR3)
            b = -b
        b_prev = b
        w = w_flank
        for k, off in enumerate((-0.5 * w, -0.5 * w + WALL, 0.5 * w - WALL, 0.5 * w)):
            p_ = proj(c + b * off)[0]
            out[k].append(tuple(p_))
    # (fair in the surface: the sections' outermost points kink where the surface under the line changes fast)
    curves = [smooth3(Curve3(np.array(o), "UR%d" % k), 0.04) for k, o in enumerate(out)]
    for c_ in curves:                                               # the crossing ends on the centre line
        c_.P[-1, 0] = 0.0
    return curves, smooth3(Curve3(full_c, "UR_c"), 0.04), w_flank


UR3D, URC3D, UR_W = build_upper_groove()
REPORT["upper_groove"] = {"width_m": round(UR_W, 4), "turn_from_x": X_TURN0, "cross_x": X_CROSS}
print("SKIN part 1-3", round(time.time() - SKIN_T0, 2), flush=True)