"""Aethel-9 hull: one quad grid (rows along the hull, columns round it) on the surface of ae_lib, every opening (seven
windows and the hatch) a rectangle of the grid. The grid starts as rows spaced by need and even columns; each opening's
rectangle of indices is pinned onto its drawn outline (4 corners, sides spaced by arc length, inside by a Coons patch) and
the rest of the grid follows by a Laplace solve of the displacement from the start grid, so it relaxes back away from the
openings. Pure numpy; build_mesh() returns world-space arrays for Blender.
"""
import numpy as np
import ae_lib as L

GUIDE = False
ROOF_WINDOW = False      # the concept (2026-10-22) has a vented hood where the blueprint's TOP draws a roof window
WINSLOW = 400
NJ = L.NJ_PROFILE                       # columns round the hull (even: phi 0 and 0.5 are columns, the grid is mirror-symmetric)


def base_rows():
    """Row stations (SIDE x): ~2.4 px apart over the nose and the windows, opening to ~6 px aft."""
    v = [L.V_TIP]
    while v[-1] < L.V_END:
        x = v[-1]
        step = 1.0 if x < 241 else 2.4 if x < 300 else (3.0 if x < 410 else (3.2 if 440 < x < 495 else 6.0))
        v.append(x + step)
    v = np.array(v)
    return L.V_TIP + (v - v[0]) * (L.V_END - L.V_TIP) / (v[-1] - v[0])


# ---- outlines -----------------------------------------------------------------------------------------------------------
def fillet_outline(poly, radii, corners, n_arc=10, n_edge=60):
    """Closed polygon with a fillet radius per vertex -> dense outline and the outline indices of the corner vertices
    (the middle of their fillets). corners: the polygon indices that are the grid rectangle's corners (4, in order)."""
    P = np.asarray(poly, float); n = len(P)
    pts, marks = [], {}
    arcs = []
    for k in range(n):
        a, b, c = P[k - 1], P[k], P[(k + 1) % n]
        r = radii[k]
        da, dc = a - b, c - b
        la, lc = np.linalg.norm(da), np.linalg.norm(dc)
        r = min(r, 0.45 * la, 0.45 * lc)
        p0, p2 = b + da / la * r, b + dc / lc * r
        arcs.append((p0, b, p2))
    for k in range(n):
        p0, b, p2 = arcs[k]
        t = np.linspace(0, 1, n_arc)
        q = ((1 - t) ** 2)[:, None] * p0 + (2 * (1 - t) * t)[:, None] * b + (t ** 2)[:, None] * p2
        start = len(pts)
        pts += list(q)
        if k in corners:
            marks[k] = start + n_arc // 2
        nxt = arcs[(k + 1) % n][0]
        t = np.linspace(0, 1, n_edge + 2)[1:-1]
        pts += list(p2[None] * (1 - t)[:, None] + nxt[None] * t[:, None])
    return np.array(pts), [marks[k] for k in corners]


def to_param(view, pts):
    if view == 'side':
        v, phi = L.side_to_param(pts[:, 0], pts[:, 1])
        return np.stack([v, phi], 1)
    out = []
    for a, b in pts:
        if view == 'front':
            out.append(L.view_to_param('front', a, b))
        else:   # top: (SIDE x, X)
            out.append(L.view_to_param('top', b, a))
    return np.array(out)


def mirror(par):
    return np.stack([par[:, 0], 1.0 - par[:, 1]], 1)


def window_defs():
    """Each opening: name, param outline (v, phi) with phi unwrapped near its centre, the 4 corner indices in the order
    (i0/j0, i1/j0, i1/j1, i0/j1) -- i along v (rows), j along phi (columns)."""
    W = []
    # FRONT upper centre: corners in FRONT terms: (-28,219) stbd-bottom, (-28,164) stbd-top, (28,164), (28,219)
    # stbd is phi > 0.5 on top; i0 = the bottom edge (near the tip), j0 = port edge (phi < 0.5)
    poly = [(28, 219), (28, 164), (-28, 164), (-28, 219)]
    o, c = fillet_outline(poly, [6] * 4, [0, 1, 2, 3])
    W.append(('win_front_upper', to_param('front', o), c))
    # FRONT lower centre on the belly: port edge phi > 0 (j1), stbd phi < 0 (j0); top edge (near the tip) = i0
    poly = [(-33, 254), (-33, 317), (33, 317), (33, 254)]
    o, c = fillet_outline(poly, [6] * 4, [0, 1, 2, 3])
    par = to_param('front', o); par[:, 1] = np.where(par[:, 1] > 0.5, par[:, 1] - 1.0, par[:, 1])
    W.append(('win_front_lower', par, c))
    # SIDE upper (port): order the corners i0/j0 = front-bottom (265,218), i1/j0 = aft-bottom (376.5,182),
    # i1/j1 = aft-top (376,156.5), i0/j1 = front-top (259,211)
    poly = [(265, 218), (311, 214), (376.5, 182), (376, 156.5), (320, 160), (300, 169), (272, 189), (259, 211)]
    o, c = fillet_outline(poly, [2.5, 6, 3, 3, 10, 12, 10, 2.5], [0, 2, 3, 7])
    par = to_param('side', o)
    W.append(('win_side_upper_L', par, c))
    W.append(('win_side_upper_R', mirror(par)[::-1].copy(), [len(o) - 1 - k for k in (c[3], c[2], c[1], c[0])]))
    # SIDE lower (port): i0/j0 = front-bottom (269,307), i1/j0 = aft-bottom (283,307), i1/j1 = aft-top (292,251),
    # i0/j1 = front-top (248,255)
    poly = [(269, 307), (275, 311), (283, 307), (296, 262), (292, 251), (248, 255), (248, 262)]
    o, c = fillet_outline(poly, [3, 4, 3, 5, 3, 3, 4], [0, 2, 4, 5])
    par = to_param('side', o)
    W.append(('win_side_lower_L', par, c))
    W.append(('win_side_lower_R', mirror(par)[::-1].copy(), [len(o) - 1 - k for k in (c[3], c[2], c[1], c[0])]))
    # TOP roof window: corners (SIDE x, X): i0/j0 = (316, +21) front-port, i1/j0 = (401, +33), i1/j1 = (401, -33), i0/j1
    poly = [(316, 21), (401, 33), (401, -33), (316, -21)]
    o, c = fillet_outline(poly, [5, 4, 4, 5], [0, 1, 2, 3])
    if ROOF_WINDOW:
        W.append(('win_roof', to_param('top', o), c))
    # hatch: a circle on top, corners at 45 deg
    v0, x0, r = L.HATCH
    # start at (front-port) = angle 135 deg from +v: v0 + r cos, X = r sin
    th = np.deg2rad(135) + np.linspace(0, -2 * np.pi, 400, endpoint=False)
    o = np.stack([v0 + r * np.cos(th), x0 + r * np.sin(th)], 1)
    c = [0, 100, 200, 300]   # 135 (front-port) -> 45 (aft-port) -> -45 (aft-stbd) -> -135 (front-stbd)
    W.append(('hatch', to_param('top', o), c))
    return W


# ---- grid planning ------------------------------------------------------------------------------------------------------
def arc_resample(curve3, curvep, n):
    """Spread n points along a param curve by 3D arc length; returns their param values."""
    s = np.concatenate([[0], np.cumsum(np.linalg.norm(np.diff(curve3, axis=0), axis=1))])
    t = np.linspace(0, s[-1], n)
    return np.stack([np.interp(t, s, curvep[:, 0]), np.interp(t, s, curvep[:, 1])], 1)


def p3(par):
    return L.surf_world(par[:, 0], par[:, 1])


def plan(rows=None):
    rows = base_rows() if rows is None else rows
    NI = len(rows)
    phi = np.arange(NJ) / NJ
    BV, BP = np.meshgrid(rows, phi, indexing='ij')
    pinV = np.zeros_like(BV) * np.nan
    pinP = np.zeros_like(BV) * np.nan
    rects = {}
    for name, par, c in window_defs():
        n = len(par)
        # side curves between corners, in outline order: c0->c1 (j0 edge, along v), c1->c2 (i1 edge), c2->c3 (j1), c3->c0 (i0)
        def seg(a, b):
            idx = np.arange(a, b + 1) if b >= a else np.concatenate([np.arange(a, n), np.arange(0, b + 1)])
            return par[idx]
        sides = [seg(c[0], c[1]), seg(c[1], c[2]), seg(c[2], c[3]), seg(c[3], c[0])]
        corner = par[c]
        vi = lambda vv: int(np.argmin(np.abs(rows - vv)))
        pj = lambda pp: int(np.round(pp * NJ))
        i0 = vi((corner[0, 0] + corner[3, 0]) / 2); i1 = vi((corner[1, 0] + corner[2, 0]) / 2)
        j0 = pj((corner[0, 1] + corner[1, 1]) / 2); j1 = pj((corner[2, 1] + corner[3, 1]) / 2)
        if i1 - i0 < 2 or j1 - j0 < 2:
            raise ValueError(f'{name}: rect too small {i0, i1, j0, j1}')
        ni, nj = i1 - i0 + 1, j1 - j0 + 1
        bj0 = arc_resample(p3(sides[0]), sides[0], ni)            # j = j0, i0 -> i1
        bi1 = arc_resample(p3(sides[1]), sides[1], nj)            # i = i1, j0 -> j1
        bj1 = arc_resample(p3(sides[2]), sides[2], ni)[::-1]      # j = j1, i0 -> i1
        bi0 = arc_resample(p3(sides[3]), sides[3], nj)[::-1]      # i = i0, j0 -> j1
        # Coons patch in param space
        s = np.linspace(0, 1, ni)[:, None, None]; t = np.linspace(0, 1, nj)[None, :, None]
        Cj0 = bj0[:, None, :]; Cj1 = bj1[:, None, :]; Ci0 = bi0[None, :, :]; Ci1 = bi1[None, :, :]
        P00, P10, P01, P11 = bj0[0], bj0[-1], bj1[0], bj1[-1]
        patch = (1 - t) * Cj0 + t * Cj1 + (1 - s) * Ci0 + s * Ci1 - (
            (1 - s) * (1 - t) * P00 + s * (1 - t) * P10 + (1 - s) * t * P01 + s * t * P11)
        jj = np.arange(j0, j1 + 1) % NJ
        for a, i in enumerate(range(i0, i1 + 1)):
            for b, j in enumerate(jj):
                if not np.isnan(pinV[i, j]):
                    raise ValueError(f'{name} overlaps another opening at row {i} col {j}')
                pinV[i, j] = patch[a, b, 0]
                pinP[i, j] = patch[a, b, 1]
        rects[name] = (i0, i1, j0, j1)
        if name.startswith('win_front'):
            # the front windows' near edge crosses the rows (a level line on a steep nose): fill the columns between the
            # nose cap (row 0) and that edge by straight interpolation, so the rows there cannot fold
            # (columns kept at the edge's phi, the nose cap's ring follows them)
            for b, j in enumerate(jj):
                for i in range(0, i0):
                    f = i / i0
                    pinV[i, j] = (1 - f) * rows[0] + f * patch[0, b, 0]
                    pinP[i, j] = patch[0, b, 1]
    # guide: the column between the lower front window and the lower side window runs diagonally between them (their
    # facing edges are parallel in (v, phi) but cross the grid at right angles); pinned along a path from the nose cap,
    # above the front window's corner, to the side window's front-bottom corner. Mirrored for starboard.
    jf = rects['win_front_lower'][3]; isl, _, jsl, _ = rects['win_side_lower_L']
    ic = rects['win_front_lower'][0]
    corner = (pinV[ic, jf % NJ], pinP[ic, jf % NJ])
    fb = (pinV[isl, jsl], pinP[isl, jsl])
    for jg, sgn in (((jf + 1, 1), (NJ - jf - 1, -1)) if GUIDE else ()):
        for i in range(0, isl):
            if i <= ic:
                f = i / ic
                v_ = (1 - f) * rows[0] + f * (corner[0] + 0.7)
                p_ = corner[1] + 0.011
            else:
                f = (i - ic) / (isl - ic)
                v_ = (1 - f) * (corner[0] + 0.7) + f * fb[0]
                p_ = (1 - f) * (corner[1] + 0.011) + f * fb[1]
            pinV[i, jg] = v_
            pinP[i, jg] = p_ if sgn > 0 else 1.0 - p_
    # unwrap pinned phi to be near the base column's phi
    d = pinP - BP
    d = np.where(d > 0.5, d - 1, np.where(d < -0.5, d + 1, d))
    pinned = ~np.isnan(pinV)
    DV = np.where(pinned, pinV - BV, 0.0)
    DP = np.where(pinned, d, 0.0)
    fixedV = pinned.copy(); fixedV[0] = fixedV[-1] = True   # the end rings keep their station, slide round freely
    fixedP = pinned
    for _ in range(6000):
        nV = 0.25 * (np.roll(DV, 1, 1) + np.roll(DV, -1, 1) + np.vstack([DV[:1], DV[:-1]]) + np.vstack([DV[1:], DV[-1:]]))
        nP = 0.25 * (np.roll(DP, 1, 1) + np.roll(DP, -1, 1) + np.vstack([DP[:1], DP[:-1]]) + np.vstack([DP[1:], DP[-1:]]))
        DV = np.where(fixedV, DV, nV); DP = np.where(fixedP, DP, nP)
    GV = BV + DV; GP = BP + DP
    GV[0] = L.V_TIP; GV[-1] = L.V_END
    if WINSLOW:
        GV, GP = winslow(GV, GP, fixedV, fixedP, WINSLOW)
    return dict(rows=rows, V=GV, P=GP, rects=rects, NI=NI)


def winslow(GV, GP, fixedV, fixedP, iters, k=500.0):
    """Winslow (elliptic) smoothing of the free vertices in (v, phi*k): unlike Laplace it keeps cells from folding where
    pinned edges lie at an angle to the grid. Columns wrap; the end rings keep their v."""
    x = GV.copy(); y = GP * k
    for _ in range(iters):
        def sh(a, di, dj):
            b = np.roll(a, -dj, 1)
            if di == 1: b = np.vstack([b[1:], b[-1:]])
            if di == -1: b = np.vstack([b[:1], b[:-1]])
            return b
        out = []
        xe = (sh(x, 0, 1) - sh(x, 0, -1)) / 2; ye = (sh(y, 0, 1) - sh(y, 0, -1)) / 2
        # unwrap phi differences across the seam
        ye = np.where(ye > k / 4, ye - k / 2, np.where(ye < -k / 4, ye + k / 2, ye))
        xx = (sh(x, 1, 0) - sh(x, -1, 0)) / 2; yx = (sh(y, 1, 0) - sh(y, -1, 0)) / 2
        a = xe * xe + ye * ye; b = xx * xe + yx * ye; c = xx * xx + yx * yx
        for f in (x, y):
            fpi, fmi, fpj, fmj = sh(f, 1, 0), sh(f, -1, 0), sh(f, 0, 1), sh(f, 0, -1)
            fpp, fpm, fmp, fmm = sh(f, 1, 1), sh(f, 1, -1), sh(f, -1, 1), sh(f, -1, -1)
            if f is y:   # seam: bring neighbours within half a turn of the centre
                fix = lambda n: np.where(n - f > k / 2, n - k, np.where(n - f < -k / 2, n + k, n))
                fpi, fmi, fpj, fmj, fpp, fpm, fmp, fmm = map(fix, (fpi, fmi, fpj, fmj, fpp, fpm, fmp, fmm))
            cross = (fpp - fpm - fmp + fmm) / 4
            new = (a * (fpi + fmi) + c * (fpj + fmj) - 2 * b * cross) / (2 * (a + c) + 1e-12)
            out.append(new)
        x = np.where(fixedV, x, 0.5 * x + 0.5 * out[0])
        y = np.where(fixedP, y, 0.5 * y + 0.5 * out[1])
    return x, y / k


def hole_mask(g):
    """Faces (i, j) -> True where the face lies inside an opening (face (i,j) spans rows i..i+1, columns j..j+1)."""
    NI = g['NI']
    m = np.zeros((NI - 1, NJ), bool)
    for name, (i0, i1, j0, j1) in g['rects'].items():
        for i in range(i0, i1):
            for j in range(j0, j1):
                m[i, j % NJ] = True
    return m


def flipped(g):
    """Faces whose param-space orientation is reversed (fold-overs)."""
    V, P = g['V'], g['P']
    a = np.stack([V, P * 600.0], -1)
    p00 = a[:-1]; p10 = a[1:]; p01 = np.roll(a, -1, 1)[:-1]
    d1 = p10 - p00; d2 = p01 - p00
    d2[..., 1] = np.where(d2[..., 1] < -300, d2[..., 1] + 600, d2[..., 1])
    cr = d1[..., 0] * d2[..., 1] - d1[..., 1] * d2[..., 0]
    return cr <= 0


if __name__ == '__main__':
    g = plan()
    print('rows', g['NI'], 'cols', NJ)
    for k, r in g['rects'].items():
        print(k, r)
    f = flipped(g) & ~hole_mask(g)
    print('flipped faces (outside openings)', int(f.sum()), np.argwhere(f)[:10])