# Mako: everything on the hull (exec'd by build_mako.py after mk_hull.py). Quads only: lofts and foils with three-quad end
# caps, revolved parts with Coons caps, cube spheres, ladder patches (mk_topo.py); parts that sink into the hull are pruned.
# Numbers are reference pixels from mk_measure.py unless they are in metres.
from mathutils.bvhtree import BVHTree

att = MB("_attached", [PAINT, DARK, GRAPH, FIN, SENSOR, STEEL, RUBBER, YELLOW, LENS, LAMP])
hull.bm.normal_update()
HBVH = BVHTree.FromBMesh(hull.bm)

VIEW = {"top": (lambda u, v: Vector((u, v, 5.0)), Vector((0, 0, -1))),
        "bot": (lambda u, v: Vector((u, v, -5.0)), Vector((0, 0, 1))),
        "sideP": (lambda u, v: Vector((5.0, u, v)), Vector((-1, 0, 0))),
        "sideS": (lambda u, v: Vector((-5.0, u, v)), Vector((1, 0, 0)))}

def proj(view, u, v):
    """(u, v) of a view's 2D coordinates (metres) -> (point on the hull, outward normal)"""
    of, d = VIEW[view]
    hit, nrm, fi, dist = HBVH.ray_cast(of(u, v), d)
    if hit is None:
        return None, None
    n = normal_at(hit) if abs(nrm.dot(normal_at(hit))) > 0.5 else nrm.normalized()
    if n.dot(-d) < 0:
        n = -n
    return hit, n

def region(outline, axis):
    o = 1 - axis
    ts = [p[axis] for p in outline]
    t0, t1 = min(ts), max(ts)
    pad = (t1 - t0) * 0.02
    def lo_hi(t):
        xs = []
        n = len(outline)
        for i in range(n):
            a, b = outline[i], outline[(i + 1) % n]
            if (a[axis] - t) * (b[axis] - t) <= 0 and a[axis] != b[axis]:
                xs.append(a[o] + (b[o] - a[o]) * (t - a[axis]) / (b[axis] - a[axis]))
        return (min(xs), max(xs)) if xs else (0.0, 0.0)
    return lo_hi, t0 + pad, t1 - pad

def patch(mb, outline, view, mat, axis=0, lift=0.004, bulge=0.0, rim=None, rim_r=0.008, nt=8, nacross=3, cell=0.035):
    """a thin plate on the hull over a 2D outline, as a ladder of quads (no fan); optional rim tube round it"""
    lo_hi, t0, t1 = region(outline, axis)
    wmax = max(lo_hi(t0 + (t1 - t0) * i / 8)[1] - lo_hi(t0 + (t1 - t0) * i / 8)[0] for i in range(9))
    nt = max(nt, int(math.ceil((t1 - t0) / cell)))
    nacross = max(nacross, int(math.ceil(wmax / cell)))
    def to3d(t, s):
        uv = (t, s) if axis == 0 else (s, t)
        p, n = proj(view, *uv)
        return p
    def nrm(p):
        _, n = proj(view, *((p.x, p.y) if view in ("top", "bot") else (p.y, p.z)))
        return n if n is not None else normal_at(p)
    ladder_patch(mb, lo_hi, t0, t1, to3d, mat, nt=nt, nacross=nacross, lift=lift, bulge=bulge, normal=nrm)
    if rim is not None:
        chain = [(t0 + (t1 - t0) * i / nt, lo_hi(t0 + (t1 - t0) * i / nt)) for i in range(nt + 1)]
        loop2 = [(t, lh[0]) for t, lh in chain] + [(t, lh[1]) for t, lh in reversed(chain)]
        pts = []
        for (t, s) in loop2 + [loop2[0]]:
            uv = (t, s) if axis == 0 else (s, t)
            p, n = proj(view, *uv)
            pts.append(p + n * rim_r * 0.35)
        mb.sweep(pts, rim_r, rim, sides=6, caps=False)

def band_outline(pts2, half):
    """a polyline (2D) thickened to a closed band of half width `half` with round ends"""
    n = len(pts2)
    L, R = [], []
    for j in range(n):
        a = Vector(pts2[max(j - 1, 0)]); b = Vector(pts2[min(j + 1, n - 1)])
        t = (b - a).normalized()
        nn = Vector((-t.y, t.x))
        p = Vector(pts2[j])
        L.append(p + nn * half); R.append(p - nn * half)
    def cap(c, t, sgn):
        nn = Vector((-t.y, t.x))
        return [c + (nn * math.cos(a) + t * math.sin(a) * sgn) * half for a in [math.pi * q / 4 for q in (1, 2, 3)]]
    t_end = (Vector(pts2[-1]) - Vector(pts2[-2])).normalized()
    t_beg = (Vector(pts2[1]) - Vector(pts2[0])).normalized()
    out = L + cap(Vector(pts2[-1]), t_end, 1) + R[::-1] + cap(Vector(pts2[0]), -t_beg, 1)
    return [(p.x, p.y) for p in out]

def fin_from(le, te, tip, th0, th1, mat, n=14, s0=-0.06, s1=0.975, side_normal=None):
    """a triangular fin: root chord le..te, both edges running to the tip; lens sections, quads, three-quad end caps"""
    pl = (te - le).cross(tip - le).normalized()
    if side_normal is not None and pl.dot(side_normal) < 0:
        pl = -pl
    secs = []
    for j in range(n + 1):
        s = s0 + (s1 - s0) * j / n
        a = le + (tip - le) * s
        b = te + (tip - te) * s
        c = (a + b) * 0.5
        ch = (b - a).length
        secs.append((c, (b - a).normalized(), pl, ch, th0 + (th1 - th0) * max(0.0, s)))
    foil(att, secs, mat)
    return secs

def ysec(le_tab, te_tab, z):
    """the chord's ends at height z (side view) from a leading edge and a trailing edge polyline [(img x, img y)]"""
    def at(tab):
        pts = sorted(tab, key=lambda p: -p[1])           # by height
        for (xa, va), (xb, vb) in zip(pts, pts[1:]):
            za, zb = Zi(va), Zi(vb)
            if min(za, zb) - 1e-9 <= z <= max(za, zb) + 1e-9 and za != zb:
                return Yi(xa + (xb - xa) * (z - za) / (zb - za))
        p = pts[0] if abs(z - Zi(pts[0][1])) < abs(z - Zi(pts[-1][1])) else pts[-1]
        return Yi(p[0])
    return at(le_tab), at(te_tab)

def vfin(le_tab, te_tab, z0, z1, th0, th1, mat, n=12, root_extra=0.12, down=False):
    """a fin in the centre plane: chords at heights z0 (root, extended into the hull) .. z1 (tip)"""
    secs = []
    for j in range(n + 1):
        t = j / n
        z = z0 + (z1 - z0) * t
        a, b = ysec(le_tab, te_tab, z)
        if j == 0:                                          # carry the root's chord down into the hull
            zr = z0 + (root_extra if down else -root_extra)
            secs.append((Vector((0, (a + b) / 2, zr)), Vector((0, 1, 0)), Vector((1, 0, 0)), abs(b - a), th0))
        if j == n:                                          # a short rounded tip, not a needle
            ch = max(abs(b - a), 0.035)
        else:
            ch = abs(b - a)
        secs.append((Vector((0, (a + b) / 2, z)), Vector((0, 1, 0)), Vector((1, 0, 0)), ch, th0 + (th1 - th0) * t))
    foil(att, secs, mat)

# ------------------------------------------------------------------ side pods: water-jet nacelles (TOP: grey outboard; FRONT: dark pods)
PODS = {}
def pod(side):
    x0, x1, lat, yc, hw, hh = POD["x0"], POD["x1"], POD["lat"], POD["y"], POD["hw"], POD["hh"]
    cx, cz = side * lat * S, Zi(yc)
    NP = 16
    def ring(t, k_w, k_h, dy=0.0):
        """rounded-rectangle section (superellipse n 2.6) at a share t of the pod's length"""
        x = x0 + (x1 - x0) * t
        y = Yi(x) + dy
        out = []
        for q in range(NP):
            a = 2 * math.pi * (q + 0.5) / NP
            c, s = math.cos(a), math.sin(a)
            out.append(Vector((cx + hw * S * k_w * math.copysign(abs(c) ** (2 / 2.6), c),
                               y, cz + hh * S * k_h * math.copysign(abs(s) ** (2 / 2.6), s))))
        return out
    # outer skin: full from 18 % to 55 %, rounding into the intake lip at the front and the nozzle at the back
    prof = [(0.00, 0.84, 0.84), (0.015, 0.95, 0.95), (0.05, 1.0, 1.0), (0.20, 1.0, 1.0), (0.45, 0.98, 0.98), (0.65, 0.88, 0.86),
            (0.82, 0.70, 0.68), (0.94, 0.56, 0.54), (1.00, 0.52, 0.50)]
    rings = [ring(t, kw, kh) for t, kw, kh in prof]
    # intake: the lip turns in and runs 0.22 m back to a dark stator face
    lip = [ring(0.0, 0.72, 0.72, 0.012), ring(0.0, 0.64, 0.64, 0.05), ring(0.0, 0.62, 0.62, 0.22)]
    noz = [ring(1.0, 0.44, 0.42, -0.012), ring(1.0, 0.40, 0.38, -0.16)]
    grid_rings(att, lip[::-1] + rings + noz, [DARK, DARK] + [GRAPH] + [DARK] * (len(rings) - 1) + [GRAPH, GRAPH], closed_u=True,
               expect=(4, 0, (1, 0, 0)))
    # stator faces deep in the intake and the nozzle: Coons caps (quads)
    for r_, outward in ((lip[-1], -1), (noz[-1], 1)):
        bm2 = bmesh.new()
        vs = [bm2.verts.new(p) for p in r_]
        coons_cap(bm2, vs if outward < 0 else vs[::-1])
        bm2.normal_update()
        if sum(f.normal.y for f in bm2.faces) * outward < 0:
            bmesh.ops.reverse_faces(bm2, faces=bm2.faces[:])
        att.merge(bm2, SENSOR)
    # stator vanes across the intake and the nozzle
    for yv, kk in ((Yi(x0) + 0.20, 0.62), (Yi(x1) - 0.14, 0.40)):
        for a in (0.0, math.pi / 2):
            qbox(att, (2 * hw * S * kk * abs(math.cos(a)) + 0.02, 0.02, 2 * hh * S * kk * abs(math.sin(a)) + 0.02),
                 (cx, yv, cz), GRAPH)
    PODS[side] = (cx, cz)
for side in (1, -1):
    pod(side)

# ------------------------------------------------------------------ big fins: long grey spikes, low on the flank, drooped about 60 deg
def fin_pts(side):
    (lx, ly), (tx, ty), (px_, plat, py) = FIN_LE_ROOT, FIN_TE_ROOT, FIN_TIP
    le = on_side(side, Yi(lx), Zi(ly)) - Vector((side * 0.03, 0, 0))
    te = on_side(side, Yi(tx), Zi(ty)) - Vector((side * 0.03, 0, 0))
    tip = Vector((side * plat * S, Yi(px_), Zi(py)))
    return le, te, tip
FIN_SECS = {}
for side in (1, -1):
    le, te, tip = fin_pts(side)
    FIN_SECS[side] = fin_from(le, te, tip, 0.12, 0.022, FIN, n=16, side_normal=Vector((side, 0, 0.3)))
    # root fairing: a fillet along the root (the SIDE view's lighter "elbow" panel)
    pts = [le.lerp(te, t) for t in [j / 10 for j in range(11)]]
    Pp = [on_side(side, p.y, p.z) for p in pts]
    bm2 = bmesh.new()
    R = []
    n = len(Pp)
    for j in range(n):
        t = (Pp[min(j + 1, n - 1)] - Pp[max(j - 1, 0)]).normalized()
        N = normal_at(Pp[j]); N = (N - N.dot(t) * t).normalized()
        b = N.cross(t)
        u = j / (n - 1)
        w = 0.05 + 0.07 * math.sin(math.pi * min(1.0, u * 1.1))
        h = 0.03 + 0.05 * math.sin(math.pi * min(1.0, u * 1.1))
        sec_pts = [(-w, -0.02), (-0.7 * w, 0.6 * h), (-0.25 * w, h), (0.25 * w, h), (0.7 * w, 0.6 * h), (w, -0.02),
                   (0.4 * w, -0.03), (-0.4 * w, -0.03)]
        R.append([bm2.verts.new(Pp[j] + b * uu + N * vv) for uu, vv in sec_pts])
    ring_rows(bm2, R)
    for r in (R[0], R[-1]):
        bm2.faces.new((r[0], r[1], r[2], r[7])); bm2.faces.new((r[7], r[2], r[3], r[6])); bm2.faces.new((r[6], r[3], r[4], r[5]))
    bm2.normal_update()
    bmesh.ops.recalc_face_normals(bm2, faces=bm2.faces[:])
    att.merge(bm2, FIN)

# ------------------------------------------------------------------ dorsal fin, tail fin (orange) and the ventral stern fin (grey)
vfin(DORSAL_LE, DORSAL_TE, Zi(670.0), Zi(549.0), 0.085, 0.018, PAINT, n=12)
vfin(TAIL_LE, TAIL_TE, Zi(760.0), Zi(659.0), 0.060, 0.015, PAINT, n=9)
VLE = [VENTRAL[0], VENTRAL[1], VENTRAL[2], VENTRAL[3]]
VTE = [VENTRAL[3], VENTRAL[4], VENTRAL[5], VENTRAL[6]]
vfin(VLE, VTE, Zi(850.0), Zi(940.0), 0.055, 0.014, FIN, n=9, down=True)

# ------------------------------------------------------------------ the jaw: two curved grooves between its slats, two prongs per side
for sg, vw in ((1, "sideP"), (-1, "sideS")):
    for g in GROOVES:
        pts2 = [(Yi(x), Zi(v)) for x, v in g]
        patch(att, band_outline(pts2, 0.020), vw, SENSOR, axis=0, lift=0.003, rim=GRAPH, rim_r=0.007, nt=10, nacross=2)
for side in (1, -1):
    for (rx, ry), (tx, ty), lat in PRONGS:                 # a pair of flat slats, the jaw's flank carried on below its keel
        p_side = Vector((side * lat, Yi(rx), Zi(ry)))       # where the SIDE view shows the slat leave the jaw
        tip = Vector((side * (lat + 0.02), Yi(tx), Zi(ty)))
        ax = (tip - p_side).normalized()
        # the slat runs nearly along the jaw's flank, so its root bends in: one section well inside the jaw (in and up),
        # one where it leaves the skin, then straight to the tip
        p_in = p_side + Vector((-side * 0.12, -0.10, 0.18))
        def frame(a):
            cd = Vector((0, 1, 0)); cd = (cd - cd.dot(a) * a).normalized()
            nd = Vector((side, 0, 0)); nd = (nd - nd.dot(a) * a - nd.dot(cd) * cd).normalized()
            return cd, nd
        cd0, nd0 = frame((p_side - p_in).normalized())
        cd, nd = frame(ax)
        secs = [(p_in, cd0, nd0, 0.090, 0.040)]
        for j in range(10):
            t = j / 9
            c = p_side.lerp(tip, t)
            ch = 0.10 - 0.035 * smoothstep(0.35, 1.0, t)
            secs.append((c, cd, nd, ch if j < 9 else 0.05, 0.042 - 0.012 * t))
        foil(att, secs, PAINT)

# ------------------------------------------------------------------ belly pack (SIDE: the dark box under the rear bay), ribbed, on two skids
bx0, bx1, bhw, bbot = BELLY
byc, blen = Yi((bx0 + bx1) / 2), (bx1 - bx0) * S
bz0 = Zi(bbot)
bz1 = hull_z_at(0.0, byc, top=False) + 0.10
qbox(att, (2 * bhw * S, blen, bz1 - bz0), (0, byc, (bz0 + bz1) / 2), DARK, bevel=0.04)
for t in (0.18, 0.5, 0.82):
    qbox(att, (2 * bhw * S + 0.03, 0.035, 0.08), (0, Yi(bx0) + blen * t, bz0 + 0.12), GRAPH, bevel=0.01)
for side in (1, -1):
    qbox(att, (0.05, blen + 0.10, 0.045), (side * (bhw * S - 0.06), byc, bz0 - 0.02), RUBBER, bevel=0.015)

# ------------------------------------------------------------------ stern housing: the tail's last section carried on, with the thruster duct
sx0, sx1 = STERN_HOUSING
yE = Y_TAIL
hrows = []
for t, k in ((0.0, 1.02), (0.25, 1.04), (0.7, 1.0), (1.0, 0.94)):
    y = Yi(sx0) + (Yi(sx1) - Yi(sx0)) * t
    yy = min(y, Y_TAIL)
    row = []
    zm = ZMf(Y_TAIL)
    for q in range(16):
        a = 2 * math.pi * (q + 0.5) / 16
        p = sec(Y_TAIL, a if a <= math.pi / 2 or a >= 3 * math.pi / 2 else a)
        row.append(Vector((p.x * k, y, zm + (p.z - zm) * k)))
    hrows.append(row)
grid_rings(att, hrows, DARK, closed_u=True, expect=(1, 4, (0, 0, 1)))
bm2 = bmesh.new()
vs = [bm2.verts.new(p) for p in hrows[-1]]
coons_cap(bm2, vs[::-1])
bm2.normal_update()
if sum(f.normal.y for f in bm2.faces) < 0:
    bmesh.ops.reverse_faces(bm2, faces=bm2.faces[:])
att.merge(bm2, GRAPH)
TH_C = Vector((0, Yi(sx1), ZMf(Y_TAIL) - 0.02))          # thruster axis on the housing's aft face
def duct_ring(r, dy, n=24):
    return [TH_C + Vector((r * math.cos(2 * math.pi * q / n), dy, r * math.sin(2 * math.pi * q / n))) for q in range(n)]
# one profile, outer wall -> lip -> inner wall, skinned in order so the inner wall faces the axis
grid_rings(att, [duct_ring(0.150, -0.01), duct_ring(0.150, 0.11), duct_ring(0.138, 0.125), duct_ring(0.122, 0.11),
                 duct_ring(0.122, -0.03)], [GRAPH, GRAPH, STEEL, DARK], closed_u=True, expect=(0, 6, (0, 0, 1)))
qcyl(att, 0.045, 0.05, TH_C + Vector((0, 0.005, 0)), GRAPH, rot=(math.pi / 2, 0, 0), segs=12)   # hub stator, ends 3 cm before the rotor
for a in range(4):
    ang = math.pi / 4 + a * math.pi / 2
    qbox(att, (0.012, 0.03, 0.115), TH_C + Vector((0.07 * math.cos(ang), 0.02, 0.07 * math.sin(ang))), GRAPH, rot=(0, -ang + math.pi / 2, 0))
# small fittings on the housing: a sensor lens above, a light below
qsphere(att, Vector((0, Yi(sx1) - 0.02, ZMf(Y_TAIL) + 0.22)), (0.045, 0.03, 0.045), SENSOR, cuts=2)
qcyl(att, 0.035, 0.02, Vector((0, Yi(sx1) + 0.005, ZMf(Y_TAIL) - 0.22)), LAMP, rot=(math.pi / 2, 0, 0), segs=12)

# ------------------------------------------------------------------ canopy mullions: low graphite bars over the glass, flush with the frame rims
frames = MB("SM_Mako_Frames", [GRAPH])
def mullion(pts3):
    rows = []
    for (w, h) in ((-0.014, -0.022), (-0.014, 0.010), (0.014, 0.010), (0.014, -0.022)):
        row = []
        for j, p in enumerate(pts3):
            t = (pts3[min(j + 1, len(pts3) - 1)] - pts3[max(j - 1, 0)]).normalized()
            nrm = normal_at(p)
            bn = nrm.cross(t).normalized()
            row.append(p + bn * w + nrm * h)
        rows.append(row)
    grid_rings(frames, [list(r) for r in zip(*rows)], GRAPH, closed_u=False,
               expect=(len(pts3) // 2, 1, tuple(normal_at(pts3[len(pts3) // 2]))))
ytop0, ytop1 = Yi(75.0), YS[I_HOOD] + 0.004
mullion([sec(ytop0 + (ytop1 - ytop0) * t, math.pi / 2) for t in [j / 24 for j in range(25)]])
TH_SIDE = math.radians(34.0)
yside1 = next(YS[i] for i in range(I_HOOD, I_C0 + 1) if HOODth(YS[i]) < TH_SIDE + 0.10)
for sg in (1, -1):
    ys_ = [Yi(92.0) + (yside1 - Yi(92.0)) * t for t in [j / 24 for j in range(25)]]
    mullion([sec(y, TH_SIDE if sg > 0 else math.pi - TH_SIDE) for y in ys_])

# ------------------------------------------------------------------ lamps: a pair low on the jaw's front, looking ahead and down
LAMPS = []
for side in (1, -1):
    y = Yi(178.0)
    z = CHINz(y) - 0.20
    p = on_side(side, y, z)
    n = normal_at(p)
    aim = n                                              # flush with the jaw: the bezel sits as a ring round the lens
    rot = Vector((0, 0, 1)).rotation_difference(aim).to_euler()
    revolve(att, [(0.062, -0.025), (0.062, 0.010), (0.055, 0.018), (0.046, 0.010)], [GRAPH, GRAPH, GRAPH], segs=16,
            xf=Matrix.Translation(p) @ rot.to_matrix().to_4x4(), cap0=False, cap1=True, cap1_mat=LENS)
    qcyl(att, 0.040, 0.004, p + aim * 0.006, LAMP, rot=tuple(rot), segs=16)
    LAMPS.append((p + aim * 0.03, (n + Vector((0, -1.6, -0.3))).normalized()))   # the beam looks ahead and down
REPORT["stage"] = "fittings"