"""The canopy: a glass bubble cut along the sill, its dark frame (round the rear rim and along both sills), as numpy grids.

canopy(prm) -> (glass verts/faces, frame verts/faces); prm is a dict (see DEFAULT); tg_fit_canopy.py fits it to the picture
and writes build/canopy.json, which tg_geom.py reads.

The bubble's section at station y is an ellipse-like ring: crest z_c above, bottom z_b below, widest line z_m, half width w;
the glass is the part above the sill: from the sill on the port flank over the crest to the sill on the starboard flank. Ahead
of the jaw's front (y_j) the sill runs down to the bubble's bottom, so the glass closes under the nose there.
"""
import math
import json
import numpy as np
from pathlib import Path

HERE = Path(__file__).resolve().parent

# the Mako sheet's TOP plan of the bubble (half width px over image x 41..286), normalised: (t along the canopy, w / w_max)
_PLAN = [(41, 0.0), (43, 1.5), (44, 2.5), (46, 4.5), (48, 7.5), (50, 10.0), (54, 15.0), (58, 19.0), (62, 23.0), (66, 27.0),
         (70, 31.0), (80, 39.5), (90, 47.0), (100, 54.5), (110, 61.0), (120, 67.0), (130, 73.0), (140, 78.0), (150, 83.0),
         (160, 87.0), (170, 91.0), (180, 95.0), (190, 98.0), (200, 100.0), (220, 102.0), (240, 103.0), (260, 101.0),
         (286, 99.0)]
PLAN_T = np.array([(x - 41) / 245.0 for x, _ in _PLAN])
PLAN_W = np.array([w / 103.0 for _, w in _PLAN])

DEFAULT = dict(
    y_nose=-2.105, z_nose=0.17,
    y_rear=-0.30,
    w_max=0.62, w_pow=1.0,                   # plan: the sheet's shape scaled to w_max, sharpened by w_pow
    crest=[0.40, 0.56, 0.68, 0.74, 0.74],    # crest z at t = 0.12, 0.3, 0.5, 0.75, 1.0 (t: nose 0 .. rear 1)
    bot=[-0.02, -0.15, -0.22],               # bubble bottom z at t = 0.12, 0.3, 0.5 (then flat)
    mid_frac=0.55,                           # widest line: crest - mid_frac * (crest - bottom)
    y_jaw=-1.45,                             # where the sill comes down to the bubble's bottom (the jaw's front)
    sill_rear=0.05,                          # sill z at the rear rim
    sill_mid=-0.08,                          # sill z half way between the jaw's front and the rear
    ex_top=1.0, ex_bot=0.8, ex_side=0.85,    # section superellipse exponents
)
T_CREST = [0.0, 0.12, 0.3, 0.5, 0.75, 1.0]
T_BOT = [0.0, 0.12, 0.3, 0.5, 1.0]


def load(name="canopy.json"):
    p = dict(DEFAULT)
    f = HERE / name
    if f.exists():
        p.update(json.loads(f.read_text()))
    return p


def _lerp_tab(ts, vs, t):
    return np.interp(t, ts, vs)


def profiles(p):
    y0, y1 = p["y_nose"], p["y_rear"]

    def T(y):
        return np.clip((y - y0) / (y1 - y0), 0, 1)

    def crest(y):
        t = T(y)
        vs = [p["z_nose"]] + list(p["crest"])
        # smooth: cubic through the points (Catmull-Rom on the uniform-ish t grid)
        return _cr(T_CREST, vs, float(t))

    def bot(y):
        t = T(y)
        vs = [p["z_nose"]] + list(p["bot"]) + [p["bot"][-1]]
        return _cr(T_BOT, vs, float(t))

    def halfw(y):
        t = float(T(y))
        return p["w_max"] * float(np.interp(t, PLAN_T, PLAN_W)) ** p["w_pow"]

    def mid(y):
        c, b = crest(y), bot(y)
        return c - p["mid_frac"] * (c - b)

    def sill(y):
        # from the bubble's bottom at y_jaw up to sill_rear at the rear rim, through sill_mid half way
        yj, yr = p["y_jaw"], y1
        if y <= yj:
            return -9.0
        s = (y - yj) / (yr - yj)
        zb = bot(yj)
        return _cr([0, 0.5, 1.0], [zb, p["sill_mid"], p["sill_rear"]], s)
    return crest, bot, mid, halfw, sill


def _cr(ts, vs, t):
    """Catmull-Rom through (ts, vs) at t (ts ascending), clamped"""
    ts, vs = list(ts), list(vs)
    if t <= ts[0]:
        return vs[0]
    if t >= ts[-1]:
        return vs[-1]
    i = max(0, min(len(ts) - 2, int(np.searchsorted(ts, t) - 1)))
    t0, t1 = ts[i], ts[i + 1]
    v0, v1 = vs[i], vs[i + 1]
    vm = vs[i - 1] if i > 0 else v0 - (v1 - v0)
    tm = ts[i - 1] if i > 0 else t0 - (t1 - t0)
    vp = vs[i + 2] if i + 2 < len(vs) else v1 + (v1 - v0)
    tp = ts[i + 2] if i + 2 < len(ts) else t1 + (t1 - t0)
    m0 = (v1 - vm) / (t1 - tm)
    m1 = (vp - v0) / (tp - t0)
    h = t1 - t0
    s = (t - t0) / h
    return ((2 * s**3 - 3 * s**2 + 1) * v0 + (s**3 - 2 * s**2 + s) * h * m0 + (-2 * s**3 + 3 * s**2) * v1
            + (s**3 - s**2) * h * m1)


def ring_point(p, y, a, crest, bot, mid, halfw):
    """point on the bubble at station y and angle a from the crest (a > 0 to port)"""
    zc, zb, zm, w = crest(y), bot(y), mid(y), halfw(y)
    s, c = math.sin(a), math.cos(a)
    x = w * math.copysign(abs(s) ** p["ex_side"], s)
    if c >= 0:
        z = zm + (zc - zm) * c ** p["ex_top"]
    else:
        z = zm - (zm - zb) * (-c) ** p["ex_bot"]
    return x, z


def angle_of_z(p, y, z, crest, bot, mid):
    """angle from the crest at which the ring reaches height z (on either flank), pi if below the bottom"""
    zc, zb, zm = crest(y), bot(y), mid(y)
    if z >= zc:
        return 0.0
    if z >= zm:
        return math.acos(((z - zm) / (zc - zm)) ** (1 / p["ex_top"]))
    if z <= zb:
        return math.pi
    return math.pi - math.acos(((zm - z) / (zm - zb)) ** (1 / p["ex_bot"]))


def canopy(p, n_st=34, n_around=56):
    crest, bot, mid, halfw, sill = profiles(p)
    y0, y1 = p["y_nose"], p["y_rear"]
    ys = y0 + (y1 - y0) * (1 - np.cos(np.linspace(0, math.pi, n_st + 1)[1:])) / 2
    rows, amax = [], []
    for y in ys:
        a = min(math.pi, angle_of_z(p, y, sill(y), crest, bot, mid))
        amax.append(a)
        row = [ring_point(p, y, ang, crest, bot, mid, halfw) for ang in np.linspace(a, -a, n_around)]
        rows.append([(x, y, z) for x, z in row])
    V = [(0.0, y0, p["z_nose"])] + [v for r in rows for v in r]
    V = np.array(V)
    F = [(0, 1 + j + 1, 1 + j) for j in range(n_around - 1)]
    for i in range(n_st - 1):
        for j in range(n_around - 1):
            a, b = 1 + i * n_around + j, 1 + i * n_around + j + 1
            c, d = 1 + (i + 1) * n_around + j + 1, 1 + (i + 1) * n_around + j
            F.append((a, d, c, b))
    # the frame: a band along the open edge: the port sill from the jaw's front to the rear, over the rear rim, the
    # starboard sill forward again (rows run from port, +a, to starboard, -a)
    opened = [i for i in range(n_st) if amax[i] < math.pi - 1e-6]
    path = [rows[i][0] for i in opened] + rows[-1][1:-1] + [rows[i][-1] for i in reversed(opened)]
    return V, F, np.array(path), dict(ys=ys, amax=amax)


def frame_band(path, p, width=0.07, depth=0.05, n_prof=6, rear_width=0.16, y_rear=None):
    """a rounded bar swept along path (an open polyline round the glass's edge), standing outward of the glass; round the rear
    rim (points at y_rear) it is rear_width wide and lies forward of the rim, over the glass's back (the picture's dark ring)"""
    P = np.asarray(path)
    m = len(P)
    y_rear = p["y_rear"] if y_rear is None else y_rear
    rear = np.abs(P[:, 1] - y_rear) < 1e-6
    k_rear = np.convolve(rear.astype(float), np.ones(5) / 5, mode="same")
    W = width + (rear_width - width) * k_rear
    shift = -0.42 * (rear_width - 0.0) * k_rear
    T = np.gradient(P, axis=0)
    T /= np.linalg.norm(T, axis=1, keepdims=True) + 1e-12
    # outward: away from the canopy's long axis (x = 0, z = mid)
    out = P.copy()
    out[:, 1] = 0
    out[:, 0] = P[:, 0]
    out[:, 2] = P[:, 2] - 0.2
    out -= (out * T).sum(1, keepdims=True) * T
    out /= np.linalg.norm(out, axis=1, keepdims=True) + 1e-12
    side = np.cross(T, out)
    prof = [(-0.5, 0.0), (-0.5, 0.6), (-0.3, 1.0), (0.3, 1.0), (0.5, 0.6), (0.5, 0.0)]
    V = []
    for k in range(m):
        for s, o in prof:
            V.append(P[k] + side[k] * s * W[k] + out[k] * (o * depth - 0.012) + np.array([0, shift[k], 0]))
    V = np.array(V)
    np_ = len(prof)
    F = []
    for k in range(m - 1):
        for j in range(np_ - 1):
            a, b = k * np_ + j, k * np_ + j + 1
            F.append((a, b, (k + 1) * np_ + j + 1, (k + 1) * np_ + j))
    return V, F