# Grouper, round 3: the skin's columns, the cap, the hole rings, the relaxation and the mesh (exec'd by gr_hull.py after
# gr_skin_layout.py; same namespace). See gr_skin.py's head.

# ================================================================== 6. columns on the planar stations
RELAX_L = 0.30                   # a pin's pull fades over this length past the end of its feature
COL_K_LEN = 0.8                  # the free columns crowd where the section curves: density sqrt(1 + (k * COL_K_LEN)^2)


def plane_cross(C3, y):
    P = C3.P
    d = P[:, 1] - y
    idx = np.nonzero(d[:-1] * d[1:] <= 0)[0]
    if not len(idx):
        return None
    i = idx[0]
    t = d[i] / (d[i] - d[i + 1]) if d[i] != d[i + 1] else 0.0
    return P[i] + (P[i + 1] - P[i]) * t


def _arc_of_pt(S_, p):
    return S_.nearest_arc(np.array((abs(p[0]), p[2])))


def _ray_section(S_, c2, d2):
    """the farthest crossing of a ray from c2 along d2 with the section's polyline (lat, z), or None"""
    p = S_.p
    A, B = p[:-1] - c2, p[1:] - c2
    E = B - A
    den = d2[0] * E[:, 1] - d2[1] * E[:, 0]
    with np.errstate(divide="ignore", invalid="ignore"):
        t = (A[:, 0] * E[:, 1] - A[:, 1] * E[:, 0]) / den
        u = (A[:, 0] * d2[1] - A[:, 1] * d2[0]) / den
    ok = (np.abs(den) > 1e-12) & (u >= 0) & (u <= 1) & (t > 0)
    if not ok.any():
        return None
    return c2 + d2 * float(t[ok].max())


def door_outline_px(inset_px=0.0, n=10, expo=3.0):
    """the doorway (SIDE px): a rounded rectangle whose corners are quarter superellipses (exponent 3: the curvature starts at
    zero at the straight sides, G2), inset by inset_px (negative: outset)"""
    x0, x1 = DOOR["x"][0] + inset_px, DOOR["x"][1] - inset_px
    y0, y1 = DOOR["y"][0] + inset_px, DOOR["y"][1] - inset_px
    r = DOOR["corner_r"] * 1.25 - inset_px                       # (a superellipse corner reads smaller than a circle's)
    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)
            c, s_ = math.cos(a), math.sin(a)
            pts.append((cx + r * math.copysign(abs(c) ** (2 / expo), c), cy + r * math.copysign(abs(s_) ** (2 / expo), s_)))
    return pts


def _door_bottom_y(x):
    """the door ring's rectangle bottom (SIDE y px): between the lower groove's band edge and the doorway's hole edge"""
    ol = door_outline_px(-1.4, 24)
    ys = [p[1] for p in ol if abs(p[0] - x) < 0.8]
    yh = max(ys) if ys else DOOR["y"][1] + 1.4
    yg = LR_U(x)
    return 0.45 * yh + 0.55 * yg if yg - yh < 6.0 else yh + 2.7


def station_pins(x, S_):
    """{j: (arc, weight)} at a planar station"""
    y = Yi(x)
    L = S_.L
    pins = {0: (0.0, 1.0), J: (L, 1.0)}
    for k in range(4):
        p = plane_cross(UR3D[k], y)
        if p is not None:
            pins[J_UR[k]] = (_arc_of_pt(S_, p), 1.0)
    if x <= X_LR_END:
        for k in range(4):
            p = plane_cross(LR3D[k], y)
            if p is not None:
                pins[J_LR[k]] = (_arc_of_pt(S_, p), 1.0)
    else:
        S_e = sec_at(Yi(X_LR_END))
        w = 1.0 - smoothstep(0.0, RELAX_L, Yi(X_LR_END) - y)
        for k in range(4):
            p = plane_cross(LR3D[k], Yi(X_LR_END - 0.01))
            if p is not None and w > 0:
                pins[J_LR[k]] = (_arc_of_pt(S_e, p) / S_e.L * L, w)
    if x >= X_BF0 - 1e-6:
        for jj, C3 in ((J_FB[0], FB0_3D), (J_FB[1], FB1_3D)):
            p = plane_cross(C3, y)
            if p is not None:
                pins[jj] = (_arc_of_pt(S_, p), 1.0)
    else:
        S_e = sec_at(Y_BF0)
        w = 1.0 - smoothstep(0.0, RELAX_L, y - Y_BF0)
        for jj, C3 in ((J_FB[0], FB0_3D), (J_FB[1], FB1_3D)):
            p = plane_cross(C3, Y_BF0)
            if p is not None and w > 0:
                pins[jj] = (_arc_of_pt(S_e, p) / S_e.L * L, w)
    # the pod's lines (the lower inner side, the bottom): where a ray from the nacelle's axis meets the section, full on the
    # visible nacelle (x 768-834), easing out inside its tail cone (752-768) and nose (834-846) pieces
    wp = smoothstep(752.0, 768.0, x) * (1.0 - smoothstep(834.0, 846.0, x))
    if wp > 0.0:
        for jj, th in zip(J_POD, (-145.0, -90.0)):
            a = math.radians(th)
            d2 = np.array((math.cos(a), math.sin(a)))
            c2 = np.array((POD_L, POD_Z))
            q = _ray_section(S_, c2, d2)
            if q is not None:
                pins[jj] = (S_.nearest_arc(q), wp)
    # the door ring's rectangle bottom, over the door's stations
    xa, xb = X_DA, X_DB
    if xa - 4.0 <= x <= xb + 4.0:
        xc = min(max(x, xa), xb)
        S_c = sec_at(Yi(xc))
        a = S_c.arc_z_outer(Zi(_door_bottom_y(xc)))
        w = 1.0 if xa <= x <= xb else 1.0 - smoothstep(0.0, 4.0 * S, abs(Yi(x) - Yi(xc)))
        pins[J_DB] = (a / S_c.L * L, w)
    # the collar's square (plan) over the clamp axis
    yc0, yc1 = Y_CL - A_SQ, Y_CL + A_SQ
    if yc0 - RELAX_L <= y <= yc1 + RELAX_L:
        yc = min(max(y, yc0), yc1)
        S_c = sec_at(yc)
        a = S_c.arc_lat_top(A_SQ)
        w = 1.0 if yc0 <= y <= yc1 else 1.0 - smoothstep(0.0, RELAX_L, abs(y - yc))
        pins[J_HC] = (a / S_c.L * L if not (yc0 <= y <= yc1) else S_.arc_lat_top(A_SQ), w)
    return pins


def _density(S_):
    """cumulative curvature-weighted arc length along the section"""
    p, s = S_.p, S_.s
    t = np.diff(p, axis=0)
    ang = np.arctan2(t[:, 1], t[:, 0])
    dang = np.abs((np.diff(ang) + np.pi) % (2 * np.pi) - np.pi)
    seg = np.diff(s)
    k = np.zeros(len(p))
    k[1:-1] = dang / np.maximum(1e-6, 0.5 * (seg[:-1] + seg[1:]))
    k = np.convolve(k, np.ones(9) / 9, mode="same")
    rho = np.sqrt(1 + (k * COL_K_LEN) ** 2)
    cum = np.concatenate([[0.0], np.cumsum(0.5 * (rho[:-1] + rho[1:]) * seg)])
    return cum


def _solve(ctrl, S_, cum):
    keys = sorted(ctrl)
    vals = [ctrl[k] for k in keys]
    for q in range(1, len(vals)):
        vals[q] = max(vals[q], vals[q - 1] + 0.002)
    out = [0.0] * (J + 1)
    for (ka, va), (kb, vb) in zip(zip(keys, vals), zip(keys[1:], vals[1:])):
        ca, cb = np.interp(va, S_.s, cum), np.interp(vb, S_.s, cum)
        for j in range(ka, kb + 1):
            cj = ca + (cb - ca) * (j - ka) / (kb - ka)
            out[j] = float(np.interp(cj, cum, S_.s))
    return out


def column_arcs(x):
    S_ = sec_at(Yi(x))
    cum = _density(S_)
    pins = station_pins(x, S_)
    strong = {j: a for j, (a, w) in pins.items() if w >= 1.0 - 1e-9}
    soft = {j: (a, w) for j, (a, w) in pins.items() if 0.0 < w < 1.0 - 1e-9}
    free = _solve(strong, S_, cum)
    if soft:
        c2 = dict(strong)
        for j, (a, w) in soft.items():
            c2[j] = free[j] + (a - free[j]) * w
        out = _solve(c2, S_, cum)
    else:
        out = free
    return out, set(strong)


_RAW = [column_arcs(x) for x in STX]
# the free columns smoothed along the hull (fractions of their section's length, a Gaussian over the neighbouring stations)
SMOOTH_SIG = 0.05
SMOOTH_SIG_LOW = 0.16           # below the lower groove: the pod's lobe adds girth there; its columns ease round it gently
_Ls = [sec_at(Yi(x)).L for x in STX]
_F = np.array([[a / L for a in r[0]] for r, L in zip(_RAW, _Ls)])
_yst = np.array([Yi(x) for x in STX])
COLS = []
_W = {}
for sg_ in (SMOOTH_SIG, SMOOTH_SIG_LOW):
    D_ = _yst[:, None] - _yst[None, :]
    W_ = np.exp(-0.5 * (D_ / sg_) ** 2)
    W_[np.abs(D_) > 3.5 * sg_] = 0.0
    _W[sg_] = W_
for i in range(len(STX)):
    pins = _RAW[i][1]
    f = _F[i].copy()
    for j in range(1, J):
        if j in pins:
            continue
        w = _W[SMOOTH_SIG_LOW if j < J_LR[0] else SMOOTH_SIG][i]
        f[j] = float((w * _F[:, j]).sum() / w.sum())
    L = _Ls[i]
    gap = 0.004 / L
    for j in range(1, J + 1):
        if j not in pins:
            f[j] = max(f[j], f[j - 1] + gap)
    for j in range(J - 1, 0, -1):
        if j not in pins:
            f[j] = min(f[j], f[j + 1] - gap)
    COLS.append([fj * L for fj in f])
PINNED = [r[1] for r in _RAW]

# ================================================================== 7. the mesh's topology: vertices with their constraints
SKV = []                        # base positions (Vector, on the analytic surface)
SKC = []                        # constraints: ("fix",) ("sec", y) ("crv", Curve3, t) ("srf",)
SKTW = []                       # the port twin of a starboard vertex (-1: port or centre-line vertex itself)
SKROLE = []                     # (kind, j or b, x) for the levels: ("col", j, x) ("hole", region, x)


def _nv(p, con, role, twin=-1):
    SKV.append(Vector(p))
    SKC.append(con)
    SKROLE.append(role)
    SKTW.append(twin)
    return len(SKV) - 1


def mirror(p):
    return Vector((-p[0], p[1], p[2]))


# ---- CYL: planar stations (port columns, then their starboard twins)
NI_P = len(STX)
I_C = I_W + N_TR                # the ring C (the cap's boundary)
CV = [[None] * M for _ in range(I_C + 1)]
# Free zones: where a design line runs steeply across the stations (the grooves' hooks into the duct ring, the lower groove's
# plunge onto the pod, the pod's tail where its lobe grows out of the hull) a planar station shears the quads beside the line
# into twisted parallelograms. There the stations may bend: the free vertices relax on the surface and the lines' vertices
# slide along their lines, so the quads beside a line square up to it.
FREE_ZONES = [z for z in os.environ.get("GROUPER_ZONES", "stern,plunge,podtail").split(",") if z]


def free_zone(x, j):
    if "stern" in FREE_ZONES and x < 562.0:
        return True
    if "plunge" in FREE_ZONES and X_DB + 2.0 < x < X_W - 1.0 and j < J_FB[0]:
        return True
    if "podtail" in FREE_ZONES and 742.0 < x < 790.0 and j < J_LR[0]:
        return True
    return False


def _pin_curve(j, x):
    if j in J_UR:
        return UR3D[J_UR.index(j)]
    if j in J_LR and x <= X_LR_END:
        return LR3D[J_LR.index(j)]
    if j == J_FB[0] and x >= X_BF0:
        return FB0_3D
    if j == J_FB[1] and x >= X_BF0:
        return FB1_3D
    return None


for i, x in enumerate(STX):
    S_ = sec_at(Yi(x))
    for j in range(J + 1):
        p = S_.pt3(COLS[i][j], 1.0)
        if j == 0 or j == J:
            p.x = 0.0
        zone = i > 0 and 0 < j < J and free_zone(x, j)
        if i == 0 or j in (0, J):
            con = ("fix",)
        elif j in PINNED[i]:
            C3 = _pin_curve(j, x) if zone and x < X_DA else None      # (the lower groove's plunge keeps its pins)
            if C3 is not None:
                q, tt = C3.nearest(p)
                p = q
                con = ("crv", C3, tt)
            elif zone and j in J_POD:
                con = ("srf",)
            else:
                con = ("fix",)
        else:
            con = ("srf",) if zone else ("sec", Yi(x))
        CV[i][j] = _nv(p, con, ("col", j, x))
    for j in range(1, J):
        CV[i][M - j] = _nv(mirror(SKV[CV[i][j]]), ("slave",), ("col", j, x), CV[i][j])

# ---- C: the cap's boundary ring (port half: the deck part, the side part, the underside part)
_ST = column_arcs(X_CT)[0]
_SB = column_arcs(X_CB)[0]
S_CT, S_CB = sec_at(Yi(X_CT)), sec_at(Yi(X_CB))
P_TOP = S_CT.pt3(_ST[J_D])                      # the top pole (port)
P_BOT = S_CB.pt3(_SB[J_P])                      # the bottom pole (port)


def _catmull(pts, n_per=40):
    P = [np.array(tuple(p), float) for p in pts]
    P = [2 * P[0] - P[1]] + P + [2 * P[-1] - P[-2]]
    out = []
    for i in range(1, len(P) - 2):
        p0, p1, p2, p3 = P[i - 1], P[i], P[i + 1], P[i + 2]
        for k in range(n_per):
            t = k / n_per
            out.append(0.5 * ((2 * p1) + (-p0 + p2) * t + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t * t + (-p0 + 3 * p1 - 3 * p2 + p3) * t ** 3))
    out.append(P[-2])
    return np.array(out)


def _turn_mid(C3, x_from):
    """the point of a design line where it has turned half way round the nose's corner (plan heading 45 deg between forward
    and inboard), searched forward of x_from"""
    P = C3.P
    T = np.gradient(P, axis=0)
    ang = np.degrees(np.arctan2(-T[:, 0], -T[:, 1]))          # 0 = forward (-Y), 90 = inboard (-X) for the port side
    i = int(np.argmin(np.abs(ang - 45.0) + 1000.0 * (P[:, 1] > Yi(x_from))))
    return Vector(tuple(P[i]))


W_G = _turn_mid(URC3D, 872.0)
W_F = _turn_mid(FB0_3D, 896.0)
_hm, _nm_ = SKIN.ray_cast(Vector((0.585, Y_NOSE - 2.0, Zi(368.0))), Vector((0, 1, 0)))[:2]
W_M = _hm
_side_raw = _catmull([P_BOT, W_F, W_M, W_G, P_TOP], 50)        # from the bottom pole up (the columns' order)
C_SIDE = smooth3(Curve3(np.array([tuple(proj(Vector(tuple(p)))[0]) for p in _side_raw]), "C_side"), 0.05, iters=2)


def _cross_param(Cside, C3):
    """the parameter along Cside where the design line C3 crosses it (their closest approach)"""
    best = (9.0, 0.0)
    for p in C3.P[::2]:
        q, t = Cside.nearest(Vector(tuple(p)))
        d = (q - Vector(tuple(p))).length
        if d < best[0]:
            best = (d, t)
    return best[1]


_cp = {J_P: 0.0, J_D: C_SIDE.L}
for jj, C3 in ((J_FB[0], FB0_3D), (J_FB[1], FB1_3D)):
    _cp[jj] = _cross_param(C_SIDE, C3)
for k in range(4):
    _cp[J_UR[k]] = _cross_param(C_SIDE, UR3D[k])
_kk = sorted(_cp)
_tt = [_cp[k] for k in _kk]
for q in range(1, len(_tt)):
    _tt[q] = max(_tt[q], _tt[q - 1] + 0.004)
C_PTS = {}
for (ka, ta), (kb, tb) in zip(zip(_kk, _tt), zip(_kk[1:], _tt[1:])):
    for j in range(ka, kb + 1):
        C_PTS[j] = C_SIDE.at(ta + (tb - ta) * (j - ka) / (kb - ka))
for j in range(J_D, J + 1):
    C_PTS[j] = S_CT.pt3(_ST[j])
for j in range(0, J_P + 1):
    C_PTS[j] = S_CB.pt3(_SB[j])
C_PTS[0].x = 0.0
C_PTS[J].x = 0.0

# ---- the transition rings (I_W+1 .. I_C): each column from its point on the last planar station to its point on C
_XT = {}
for j in range(J + 1):
    if j >= J_D:
        x_end = X_CT
    elif j <= J_P:
        x_end = X_CB
    else:
        x_end = None
    _XT[j] = x_end


def _trans_pt(j, t):
    """column j's point on transition ring t (0 = the last planar station, N_TR = C): along its own planar column for the deck
    core and the belly, a straight blend projected onto the surface for the flank's columns (relaxed afterwards)"""
    g = t / N_TR
    if _XT[j] is not None:
        x = X_W + (_XT[j] - X_W) * g
        if t == N_TR:
            return C_PTS[j]
        p = sec_at(Yi(x)).pt3(_col_cached(x)[j])
        if j in (0, J):
            p.x = 0.0
        return p
    a = SKV[CV[I_W][j]]
    b = C_PTS[j]
    return proj(a.lerp(b, g))[0]


_CA_CACHE = {}


def _col_cached(x):
    key = round(x, 6)
    if key not in _CA_CACHE:
        _CA_CACHE[key] = column_arcs(x)[0]
    return _CA_CACHE[key]


_TRCOLS = {}
for t in range(1, N_TR + 1):
    for j in range(J + 1):
        _TRCOLS[(t, j)] = _trans_pt(j, t)


def _crv_for(j):
    if j in J_UR:
        return UR3D[J_UR.index(j)]
    if j == J_FB[0]:
        return FB0_3D
    if j == J_FB[1]:
        return FB1_3D
    return None


CREST3 = Curve3(np.array([(0.0, Yi(x), TOPf(Yi(x))) for x in np.arange(X_W, 913.0, 0.2)]), "crest")
KEEL3 = Curve3(np.array([(0.0, Yi(x), BOTf(Yi(x))) for x in np.arange(X_W, 913.0, 0.2)]), "keel")
for t in range(1, N_TR + 1):
    i = I_W + t
    for j in range(J + 1):
        p = _TRCOLS[(t, j)]
        C3 = _crv_for(j)
        if j == 0:
            con = ("crv", KEEL3, KEEL3.nearest(p)[1])
        elif j == J:
            con = ("crv", CREST3, CREST3.nearest(p)[1])
        elif C3 is not None:
            q, tt = C3.nearest(p)
            p = q
            con = ("crv", C3, tt)
        else:
            con = ("srf",)
        CV[i][j] = _nv(p, con, ("col", j, X0 - p.y / S))
    for j in range(1, J):
        CV[i][M - j] = _nv(mirror(SKV[CV[i][j]]), ("slave",), ("col", j, X0 - SKV[CV[i][j]].y / S), CV[i][j])

# ---- F: the cap (a = 0 .. 2 N_TOP across, starboard to port; b = 0 .. N_SIDE up, from the bottom side to the top side)
NA = 2 * N_TOP
FV = [[None] * (N_SIDE + 1) for _ in range(NA + 1)]
for b in range(N_SIDE + 1):
    FV[NA][b] = CV[I_C][J_P + b]                    # port side
    FV[0][b] = CV[I_C][M - (J_P + b)]               # starboard side
for q in range(N_TOP + 1):
    FV[NA - q][N_SIDE] = CV[I_C][J_D + q]           # top side (port half; q = N_TOP: the crest)
    FV[q][N_SIDE] = CV[I_C][(M - (J_D + q)) % M]
    FV[NA - q][0] = CV[I_C][J_P - q]                # bottom side (port half; q = N_TOP: the keel)
    FV[q][0] = CV[I_C][(M - (J_P - q)) % M]
def build_nose_centre(x_top, x_bot):
    """the nose's centre line (x = 0) forward of the deck: over the nose top, down the raked face and the flat front, aft under
    the chin (the cap's middle v-line slides on it)"""
    pts = []
    x = x_top
    while x < 912.99:
        pts.append((0.0, Yi(x), TOPf(Yi(x))))
        x += 0.2
    zt, zb = TOPf(Y_NOSE), BOTf(Y_NOSE)
    for q in range(0, 41):
        pts.append((0.0, Y_NOSE, zt + (zb - zt) * q / 40.0))
    x = 912.99
    while x > x_bot - 1e-6:
        pts.append((0.0, Yi(x), BOTf(Yi(x))))
        x -= 0.2
    return Curve3(np.array(pts), "NOSE_C")


NOSE_C = build_nose_centre(X_CT, X_CB)
O_NOSE = Vector((0.0, Yi(890.0), 0.10))


def _coons(a, b):
    u, v = a / NA, b / N_SIDE
    B_, T_ = SKV[FV[a][0]], SKV[FV[a][N_SIDE]]
    L_, R_ = SKV[FV[0][b]], SKV[FV[NA][b]]
    c00, c10, c01, c11 = SKV[FV[0][0]], SKV[FV[NA][0]], SKV[FV[0][N_SIDE]], SKV[FV[NA][N_SIDE]]
    return ((1 - v) * B_ + v * T_ + (1 - u) * L_ + u * R_
            - ((1 - u) * (1 - v) * c00 + u * (1 - v) * c10 + (1 - u) * v * c01 + u * v * c11))


def _radial(p):
    d = (p - O_NOSE).normalized()
    hit = SKIN.ray_cast(O_NOSE, d)
    return hit[0] if hit[0] is not None else proj(p)[0]


for a in range(N_TOP, NA):                           # the port half and the centre column (a = N_TOP)
    for b in range(1, N_SIDE):
        p = _radial(_coons(a, b))
        jj = J_P + b
        C3 = _crv_for(jj)
        if a == N_TOP and C3 is not None:
            # (a design line's row crosses the centre line where its own curve does: its curve ends there, at x = 0; sliding
            # on the centre line left a V in FRAME_BOT, 26 deg, and in the upper groove at the nose's middle)
            p = Vector(tuple(C3.P[-1]))
            p.x = 0.0
            con = ("fix",)
        elif a == N_TOP:
            p.x = 0.0
            q, tt = NOSE_C.nearest(p)
            con = ("crv", NOSE_C, tt)
            p = q
        elif C3 is not None:
            q, tt = C3.nearest(p)
            p = q
            con = ("crv", C3, tt)
        else:
            con = ("srf",)
        FV[a][b] = _nv(p, con, ("col", jj, X0 - p.y / S))
for a in range(1, N_TOP):
    for b in range(1, N_SIDE):
        tw = FV[NA - a][b]
        FV[a][b] = _nv(mirror(SKV[tw]), ("slave",), SKROLE[tw], tw)

# ---- faces (quads, as index tuples) with a tag: ("cyl", i, k) / ("cap", a, b)
SKF = []
SKFT = []
for i in range(I_C):
    for k in range(M):
        SKF.append((CV[i][k], CV[i][(k + 1) % M], CV[i + 1][(k + 1) % M], CV[i + 1][k]))
        SKFT.append(("cyl", i, k))
for a in range(NA):
    for b in range(N_SIDE):
        SKF.append((FV[a][b], FV[a + 1][b], FV[a + 1][b + 1], FV[a][b + 1]))
        SKFT.append(("cap", a, b))
REPORT["skin_topology"] = {"verts": len(SKV), "quads": len(SKF), "stations_planar": NI_P, "transition": N_TR}
print("SKIN topology", REPORT["skin_topology"], round(time.time() - SKIN_T0, 1), flush=True)