# Grouper, round 3: the skin's layout, its relaxation and its mesh (exec'd by gr_hull.py after gr_skin.py, same namespace).
# See gr_skin.py's head for what and why.

# ---- FRAME_BOT (FB0 the drawn line, shell side; FB1 30 mm up it, the frame side): from the band's front wall along the flank
# as SIDE draws it, then (forward of X_FB_TURN0) a contour of the surface at the height it reaches the nose's front (SIDE: y
# 389.5-391), round the nose's lower corner and across the raked face just above the flat front, to the centre line.
X_FB_TURN0 = 884.0
Z_FB_CROSS = Zi(390.4)


def _fb_y(x):
    P = FB_PX
    m = P[:, 0] >= 840.0
    xs, ys = P[m, 0], P[m, 1]
    o = np.argsort(xs)
    return float(np.interp(max(x, 840.0), xs[o], ys[o]))


def _offset_curve(P, d, up_hint):
    """a curve offset by d in the surface (along the in-surface normal, turned toward up_hint at the start)"""
    out = []
    bprev = None
    for i in range(len(P)):
        c = Vector(tuple(P[i]))
        t = Vector(tuple(P[min(i + 1, len(P) - 1)] - P[max(i - 1, 0)])).normalized()
        n = proj(c)[1]
        b = n.cross(t).normalized()
        if bprev is None:
            if b.dot(up_hint) < 0:
                b = -b
        elif b.dot(bprev) < 0:
            b = -b
        bprev = b
        out.append(tuple(proj(c + b * d)[0]))
    return np.array(out)


def build_frame_bottom():
    pts = []
    xs = np.arange(X_BF0, X_FB_TURN0 + 1e-6, 0.5)
    for x in xs:
        S_ = sec_at(Yi(x))
        z = Zi(_fb_y(x))
        w = smoothstep(870.0, X_FB_TURN0, x)
        z = z + (Z_FB_CROSS - z) * w
        pts.append(tuple(S_.pt3(S_.arc_z_outer(z))))
    x = X_FB_TURN0
    while x < 912.9:
        x = min(912.9, x + 0.25)
        S_ = sec_at(Yi(x))
        p = S_.pt3(S_.arc_z_outer(Z_FB_CROSS))
        n = proj(p)[1]
        if n.y < -0.72:                                          # the surface faces forward: continue across the face
            break
        pts.append(tuple(p))
    lat0 = pts[-1][0]
    for q in range(1, 81):
        lat = lat0 * (1 - q / 80.0)
        hit, nrm, fi, dist = SKIN.ray_cast(Vector((lat, Y_NOSE - 1.0, Z_FB_CROSS)), Vector((0, 1, 0)))
        if hit is not None:
            pts.append(tuple(hit))
    P = np.array(pts)
    for _ in range(6):                                           # fair the junction of the contour and the face's crossing
        P[1:-1] = 0.25 * P[:-2] + 0.5 * P[1:-1] + 0.25 * P[2:]
        P = np.array([tuple(proj(Vector(tuple(p)))[0]) for p in P])
    P[-1, 0] = 0.0
    # (the second pass) the turn round the nose's lower corner made fair: a Gaussian of 9 cm along the line, on the surface (the
    # contour met the crossing of the face in a 1-2 cm corner, and the line's vertices turned 25 deg each there)
    P = smooth3(Curve3(P, "FB0"), 0.09, iters=2).P.copy()
    P[-1, 0] = 0.0
    fb1 = _offset_curve(P, 0.030, Vector((0, 0, 1)))
    fb1[-1, 0] = 0.0
    return Curve3(P, "FB0"), Curve3(_smooth_poly(fb1, 2), "FB1")


FB0_3D, FB1_3D = build_frame_bottom()


# ================================================================== 4. the layout's numbers
# port half columns j (0 keel .. J crest); the full ring k = 0 .. M-1 (port k = j, starboard k = M - j)
J_POD = (4, 7)                   # on the pod (its full length; easing out inside its tail cone and nose pieces): the lower inner
                                 # side next to the lower fillet (145 deg round the nacelle from its outer side) and its bottom
# (the second layout) Below the lower groove the section turns about 480 deg in all, 340 of them over the pod and its two
# fillets: the first layout's 12 columns there left 30-70 deg between neighbours and a faceted nacelle. 18 columns now, the
# extra ones taken from the deck, where neighbours differed by under 4 deg.
J_LR = (18, 19, 20, 21)          # the lower groove: shell edge, floor, floor, band edge
J_DB = 24                        # the door's ring: the rectangle's bottom (a band column; two rows between it and the groove,
                                 # which passes under the door from its aft corner, 3 cm off, to 0.45 m under its fore end)
J_FB = (25, 26)                  # FRAME_BOT: shell edge, frame edge (free aft of the band's front wall)
J_UR = (32, 33, 34, 35)          # the upper groove: band / frame edge, floor, floor, shell edge
K_DECK = 5                       # deck rows wrapped round the nose with the groove (the groove's corner turns over them)
J_D = J_UR[3] + K_DECK           # the top poles' column (40)
J_P = 8                          # the bottom poles' column: the rows above it, up to FRAME_BOT, wrap round the chin
K_CHIN = J_FB[0] - J_P
N_TOP = J_P                      # the cap's half width in v-lines = the belly's columns = the deck core's columns
J = J_D + N_TOP                  # the crest (48; M a multiple of 4: the duct's end plate is a Coons patch)
M = 2 * J                        # 96 columns round the hull
HALF = J
N_SIDE = J_D - J_P               # the cap's height in u-lines (32)
A_SQ = 0.43                      # the collar's square (plan) mapped onto the bore's circle
R_BORE_HOLE = 0.512
J_HC = J - 6                     # the collar square's columns (12 across the crest)
X_W = X_BF0 + 1.9                # the band's front wall: an ordinary interval (32 mm), its forward station the last planar one
X_CT, X_CB = 866.0, 896.0        # the cap's top side (across the deck) and bottom side (across the underside)
N_TR = 14                        # transition rings from X_W to the ring C (the cap's boundary; FRAME_BOT's turn
                                 # round the nose's lower corner runs through them)
X_DA, X_DB = 733.0, 801.0        # the door's rectangle (its stations)


def vj(k):
    k %= M
    return k if k <= HALF else M - k


def side_of(k):
    k %= M
    return 1.0 if k <= HALF else -1.0


REPORT["skin_layout"] = {"M": M, "J_D": J_D, "J_P": J_P, "N_TOP": N_TOP, "N_SIDE": N_SIDE, "K_DECK": K_DECK, "K_CHIN": K_CHIN}

# ================================================================== 5. the lower groove in 3D (port), and stations by need
# The lower groove's four lines on each section: where SIDE draws its two edges (the outermost point at each height). On the pod
# (x 780-830) the drawn lines plunge round the nacelle to its bottom; the lines keep to its lower outer side: the centre's angle
# round the nacelle's axis eases smoothly toward LR_TH_MIN (round 2 clamped it hard and squeezed the four lines to 2 cm), the
# groove keeps its width, and its depth runs out over x 824-830 (as round 2).
LR_TH_MIN = math.radians(-50.0)


def _softmax(a, b, k=0.12):
    return 0.5 * (a + b + math.sqrt((a - b) ** 2 + k * k)) - 0.5 * k


def lr_arcs(x, S_):
    zl, zu = Zi(LR_L(min(x, X_LR_END))), Zi(LR_U(min(x, X_LR_END)))
    P = S_.prof
    if P.c is not None and P.r > 0.2 and x > 785.0:              # only where the groove runs down the nacelle's outer side
        cz, r = POD_Z, pod_r_true(min(x, 836.0))
        tl = math.asin(max(-1.0, min(1.0, (zl - cz) / r)))
        tu = math.asin(max(-1.0, min(1.0, (zu - cz) / r)))
        tc = 0.5 * (tl + tu)
        tc2 = _softmax(tc, LR_TH_MIN)
        if tc2 > tc + 1e-9:                                       # lifted off the plunge: keep the width (in angle)
            d = tc2 - tc
            tl, tu = tl + d, tu + d
        zl, zu = cz + r * math.sin(tl), cz + r * math.sin(tu)
    a0, a3 = S_.arc_z_outer(zl), S_.arc_z_outer(zu)
    if a3 - a0 < 0.03:
        m = 0.5 * (a0 + a3); a0, a3 = m - 0.015, m + 0.015
    return (a0, a0 + WALL, a3 - WALL, a3)


def build_lower_groove():
    lines = [[], [], [], []]
    for x in np.arange(543.0, X_LR_END + 1e-6, 0.5):
        S_ = sec_at(Yi(x))
        arcs = lr_arcs(x, S_)
        for k in range(4):
            lines[k].append(tuple(S_.pt3(arcs[k])))
    return [smooth3(Curve3(np.array(l), "LR%d" % k), 0.04) for k, l in enumerate(lines)]


LR3D = build_lower_groove()


def geodesic_turn_need(C3, turn):
    """for a 3D line on the surface: the station spacing (m, along y) at each of its points that keeps its turning per station
    at `turn` (its geodesic curvature: the turning in the surface's tangent plane)"""
    P = C3.P
    s = C3.s
    T = np.gradient(P, s, axis=0)
    T /= np.maximum(1e-12, np.linalg.norm(T, axis=1))[:, None]
    dT = np.gradient(T, s, axis=0)
    N = np.array([tuple(proj(Vector(tuple(p)))[1]) for p in P])
    B = np.cross(N, T)
    kg = np.abs((dT * B).sum(1))
    kg = np.convolve(kg, np.ones(5) / 5, mode="same")
    dyds = np.abs(np.gradient(P[:, 1], s))
    dsdy = 1.0 / np.maximum(1e-3, dyds)
    h = turn / np.maximum(1e-6, kg * dsdy)
    return P[:, 1], h


TURN_LINE, TURN_SIL = math.radians(4.5), math.radians(5.0)
H_MAX, H_MIN, GROW = 0.11, 0.030, 0.28
_LINE_NEEDS = []
for C3 in list(UR3D) + list(LR3D):
    yy, hh = geodesic_turn_need(C3, TURN_LINE)
    m = (yy > Yi(X_W) - 1e-6)                                    # the planar part only
    _LINE_NEEDS.append((yy[m], hh[m]))


def _curv_fn(f, y, e=0.004):
    a, b, c = f(y - e), f(y), f(y + e)
    d1 = (c - a) / (2 * e)
    d2 = (c - 2 * b + a) / (e * e)
    return abs(d2) / (1 + d1 * d1) ** 1.5, math.sqrt(1 + d1 * d1)


def need_h(x):
    y = Yi(x)
    h = H_MAX
    for f in (TOPf, BOTf, HALFf):
        k, ds = _curv_fn(f, y)
        if k > 1e-6:
            h = min(h, TURN_SIL / (k * ds))
    return max(H_MIN, h)


_xs_need = np.arange(543.0, X_W + 0.26, 0.25)
_h_raw = np.array([need_h(x) for x in _xs_need])
for yy, hh in _LINE_NEEDS:                            # the design lines' needs, binned onto the same x grid
    xq = X0 - yy / S
    iq = np.clip(np.round((xq - 543.0) / 0.25).astype(int), 0, len(_xs_need) - 1)
    np.minimum.at(_h_raw, iq, np.maximum(H_MIN, hh))
_h = _h_raw.copy()
for _ in range(2):                                    # a Lipschitz envelope: the spacing grows at most GROW per metre of hull
    for q in range(1, len(_h)):
        _h[q] = min(_h[q], _h[q - 1] + GROW * 0.25 * S)
    for q in range(len(_h) - 2, -1, -1):
        _h[q] = min(_h[q], _h[q + 1] + GROW * 0.25 * S)
_dens = np.concatenate([[0.0], np.cumsum(0.25 * S / _h[1:])])


def _x_of_n(nv):
    return float(np.interp(nv, _dens, _xs_need))


def _n_of_x(x):
    return float(np.interp(x, _xs_need, _dens))


FIXED_X = sorted({543.0, round(CL_IMG - A_SQ / S, 4), round(CL_IMG + A_SQ / S, 4), X_DA, X_DB, X_BF0, X_W})
STX = []
for xa, xb in zip(FIXED_X, FIXED_X[1:]):
    na, nb = _n_of_x(xa), _n_of_x(xb)
    n = max(1, int(round(nb - na)))
    STX += [_x_of_n(na + (nb - na) * q / n) for q in range(n)]
STX.append(FIXED_X[-1])
for q, xf in enumerate(STX):                          # the fixed ones exactly
    for xf2 in FIXED_X:
        if abs(xf - xf2) < 1e-6:
            STX[q] = xf2
STX = [float(x) for x in STX]
I_W = len(STX) - 1                                    # the last planar station (X_W); stations run from the duct forward
I_BF = STX.index(X_BF0)
I_DA, I_DB = STX.index(X_DA), STX.index(X_DB)
I_CA, I_CB = STX.index(round(CL_IMG - A_SQ / S, 4)), STX.index(round(CL_IMG + A_SQ / S, 4))   # aft, fore
_dx = np.diff(STX)
REPORT["skin_stations"] = {"planar": len(STX), "min_mm": round(float(_dx.min()) * S * 1000, 1),
                           "max_mm": round(float(_dx.max()) * S * 1000, 1),
                           "max_ratio": round(float(max(max(a / b, b / a) for a, b in zip(_dx, _dx[1:]))), 3)}
print("SKIN stations", REPORT["skin_stations"], flush=True)