"""Door robot arm: the 3D fit for the scale-up and the door-centre attach (2026-10-18; plain Python + numpy).

The arm's own links (build/arm_points.npz, surface points per link, from this model's file) are posed with the rig's
inverse kinematics and checked against the real surroundings: everything in the ship round the pockets and the moon
pool (build/obstacles_ship.npz: deck, hull and lining, basin channel, aperture rim, lip and coaming, lamps, end walls,
pier and lower-deck walls; read from the ship with instances and linked objects) on a 5 cm grid of distances, and the
two doors (moved with their own animation; their inner face from door_centre.json, their outer face from the mesh
section). The design is a set of numbers (scale, rail line, link lengths, the J3 cylinder's stroke window, the tools'
reach), the carriage schedule is searched (dynamic programme along each door's slide).

    python build/dra_fit3d.py grid          build the obstacle grid (cached in build/obstacle_grid.npz)
    python build/dra_fit3d.py search        the search table
    python build/dra_fit3d.py check <json>  every quarter frame of both doors for one design -> build/fit3d_report.json
"""
import os
import sys
import json
import math
import copy
import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
ROOT = os.path.dirname(HERE)
BASE = json.load(open(os.path.join(ROOT, "design_base.json")))
DC = json.load(open(os.path.join(HERE, "door_centre.json")))
DP = json.load(open(os.path.join(HERE, "door_path.json")))
H = 0.05
ORIG = np.array([-10.2, -9.7, -11.2])
KMAX = 12
DECK = BASE["scene"]["deck_underside_z"]


# ------------------------------------------------------------------ the obstacle grid
def build_grid():
    d = np.load(os.path.join(HERE, "obstacles_ship.npz"))
    P = d["points"].astype(np.float64)
    ids = d["ids"]
    names = list(d["names"])
    dims = np.ceil((np.array([16.2, 2.2, -3.3]) - ORIG) / H).astype(int) + 1
    ix = np.floor((P - ORIG) / H).astype(int)
    ok = np.all((ix >= 0) & (ix < dims), axis=1)
    ix, ids = ix[ok], ids[ok]
    dist = np.full(dims, 255, np.uint8)
    who = np.full(dims, -1, np.int16)
    dist[ix[:, 0], ix[:, 1], ix[:, 2]] = 0
    who[ix[:, 0], ix[:, 1], ix[:, 2]] = ids
    for k in range(1, KMAX + 1):
        cur = dist <= k - 1
        new = np.zeros_like(cur)
        wnew = np.full(dims, -1, np.int16)
        for dx in (-1, 0, 1):
            for dy in (-1, 0, 1):
                for dz in (-1, 0, 1):
                    if dx == dy == dz == 0:
                        continue
                    src = [slice(max(0, -s), dims[i] - max(0, s)) for i, s in enumerate((dx, dy, dz))]
                    dst = [slice(max(0, s), dims[i] - max(0, -s)) for i, s in enumerate((dx, dy, dz))]
                    m = cur[tuple(src)] & (dist[tuple(dst)] == 255)
                    new[tuple(dst)] |= m
                    w = wnew[tuple(dst)]
                    w[m & (w < 0)] = who[tuple(src)][m & (w < 0)]
        dist[new] = k
        who[new] = wnew[new]
    np.savez_compressed(os.path.join(HERE, "obstacle_grid.npz"), dist=dist, who=who, names=np.array(names))
    print("grid", dims, "occupied", int((dist == 0).sum()))


_G = {}


def grid():
    if not _G:
        d = np.load(os.path.join(HERE, "obstacle_grid.npz"))
        _G.update(dist=d["dist"], who=d["who"], names=list(d["names"]))
    return _G


def lookup(P):
    """grid distance level (5 cm steps, Chebyshev) and obstacle index for world points P (n, 3)"""
    g = grid()
    ix = np.floor((P - ORIG) / H).astype(int)
    dims = np.array(g["dist"].shape)
    ok = np.all((ix >= 0) & (ix < dims), axis=1)
    lev = np.full(len(P), 255, np.int32)
    who = np.full(len(P), -1, np.int32)
    lev[ok] = g["dist"][ix[ok, 0], ix[ok, 1], ix[ok, 2]]
    who[ok] = g["who"][ix[ok, 0], ix[ok, 1], ix[ok, 2]]
    return lev, who


# ------------------------------------------------------------------ the arm
_A = {}


SUB = int(os.environ.get("DRA_FIT_SUB", "1500"))


def arm(sub=None):
    sub = SUB if sub is None else sub
    if not _A:
        d = np.load(os.path.join(HERE, "arm_points.npz"))
        rng = np.random.default_rng(1)
        for k in d.files:
            if k.startswith("pts__"):
                nm = k[5:].replace("_", " ")
                P = d[k].astype(np.float64)
                if sub and len(P) > sub:
                    P = P[rng.choice(len(P), sub, replace=False)]
                _A[nm] = dict(P=P, joint=d["joint__" + k[5:]], parent=str(d["parent__" + k[5:]]))
    return _A


def Rx(a):
    c, s = math.cos(a), math.sin(a)
    return np.array([[1, 0, 0, 0], [0, c, -s, 0], [0, s, c, 0], [0, 0, 0, 1.0]])


def Ry(a):
    c, s = math.cos(a), math.sin(a)
    return np.array([[c, 0, s, 0], [0, 1, 0, 0], [-s, 0, c, 0], [0, 0, 0, 1.0]])


def Rz(a):
    c, s = math.cos(a), math.sin(a)
    return np.array([[c, -s, 0, 0], [s, c, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1.0]])


def T(v):
    M = np.eye(4)
    M[:3, 3] = v
    return M


GRIP = {"proximal": math.radians(14.7), "distal": math.radians(41.0)}
FINGER_CH = {"top": (1, 1.0), "bottom": (1, -1.0), "left": (2, -1.0), "right": (2, 1.0)}


def joints_of(Dz):
    """scaled joint points of the links (arm frame), with the link length changes"""
    A = arm()
    S = Dz["S"]
    J = {nm: np.array(a["joint"]) * S for nm, a in A.items()}
    # link length changes: lower arm (J2 -> J3), upper arm (J3 -> J4), forearm (J4 -> WC)
    j2, j3, j4, wc = J["lower arm"], J["upper arm"], J["forearm"], J["wrist head"]
    j3n = j2 + (j3 - j2) * Dz.get("k_lower", 1.0)
    j4n = j3n + (j4 - j3) * Dz.get("k_upper", 1.0)
    wcn = j4n + (wc - j4) * Dz.get("k_fore", 1.0)
    shift = {"lower arm": np.zeros(3), "upper arm": j3n - j3, "forearm": j4n - j4, "wrist": wcn - wc}
    out = {}
    for nm, p in J.items():
        if nm in ("lower arm", "J3 cylinder", "J3 piston", "turret", "carriage", "carriage pinion"):
            out[nm] = p + (shift["lower arm"] if nm == "lower arm" else 0)
        elif nm == "upper arm":
            out[nm] = j3n
        elif nm == "forearm":
            out[nm] = j4n
        else:                                   # wrist head, flange, tools, fingers move with the wrist
            out[nm] = p + shift["wrist"]
    # the knob on the upper arm moves with J3 (the cylinder's far end): keep it relative to J3; Dz["knob"] = (f, rot deg)
    # moves it toward J3 by f and turns it about J3 in the arm's plane (x-z)
    k = J["J3 piston"] - j3
    f, rot = Dz.get("knob", (1.0, 0.0))
    a = math.radians(rot)
    ka, kb = k[0] * f, k[2] * f
    out["J3 piston"] = j3n + np.array([ka * math.cos(a) - kb * math.sin(a), k[1], ka * math.sin(a) + kb * math.cos(a)])
    return out


def kin(Dz):
    """IK constants of a design (dra_ik's names)"""
    J = joints_of(Dz)
    j2, j3, wc = J["lower arm"], J["upper arm"], J["wrist head"]
    L2v = j3 - j2
    L3v = wc - j3
    K = dict(j2_r=-j2[0], j2_h=j2[2], L2=math.hypot(L2v[0], L2v[2]), phi_l0=math.atan2(L2v[2], L2v[0]),
             L3=math.hypot(L3v[0], L3v[2]), beta=math.atan2(L3v[2], L3v[0]))
    # the forearm direction at rest is along A's x: beta = the angle of J3 -> WC above the forearm axis
    K["beta"] = math.atan2(L3v[2], L3v[0])
    return K, J


def ik(Dz, key, ver, loc, c):
    """joint values for door location loc and carriage c; None when out of reach"""
    K, J = kin(Dz)
    dr = Dz["doors"][key]
    wc = dr["wc_local"][ver]
    yb, zb = Dz["rail_y"], Dz["root_z"]
    dx, dy, dz = loc
    X, Y, Hh0 = dx + wc[0] - c, (yb - wc[1]) - dy, (zb - wc[2]) - dz
    q1 = math.atan2(Y, X)
    R = math.hypot(X, Y) + K["j2_r"]
    Hh = Hh0 - K["j2_h"]
    L2, L3, beta, phi_l0 = K["L2"], K["L3"], K["beta"], K["phi_l0"]
    Dd = math.hypot(R, Hh)
    cc = (L2 * L2 - L3 * L3 + Dd * Dd) / (2 * L2 * Dd)
    if abs(cc) >= 0.998:
        return None
    phi_l = math.atan2(Hh, R) + Dz.get("branch", 1) * math.acos(cc)
    phi_f = math.atan2(Hh - L2 * math.sin(phi_l), R - L2 * math.cos(phi_l)) - beta
    q2, q3 = phi_l - phi_l0, phi_f - phi_l + phi_l0
    xt, yt, zt = dr["tool_axes_C"]["x_t"], dr["tool_axes_C"]["y_t"], dr["tool_axes_C"]["z_t"]
    P = lambda v: math.cos(q1) * v[0] + math.sin(q1) * v[1]
    M11 = math.cos(phi_f) * P(xt) + math.sin(phi_f) * xt[2]
    M21 = -math.sin(q1) * xt[0] + math.cos(q1) * xt[1]
    M31 = math.cos(phi_f) * xt[2] - math.sin(phi_f) * P(xt)
    M12 = math.cos(phi_f) * P(yt) + math.sin(phi_f) * yt[2]
    M13 = math.cos(phi_f) * P(zt) + math.sin(phi_f) * zt[2]
    q5 = -math.acos(max(-1.0, min(1.0, M11)))
    q4 = math.atan2(-M21, M31)
    q6 = math.atan2(-M12, -M13)
    # the J3 cylinder: knob (on the upper arm) seen from the trunnion (on the lower arm), lower-arm frame
    tr = J["J3 cylinder"] - J["lower arm"]
    kn = J["J3 piston"] - J["upper arm"]          # the knob relative to J3 at rest
    j3l = J["upper arm"] - J["lower arm"]
    kx = (j3l[0] - tr[0]) + kn[0] * math.cos(q3) - kn[2] * math.sin(q3)
    kz = (j3l[2] - tr[2]) + kn[0] * math.sin(q3) + kn[2] * math.cos(q3)
    L = math.hypot(kx, kz)
    d0 = (J["J3 piston"] - J["J3 cylinder"])
    cyl = math.atan2(kz, kx) - math.atan2(d0[2], d0[0])
    return dict(q1=q1, q2=q2, q3=q3, q4=q4, q5=q5, q6=q6, L=L, cyl=cyl, c=c, L0=float(np.linalg.norm(d0[[0, 2]])))


def pose_points(Dz, q, ver, which=None):
    """world points of every link of this version at joint values q -> {link: (n, 3)}"""
    A = arm()
    S = Dz["S"]
    J = joints_of(Dz)
    yb, zb = Dz["rail_y"], Dz["root_z"]
    M_root = T((0.0, yb, zb)) @ Rx(math.pi)
    rot = {"carriage": np.eye(4), "carriage pinion": np.eye(4), "turret": Rz(q["q1"]), "lower arm": Ry(-q["q2"]),
           "upper arm": Ry(-q["q3"]), "forearm": Rx(q["q4"]), "wrist head": Ry(q["q5"]), "flange": Rx(q["q6"]),
           "J3 cylinder": Ry(-q["cyl"])}
    world = {}
    order = ["carriage", "carriage pinion", "turret", "lower arm", "J3 cylinder", "J3 piston", "upper arm", "forearm",
             "wrist head", "flange", "tool A hub", "tool B"] + ["finger %s %s" % (w, p) for w in ("top", "left", "bottom", "right")
                                                              for p in ("proximal", "distal")]
    out = {}
    for nm in order:
        if nm not in A:
            continue
        if ver == "A" and nm == "tool B":
            continue
        if ver == "B" and (nm == "tool A hub" or nm.startswith("finger")):
            continue
        a = A[nm]
        par = a["parent"]
        pj = J[par] if par else np.zeros(3)
        loc = J[nm] - pj
        if nm == "carriage":
            loc = loc + np.array([q["c"], 0.0, 0.0])
        if nm == "J3 piston":
            d0 = J["J3 piston"] - J["J3 cylinder"]
            u = d0 / np.linalg.norm(d0)
            loc = d0 + u * (q["L"] - q["L0"])
        R = rot.get(nm, np.eye(4))
        if nm.startswith("finger"):
            w, part = nm.split()[1], nm.split()[2]
            ch, sg = FINGER_CH[w]
            ang = sg * GRIP[part]
            R = (Rx, Ry, Rz)[ch](ang)
        Mp = world[par] if par else M_root
        M = Mp @ T(loc) @ R
        world[nm] = M
        if which and nm not in which:
            continue
        P = a["P"] * S
        if nm == "tool B" and Dz.get("plate_B"):
            P = tool_b_points(Dz)
        out[nm] = P @ M[:3, :3].T + M[:3, 3]
    return out, world


def tool_b_points(Dz):
    """version B's tool in its own frame (origin WC): the neck cone and the plate sized to the door (plate_B square)"""
    S = Dz["S"]
    neck = Dz.get("neck_B", 0.13 * S)
    pts = []
    x0 = 0.132 * S
    for k in range(13):
        x = x0 + neck * k / 12
        r = (0.077 + (0.131 - 0.077) * k / 12) * S
        for a in np.linspace(0, 2 * math.pi, 24, endpoint=False):
            pts.append((x, r * math.cos(a), r * math.sin(a)))
    s = Dz["plate_B"] / 2
    t = 0.044 * S
    xp = x0 + neck
    g = np.linspace(-s, s, 17)
    for u in g:
        for v in g:
            pts.append((xp, u, v))
            pts.append((xp + t, u, v))
    return np.array(pts)


# ------------------------------------------------------------------ the doors
def door_section(key):
    """inner and outer faces of the door's section (local y, z) as polylines; inside = between them"""
    return DC[key]["inner_curve_yz"]


def door_gap(Dz, key, loc, P):
    """signed distance (m) of points P to the door at location loc: + on the arm's side of its inner face; only points
    over the door's length count (others +9)"""
    c = np.array(DC[key]["inner_curve_yz"])
    lx0, lx1 = DC[key]["length_x"]
    cy, cz = c[:, 0] + loc[1], c[:, 1] + loc[2]
    a = np.stack([cy[:-1], cz[:-1]], 1)
    b = np.stack([cy[1:], cz[1:]], 1)
    e = b - a
    L2 = (e * e).sum(1)
    q = P[:, [1, 2]]
    t = np.clip(((q[:, None, :] - a[None]) * e[None]).sum(2) / L2[None], 0, 1)
    proj = a[None] + t[..., None] * e[None]
    dv = q[:, None, :] - proj
    dd = np.sqrt((dv * dv).sum(2))
    k = dd.argmin(1)
    n = np.stack([-e[:, 1], e[:, 0]], 1) / np.sqrt(L2)[:, None]
    nrm = np.array(DC[key]["normal_local"][1:])
    if (n[len(n) // 2] @ nrm) < 0:
        n = -n
    s = np.sign((dv[np.arange(len(q)), k] * n[k]).sum(1))
    g = dd[np.arange(len(q)), k] * np.where(s == 0, 1, s)
    over = (P[:, 0] >= loc[0] + lx0) & (P[:, 0] <= loc[0] + lx1)
    return np.where(over, g, 9.0)


# ------------------------------------------------------------------ a design and its evaluation
def make_design(S=1.5, rail_y=-5.795, plate_B=1.2, k_lower=1.0, k_upper=1.0, k_fore=1.0, neck_B=None, dt_A=None,
                stroke=(-25.0, 47.0), branch=1, attach="centre", knob=(1.0, 0.0), j1=(-170.0, 170.0)):
    Dz = dict(S=S, rail_y=rail_y, plate_B=plate_B, k_lower=k_lower, k_upper=k_upper, k_fore=k_fore, branch=branch,
              stroke=stroke, knob=tuple(knob), j1=tuple(j1))
    Dz["root_z"] = DECK + S * (BASE["carriage"]["underside_z"] - DECK)
    pad_top = 0.012 + max(v[1] for k, v in DC["AFT"]["sag_under_plate_m"].items() if float(k) <= (plate_B or 0.6) + 0.11)
    Dz["neck_B"] = neck_B if neck_B is not None else 0.13 * S
    flange_B = 0.132 * S
    dB = flange_B + Dz["neck_B"] + 0.044 * S + pad_top
    dA = dt_A if dt_A is not None else 0.63 * S
    Dz["d_t"] = {"A": dA, "B": dB}
    Dz["doors"] = {}
    for key in ("AFT", "FORE"):
        att, n = DC[key]["attach_local"], DC[key]["normal_local"]
        xt = [-n[0], -n[1], -n[2]]
        yt = [-1.0, 0.0, 0.0]
        zt = [xt[1] * yt[2] - xt[2] * yt[1], xt[2] * yt[0] - xt[0] * yt[2], xt[0] * yt[1] - xt[1] * yt[0]]
        C = lambda v: [v[0], -v[1] + 0.0, -v[2] + 0.0]
        Dz["doors"][key] = dict(wc_local={v: [att[i] + Dz["d_t"][v] * n[i] for i in range(3)] for v in ("A", "B")},
                                tool_axes_C=dict(x_t=C(xt), y_t=C(yt), z_t=C(zt)))
    return Dz


LIM = {"q2": (-95, 45), "q3": (-40, 55), "q4": (-170, 170), "q5": (-120, 120), "q6": (-170, 170)}


def evaluate(Dz, key, loc, c, vers=("A", "B"), detail=False):
    """(static gap level min, door gap min, joint margin min, piston ok margin) over both versions"""
    worst = dict(lev=99, door=9.0, joint=99.0, stroke=99.0, who="", link="", q3=[999.0, -999.0], jname="")
    for ver in vers:
        q = ik(Dz, key, ver, loc, c)
        if q is None:
            return dict(lev=-1, door=-9.0, joint=-99.0, stroke=-99.0, who="out of reach", link="", q3=[999.0, -999.0], jname="")
        q3d = math.degrees(q["q3"])
        worst["q3"] = [min(worst["q3"][0], q3d), max(worst["q3"][1], q3d)]
        for j, (lo, hi) in LIM.items():
            if j == "q3":
                if Dz.get("free_q3"):
                    continue
                lo, hi = Dz["stroke"]
            v = math.degrees(q[j])
            m = min(v - lo, hi - v)
            if m < worst["joint"]:
                worst["joint"] = m
                worst["jname"] = ver + j
        pts, _ = pose_points(Dz, q, ver)
        for nm, P in pts.items():
            if nm in ("carriage", "carriage pinion"):
                continue
            lev, who = lookup(P)
            i = int(lev.argmin())
            if lev[i] < worst["lev"]:
                worst.update(lev=int(lev[i]), who=grid()["names"][who[i]] if who[i] >= 0 else "", link=nm)
            if nm in ("tool B", "tool A hub", "turret", "J3 piston") or nm.startswith("finger"):
                continue
            for dk in (key,):
                dl = door_location(dk, loc, key)
                g = door_gap(Dz, dk, dl, P)
                m = float(g.min())
                if m < worst["door"]:
                    worst["door"] = m
                    worst["door_link"] = nm
    return worst


def door_location(dk, loc, key):
    """the other door moves the same way, mirrored about the mid plane (x 3.05)"""
    if dk == key:
        return loc
    return (2 * 3.05 - loc[0], loc[1], loc[2])


def score(w):
    return min((w["lev"] - 2) * H, w["door"] - 0.08, (w["joint"] - 3.0) / 50.0)


def schedule(Dz, key="AFT", n_c=15, bmax=1.5, verbose=False, step=8):
    rail_end = deck_edge(Dz["rail_y"], Dz["S"])
    S = Dz["S"]
    half = 0.429 * S
    if key == "AFT":
        cmax = rail_end - 0.04 - half
        cmin = -9.0348 + 0.05 + max(half, 0.53 * S)
    else:
        cmin = 2 * 3.05 - rail_end + 0.04 + half
        cmax = 15.1348 - 0.05 - max(half, 0.53 * S)
    cs = np.linspace(cmin, cmax, n_c)
    pts = [p for p in DP[key][::step]]
    xs_closed = pts[0][1]
    # stages: phase 1 (the door at its closed x, pulled in), then the slide (frames 12-48), as door x values
    ph1 = [p for p in pts if abs(p[1] - xs_closed) < 1e-4]
    slide = sorted({round(p[1], 4) for p in pts if abs(p[2] + 4.88) < 1e-4}, reverse=(key == "AFT"))
    dz = pts[0][3]
    stage = []
    for x in slide:
        stage.append([score(evaluate(Dz, key, (x, -4.88, dz), c)) for c in cs])
    p1 = [min(score(evaluate(Dz, key, (p[1], p[2], p[3]), c)) for p in ph1[::2]) for c in cs]
    stage[0] = [min(a, b) for a, b in zip(stage[0], p1)]
    best = [(stage[0][j], [j]) for j in range(n_c)]
    for k in range(1, len(slide)):
        ddx = abs(slide[k] - slide[k - 1])
        nb = []
        for j in range(n_c):
            cand = None
            for i in range(n_c):
                if abs(cs[i] - cs[j]) > bmax * ddx + 1e-9:
                    continue
                v = min(best[i][0], stage[k][j])
                if cand is None or v > cand[0]:
                    cand = (v, best[i][1] + [j])
            nb.append(cand or (-99.0, [j]))
        best = nb
    v, path = max(best, key=lambda b: b[0])
    sched = [[round(x, 3), round(float(cs[j]), 3)] for x, j in zip(slide, path)]
    q3r = [999.0, -999.0]
    for (x, c) in sched:
        w = evaluate(Dz, key, (x, -4.88, dz), c)
        q3r = [min(q3r[0], w["q3"][0]), max(q3r[1], w["q3"][1])]
    for p in ph1[::2]:
        w = evaluate(Dz, key, (p[1], p[2], p[3]), sched[0][1])
        q3r = [min(q3r[0], w["q3"][0]), max(q3r[1], w["q3"][1])]
    return dict(score=round(float(v), 3), schedule=sched, c_limits=[round(cmin, 3), round(cmax, 3)],
                q3_range=[round(q3r[0], 1), round(q3r[1], 1)])


_DE = {}


def deck_edge(y, S):
    """the deck's edge at the moon pool (AFT x) over the rail's width at line y (from pocket_meshes.json)"""
    key = (round(y, 3), round(S, 3))
    if key in _DE:
        return _DE[key]
    d = json.load(open(os.path.join(HERE, "pocket_meshes.json")))["Moonpool | deck around open well"]
    V = np.array(d["v"])
    Tt = np.array(d["t"])
    A, B, Cc = V[Tt[:, 0]], V[Tt[:, 1]], V[Tt[:, 2]]
    N = np.cross(B - A, Cc - A)
    N /= np.linalg.norm(N, axis=1)[:, None] + 1e-12
    down = N[:, 2] < -0.9
    A, B, Cc = A[down], B[down], Cc[down]
    half_w = 0.229 * S + 0.03
    best = 9.0
    for yy in np.linspace(y - half_w, y + half_w, 7):
        x = -9.0
        while x < 3.0:
            p = np.array([x + 0.025, yy])
            v0, v1 = Cc[:, :2] - A[:, :2], B[:, :2] - A[:, :2]
            v2 = p[None] - A[:, :2]
            d00 = (v0 * v0).sum(1); d01 = (v0 * v1).sum(1); d11 = (v1 * v1).sum(1)
            d20 = (v2 * v0).sum(1); d21 = (v2 * v1).sum(1)
            den = d00 * d11 - d01 * d01
            den = np.where(abs(den) < 1e-12, 1e-12, den)
            u = (d11 * d20 - d01 * d21) / den
            w = (d00 * d21 - d01 * d20) / den
            if not np.any((u >= 0) & (w >= 0) & (u + w <= 1)):
                break
            x += 0.05
        best = min(best, x)
    _DE[key] = best
    return best


if __name__ == "__main__":
    cmd = sys.argv[1] if len(sys.argv) > 1 else "search"
    if cmd == "grid":
        build_grid()
    elif cmd == "probe":
        Dz = make_design(S=float(sys.argv[2]), rail_y=float(sys.argv[3]))
        print(json.dumps(schedule(Dz), indent=0)[:1500])


def explain(Dz, r, key="AFT", step=8):
    """the components along a schedule (from schedule())"""
    pts = [p for p in DP[key][::step]]
    dz = pts[0][3]
    rows = []
    for (x, c) in r["schedule"]:
        w = evaluate(Dz, key, (x, -4.88, dz), c)
        rows.append((x, c, w))
    xs_closed = pts[0][1]
    c0 = r["schedule"][0][1]
    for p in pts:
        if abs(p[1] - xs_closed) < 1e-4 and p[2] < -4.881:
            rows.insert(0, (p[1], c0, evaluate(Dz, key, (p[1], p[2], p[3]), c0)))
    for x, c, w in rows:
        print("   door x %6.2f c %6.2f | lev %3d (%s / %s) door %6.2f (%s) joints %6.1f" % (
            x, c, w["lev"], w.get("link", ""), str(w.get("who", ""))[:34], w["door"], w.get("door_link", ""), w["joint"]))
    return rows