"""Procedural brass orrery: one Python file, numpy only, writes a glTF binary (.glb).

  python orrery.py [out.glb]

Everything is built from math, no assets:
- involute gear teeth (20° pressure angle) placed so the gear train really meshes;
- a walnut pedestal and fluted columns turned on a "lathe" (surface of revolution, flutes as a radius modulation);
- planets coloured per vertex by procedural noise (continents, Jupiter bands and spot, Saturn ring divisions);
- a zodiac ring with 72 ticks and 12 enamel plaques with star emblems, an armillary around the Sun, a hand crank.
Moving parts are separate nodes with extras {"spin": rad/s, "axis": "x"|"y"}, so a viewer can animate them.
"""
import json
import math
import os
import struct
import sys

import numpy as np

TAU = 2 * math.pi


# ---------------- primitives: each returns (V, F) ----------------

def rot(axis, ang):
    x, y, z = np.asarray(axis, float) / np.linalg.norm(axis)
    c, s, C = math.cos(ang), math.sin(ang), 1 - math.cos(ang)
    return np.array([[c + x * x * C, x * y * C - z * s, x * z * C + y * s],
                     [y * x * C + z * s, c + y * y * C, y * z * C - x * s],
                     [z * x * C - y * s, z * y * C + x * s, c + z * z * C]])


def tf(part, R=None, t=(0, 0, 0), s=1.0):
    V, F = part
    V = V * np.asarray(s, float)
    if R is not None:
        V = V @ R.T
    return V + np.asarray(t, float), F


def merge(*parts):
    Vs, Fs, o = [], [], 0
    for V, F in parts:
        Vs.append(V)
        Fs.append(F + o)
        o += len(V)
    return np.concatenate(Vs), np.concatenate(Fs)


def strip(A, B):
    """Closed quad strip between two loops; faces outward for A below B with angle growing x->z."""
    n = len(A)
    j = np.arange(n)
    a, b, c, d = j, (j + 1) % n, n + (j + 1) % n, n + j
    return np.concatenate([A, B]), np.concatenate([np.stack([a, c, b], 1), np.stack([a, d, c], 1)])


def loop(r, th, y):
    r = np.broadcast_to(np.asarray(r, float), np.shape(th))
    return np.stack([r * np.cos(th), np.full(len(th), y), r * np.sin(th)], 1)


def prism(r_out, th, r_in, y0, y1):
    """Extrude the ring between an outline (r_out, th) and radius r_in (0 = solid) from y0 to y1."""
    ob, ot, ib, it = loop(r_out, th, y0), loop(r_out, th, y1), loop(r_in, th, y0), loop(r_in, th, y1)
    return merge(strip(ob, ot), strip(it, ib), strip(ot, it), strip(ib, ob))


def lathe(prof, seg=96, closed=False, rmod=None):
    """Surface of revolution of a (r, y) profile around the y axis; rmod(y, theta) scales the radius."""
    p = np.asarray(prof, float)
    n = len(p)
    th = np.arange(seg) * TAU / seg
    r = np.repeat(p[:, :1], seg, 1)
    if rmod is not None:
        r = r * rmod(p[:, 1:2], th[None, :])
    y = np.repeat(p[:, 1:2], seg, 1)
    V = np.stack([r * np.cos(th), y, r * np.sin(th)], -1).reshape(-1, 3)
    i = np.arange(n if closed else n - 1)[:, None]
    j = np.arange(seg)[None, :]
    a, b, c, d = i * seg + j, i * seg + (j + 1) % seg, (i + 1) % n * seg + (j + 1) % seg, (i + 1) % n * seg + j
    return V, np.concatenate([np.stack([a, c, b], -1).reshape(-1, 3), np.stack([a, d, c], -1).reshape(-1, 3)])


def hard(*pts):
    """Profile corners: each point twice, so smooth normals keep a sharp edge there."""
    return [p for p in pts for _ in (0, 1)]


def ogee(r0, y0, r1, y1, n=10):
    t = np.linspace(0, 1, n)
    s = 3 * t ** 2 - 2 * t ** 3
    return list(zip(r0 + (r1 - r0) * s, y0 + (y1 - y0) * t))


def sphere(r, seg=48, rings=24):
    phi = np.linspace(0, math.pi, rings + 1)
    return lathe(list(zip(r * np.sin(phi), -r * np.cos(phi))), seg)


def torus(R, r, seg=96, rseg=14):
    phi = np.arange(rseg) * TAU / rseg
    return lathe(list(zip(R + r * np.cos(phi), r * np.sin(phi))), seg, closed=True)


def washer(r_in, r_out, y0, y1, seg=64):
    return lathe(hard((r_in, y0), (r_out, y0), (r_out, y1), (r_in, y1)), seg, closed=True)


def box(sx, sy, sz):
    x, y, z = sx / 2, sy / 2, sz / 2
    V = np.array([[-x, -y, -z], [x, -y, -z], [x, y, -z], [-x, y, -z], [-x, -y, z], [x, -y, z], [x, y, z], [-x, y, z]])
    F = np.array([[0, 2, 1], [0, 3, 2], [4, 5, 6], [4, 6, 7], [0, 1, 5], [0, 5, 4],
                  [3, 7, 6], [3, 6, 2], [0, 4, 7], [0, 7, 3], [1, 2, 6], [1, 6, 5]])
    return V, F


def align_y(d):
    """Rotation taking +y to direction d."""
    d = np.asarray(d, float) / np.linalg.norm(d)
    ax = np.cross([0, 1, 0], d)
    if np.linalg.norm(ax) < 1e-9:
        return np.eye(3) if d[1] > 0 else rot((1, 0, 0), math.pi)
    return rot(ax, math.acos(np.clip(d[1], -1, 1)))


def rod(p0, p1, r, seg=24):
    p0, p1 = np.asarray(p0, float), np.asarray(p1, float)
    L = np.linalg.norm(p1 - p0)
    return tf(lathe(hard((0, 0), (r, 0), (r, L), (0, L)), seg), align_y(p1 - p0), p0)


# ---------------- gears ----------------

def gear_outline(z, m, phase=0.0, pa=math.radians(20), k=7):
    """Involute spur gear outline as (radius, angle) arrays; a tooth is centred at phase + i*2pi/z."""
    rp = m * z / 2
    rb, ra, rf = rp * math.cos(pa), rp + m, rp - 1.25 * m
    inv = lambda a: math.tan(a) - a
    psi_p = math.pi / (2 * z) + inv(pa)  # half tooth angle at the base circle
    psi = lambda r: psi_p - inv(math.acos(min(1.0, rb / max(r, rb))))
    rs = np.linspace(rf, ra, k)
    R, T = [], []
    for i in range(z):
        c = phase + i * TAU / z
        R += list(rs); T += [c - psi(r) for r in rs]
        top = np.linspace(c - psi(ra), c + psi(ra), 4)[1:-1]
        R += [ra] * len(top); T += list(top)
        R += list(rs[::-1]); T += [c + psi(r) for r in rs[::-1]]
        root = np.linspace(c + psi(rf), c + TAU / z - psi(rf), 5)[1:-1]
        R += [rf] * len(root); T += list(root)
    return np.array(R), np.array(T), rf


def gear(z, m, th, bore, phase=0.0, spokes=0):
    R, T, rf = gear_outline(z, m, phase)
    y0, y1 = -th / 2, th / 2
    hub_r = max(bore * 2.0, bore + 2.5 * m)
    hub = washer(bore, hub_r, y0 - th * 0.35, y1 + th * 0.35, 48)
    if not spokes:
        return merge(prism(R, T, bore, y0, y1), hub)
    rim = rf - 2.4 * m
    parts = [prism(R, T, rim, y0, y1), hub, torus(rim, th * 0.22, 128, 8)]
    mid, L = (hub_r + rim) / 2, rim - hub_r + 2 * m
    for s in range(spokes):
        a = phase + (s + 0.5) * TAU / spokes
        parts.append(tf(box(L, th * 0.55, max(1.4 * m, rim * 0.09)), rot((0, 1, 0), -a),
                        (mid * math.cos(a), 0, mid * math.sin(a))))
    return merge(*parts)


def star(n, ro, ri, h):
    th = np.arange(2 * n) * math.pi / n
    return prism(np.where(np.arange(2 * n) % 2 == 0, ro, ri), th, 0.0, 0, h)


# ---------------- procedural colour ----------------

def lin(c):
    return [float(v) ** 2.2 for v in c]


def fbm(P, seed, octaves=5, freq=1.6):
    rng = np.random.default_rng(seed)
    v, amp = np.zeros(len(P)), 1.0
    for o in range(octaves):
        k = rng.normal(size=(8, 3)) * freq * 2 ** o
        v += amp * np.sin(P @ k.T + rng.uniform(0, TAU, 8)).mean(1)
        amp *= 0.55
    return v


def mix(a, b, t):
    t = np.clip(t, 0, 1)[:, None]
    return np.asarray(a) * (1 - t) + np.asarray(b) * t


def paint(kind, V):
    P = V / np.linalg.norm(V, axis=1, keepdims=True)
    y = P[:, 1]
    if kind == "earth":
        n = fbm(P, 3) + 0.15 * fbm(P * 3, 4)
        land = mix((0.18, 0.42, 0.14), (0.55, 0.45, 0.25), fbm(P, 9) * 2 + 0.5)
        sea = mix((0.02, 0.1, 0.35), (0.06, 0.3, 0.6), n * 4 + 0.8)
        c = np.where((n > 0.06)[:, None], land, sea)
        return np.where((np.abs(y) > 0.86 + 0.05 * fbm(P, 5))[:, None], (0.95, 0.97, 1.0), c)
    if kind == "jupiter":
        b = np.sin(y * 18 + 1.2 * fbm(P * 2, 7)) * 0.5 + 0.5
        c = mix((0.92, 0.84, 0.7), (0.66, 0.44, 0.28), b)
        spot = np.linalg.norm(P - np.array([0.75, -0.38, 0.54]) / np.linalg.norm([0.75, -0.38, 0.54]), axis=1)
        return mix(c, (0.78, 0.3, 0.18), (0.2 - spot) * 12)
    if kind == "saturn":
        return mix((0.95, 0.88, 0.66), (0.78, 0.66, 0.44), np.sin(y * 14 + 0.6 * fbm(P, 2)) * 0.5 + 0.5)
    if kind == "mars":
        c = mix((0.78, 0.36, 0.17), (0.42, 0.2, 0.12), fbm(P, 11) * 2.5 + 0.4)
        return np.where((np.abs(y) > 0.9)[:, None], (0.95, 0.93, 0.9), c)
    if kind == "venus":
        return mix((0.95, 0.86, 0.6), (0.82, 0.68, 0.42), fbm(P * 1.5 + y[:, None] * 3, 13) * 2 + 0.5)
    # mercury, moon: grey with craters-ish noise
    return mix((0.38, 0.36, 0.34), (0.75, 0.73, 0.7), fbm(P * 2, 17 if kind == "moon" else 19) * 2.5 + 0.5)


# ---------------- scene assembly ----------------

class Geo:
    """Parts grouped by material; flat parts get per-face normals, the rest smooth normals."""

    def __init__(self):
        self.items = {}

    def add(self, mat, part, color=None, flat=False):
        self.items.setdefault(mat, []).append((part, color, flat))
        return self

    def build(self):
        out = {}
        for mat, items in self.items.items():
            Vs, Ns, Fs, Cs, o = [], [], [], [], 0
            for (V, F), C, flat in items:
                if flat:
                    V = V[F].reshape(-1, 3)
                    C = None if C is None else np.asarray(C)[F].reshape(-1, 3)
                    F = np.arange(len(V)).reshape(-1, 3)
                fn = np.cross(V[F[:, 1]] - V[F[:, 0]], V[F[:, 2]] - V[F[:, 0]])
                N = np.zeros_like(V)
                for k in range(3):
                    np.add.at(N, F[:, k], fn)
                ln = np.linalg.norm(N, axis=1, keepdims=True)
                N = np.where(ln > 1e-12, N / np.where(ln > 1e-12, ln, 1), [0, 1, 0])
                Vs.append(V); Ns.append(N); Fs.append(F + o); o += len(V)
                Cs.append(None if C is None else np.asarray(C, float) ** 2.2)
            C = np.concatenate(Cs) if all(c is not None for c in Cs) else None
            out[mat] = (np.concatenate(Vs), np.concatenate(Ns), np.concatenate(Fs), C)
        return out


class GLB:
    def __init__(self):
        self.bin = bytearray()
        self.j = {"asset": {"version": "2.0", "generator": "orrery.py (numpy)"}, "scene": 0,
                  "scenes": [{"nodes": [0]}], "nodes": [], "meshes": [], "materials": [],
                  "accessors": [], "bufferViews": [], "buffers": []}
        self.mat = {}
        self.tris = 0

    def material(self, name, color, metal, rough, emissive=None):
        m = {"name": name, "doubleSided": True,
             "pbrMetallicRoughness": {"baseColorFactor": lin(color) + [1.0], "metallicFactor": metal,
                                      "roughnessFactor": rough}}
        if emissive:
            m["emissiveFactor"] = list(emissive)
        self.mat[name] = len(self.j["materials"])
        self.j["materials"].append(m)

    def _acc(self, arr, target, typ, comp, minmax=False):
        off = len(self.bin)
        self.bin += arr.tobytes()
        self.bin += b"\0" * (-len(self.bin) % 4)
        self.j["bufferViews"].append({"buffer": 0, "byteOffset": off, "byteLength": arr.nbytes, "target": target})
        a = {"bufferView": len(self.j["bufferViews"]) - 1, "componentType": comp, "count": len(arr), "type": typ}
        if minmax:
            a["min"], a["max"] = arr.min(0).tolist(), arr.max(0).tolist()
        self.j["accessors"].append(a)
        return len(self.j["accessors"]) - 1

    def mesh(self, name, geo):
        prims = []
        for mat, (V, N, F, C) in geo.build().items():
            att = {"POSITION": self._acc(V.astype(np.float32), 34962, "VEC3", 5126, True),
                   "NORMAL": self._acc(N.astype(np.float32), 34962, "VEC3", 5126)}
            if C is not None:
                att["COLOR_0"] = self._acc(C.astype(np.float32), 34962, "VEC3", 5126)
            prims.append({"attributes": att, "material": self.mat[mat],
                          "indices": self._acc(F.astype(np.uint32).ravel(), 34963, "SCALAR", 5125)})
            self.tris += len(F)
        self.j["meshes"].append({"name": name, "primitives": prims})
        return len(self.j["meshes"]) - 1

    def node(self, name, mesh=None, t=None, q=None, spin=None, axis="y", parent=0):
        n = {"name": name}
        if mesh is not None:
            n["mesh"] = mesh
        if t is not None:
            n["translation"] = [float(v) for v in t]
        if q is not None:
            n["rotation"] = [float(v) for v in q]
        if spin:
            n["extras"] = {"spin": spin, "axis": axis}
        self.j["nodes"].append(n)
        idx = len(self.j["nodes"]) - 1
        if parent is not None and idx != parent:
            self.j["nodes"][parent].setdefault("children", []).append(idx)
        return idx

    def save(self, path):
        self.j["buffers"] = [{"byteLength": len(self.bin)}]
        js = json.dumps(self.j, separators=(",", ":")).encode()
        js += b" " * (-len(js) % 4)
        total = 12 + 8 + len(js) + 8 + len(self.bin)
        with open(path, "wb") as f:
            f.write(struct.pack("<4sII", b"glTF", 2, total))
            f.write(struct.pack("<I4s", len(js), b"JSON") + js)
            f.write(struct.pack("<I4s", len(self.bin), b"BIN\0") + bytes(self.bin))


def quat(axis, ang):
    a = np.asarray(axis, float) / np.linalg.norm(axis)
    return [*(a * math.sin(ang / 2)), math.cos(ang / 2)]


def build(out):
    g = GLB()
    g.material("brass", (0.80, 0.62, 0.30), 1.0, 0.28)
    g.material("bronze", (0.62, 0.40, 0.22), 1.0, 0.38)
    g.material("steel", (0.66, 0.68, 0.72), 1.0, 0.22)
    g.material("walnut", (0.30, 0.16, 0.08), 0.0, 0.45)
    g.material("enamel", (0.06, 0.12, 0.40), 0.1, 0.15)
    g.material("ivory", (0.92, 0.88, 0.78), 0.0, 0.35)
    g.material("planet", (1, 1, 1), 0.0, 0.65)
    g.material("sun", (1.0, 0.72, 0.25), 0.3, 0.4, emissive=(1.0, 0.55, 0.12))
    root = g.node("Orrery", parent=None)

    # pedestal: walnut turned base on ball feet, brass top plate with bead, engraved rings and rivets
    ped = [(0.0, 0.0)] + hard((1.0, 0.0), (1.0, 0.035)) + ogee(1.0, 0.035, 0.92, 0.12)
    ped += [(0.90 + 0.022 * math.sin(a), 0.142 - 0.022 * math.cos(a)) for a in np.linspace(0, math.pi, 12)]
    ped += hard((0.88, 0.17)) + [(0.88, 0.25)] + ogee(0.88, 0.25, 0.975, 0.33)[1:] + hard((0.0, 0.33))
    base = Geo().add("walnut", lathe(ped, 192))
    for k in range(8):
        a = k * TAU / 8
        base.add("walnut", tf(sphere(0.07, 32, 16), None, (0.84 * math.cos(a), -0.02, 0.84 * math.sin(a)), (1, 0.55, 1)))
    base.add("brass", lathe([(0.0, 0.33)] + hard((0.975, 0.33), (0.975, 0.35)) + [(0.0, 0.35)], 192))
    base.add("brass", tf(torus(0.975, 0.011, 192, 12), None, (0, 0.34, 0)))
    for R in (0.50, 0.72, 0.905):
        base.add("bronze", tf(torus(R, 0.0035, 192, 8), None, (0, 0.35, 0)))
    for k in range(96):
        a = k * TAU / 96
        base.add("bronze", tf(sphere(0.008, 10, 6), None, (0.94 * math.cos(a), 0.35, 0.94 * math.sin(a)), (1, 0.6, 1)))
    base.add("brass", tf(torus(0.03, 0.007, 32, 10), rot((0, 0, 1), math.pi / 2), (0.885, 0.21, 0)))  # crank bezel
    # central shaft and the axles of the gear train (static)
    base.add("steel", lathe([(0, 0.35)] + hard((0.022, 0.35), (0.022, 1.14)) + [(0, 1.14)], 48))
    base.add("brass", washer(0.022, 0.05, 0.35, 0.37))
    g.node("pedestal", g.mesh("pedestal", base))

    # six fluted columns holding the zodiac ring (one mesh, six nodes)
    col = [(0, 0)] + hard((0.05, 0), (0.05, 0.02), (0.038, 0.02))
    col += [(0.038 + 0.012 * math.sin(a), 0.032 - 0.012 * math.cos(a)) for a in np.linspace(0, math.pi, 10)[1:]]
    col += hard((0.028, 0.048)) + [(0.028, 0.055), (0.026, 0.19), (0.024, 0.325), (0.024, 0.33)]
    col += ogee(0.024, 0.33, 0.048, 0.385)[1:] + hard((0.05, 0.385), (0.05, 0.40)) + [(0, 0.40)]
    flutes = lambda y, th: np.where((y > 0.054) & (y < 0.326), 1 - 0.16 * (0.5 + 0.5 * np.cos(16 * th)) ** 3, 1.0)
    col_mesh = g.mesh("column", Geo().add("brass", lathe(col, 128, rmod=flutes)))
    for k in range(6):
        a = k * TAU / 6
        g.node(f"column_{k}", col_mesh, (0.88 * math.cos(a), 0.35, 0.88 * math.sin(a)))

    # zodiac ring: 72 ticks, 12 enamel plaques, each with a star of 3..14 points
    zod = Geo().add("brass", washer(0.80, 0.96, 0.75, 0.77, 192))
    for R in (0.80, 0.96):
        zod.add("brass", tf(torus(R, 0.007, 192, 10), None, (0, 0.77, 0)))
    for k in range(72):
        a, L = k * TAU / 72, 0.045 if k % 6 == 0 else 0.022
        zod.add("bronze", tf(box(L, 0.004, 0.0035), rot((0, 1, 0), -a),
                             ((0.952 - L / 2) * math.cos(a), 0.772, (0.952 - L / 2) * math.sin(a))), flat=True)
    for k in range(12):
        a = (k + 0.5) * TAU / 12
        p = (0.865 * math.cos(a), 0.771, 0.865 * math.sin(a))
        zod.add("enamel", tf(box(0.06, 0.004, 0.11), rot((0, 1, 0), -a), p), flat=True)
        zod.add("brass", tf(star(k + 3, 0.022, 0.009, 0.004), rot((0, 1, 0), -a), (p[0], 0.773, p[2])), flat=True)
    g.node("zodiac_ring", g.mesh("zodiac_ring", zod))

    # gear train: sun gear A0 on the shaft, three idlers A1, each with a pinion B1 driving B2
    w0, mA, mB = 0.3, 0.012, 0.008
    g.node("gear_A0", g.mesh("gear_A0", Geo().add("brass", gear(72, mA, 0.025, 0.023, 0.0, spokes=6), flat=True)),
           (0, 0.42, 0), spin=w0)
    axles = Geo()
    # identical gears share one mesh; each node turns it to its own tooth phase (node rotation -phi = angle +phi)
    meshA1 = g.mesh("gear_A1", Geo().add("bronze", gear(30, mA, 0.022, 0.009, 0.0, spokes=5), flat=True))
    meshB1 = g.mesh("pinion_B1", Geo().add("steel", gear(14, mB, 0.03, 0.009, 0.0), flat=True))
    meshB2 = g.mesh("gear_B2", Geo().add("brass", gear(40, mB, 0.018, 0.009, 0.0, spokes=4), flat=True))
    for k in range(3):
        phi = math.radians(30 + 120 * k)
        d = (72 + 30) * mA / 2
        c1 = np.array([d * math.cos(phi), 0, d * math.sin(phi)])
        g.node(f"gear_A1_{k}", meshA1, c1 + (0, 0.42, 0), quat((0, 1, 0), -(phi + math.pi + math.pi / 30)),
               spin=-w0 * 72 / 30)
        dB = (14 + 40) * mB / 2
        c2 = c1 + dB * np.array([-math.sin(phi), 0, math.cos(phi)])
        psi = math.atan2(c2[2] - c1[2], c2[0] - c1[0])
        g.node(f"pinion_B1_{k}", meshB1, c1 + (0, 0.49, 0), quat((0, 1, 0), -psi), spin=-w0 * 72 / 30)
        g.node(f"gear_B2_{k}", meshB2, c2 + (0, 0.49, 0), quat((0, 1, 0), -(psi + math.pi + math.pi / 40)),
               spin=w0 * 72 / 30 * 14 / 40)
        for c, top in ((c1, 0.52), (c2, 0.52)):
            axles.add("steel", rod(c + (0, 0.35, 0), c + (0, top, 0), 0.008, 20))
            axles.add("brass", tf(lathe([(0, 0)] + ogee(0.016, 0, 0.004, 0.02) + [(0.009, 0.026), (0, 0.032)], 24),
                                  None, c + (0, top, 0)))
    g.node("axles", g.mesh("axles", axles))

    # the Sun with rays and a slowly turning armillary
    sun = Geo().add("sun", sphere(0.12, 96, 48))
    n_rays = 48
    for k in range(n_rays):  # Fibonacci directions
        yk = 1 - 2 * (k + 0.5) / n_rays
        a = k * math.pi * (3 - math.sqrt(5))
        d = np.array([math.sqrt(1 - yk * yk) * math.cos(a), yk, math.sqrt(1 - yk * yk) * math.sin(a)])
        sun.add("brass", tf(lathe([(0, 0), (0.014, 0), (0.0, 0.07)], 12), align_y(d), 0.112 * d))
    g.node("sun", g.mesh("sun", sun), (0, 1.25, 0))
    arm = Geo()
    for R, rx, ry in ((0.19, 0, 0), (0.195, math.radians(23.44), 0), (0.2, math.pi / 2, 0), (0.205, math.pi / 2, math.pi / 2)):
        arm.add("brass", tf(torus(R, 0.0045, 160, 10), rot((0, 1, 0), ry) @ rot((1, 0, 0), rx)))
    for k in range(36):
        a = k * TAU / 36
        arm.add("bronze", tf(sphere(0.007, 12, 8), rot((1, 0, 0), math.radians(23.44)),
                             (0.195 * math.cos(a), 0, 0.195 * math.sin(a))))
    g.node("armillary", g.mesh("armillary", arm), (0, 1.25, 0), spin=0.05)

    # planets: arm, counterweight, rod, planet node (with tilt and own spin), moons and rings
    planets = [("mercury", 0.30, 1.10, 0.022, 0.0), ("venus", 0.42, 1.05, 0.034, 3.0), ("earth", 0.56, 1.00, 0.036, 23.4),
               ("mars", 0.70, 0.95, 0.027, 25.2), ("jupiter", 0.86, 0.90, 0.072, 3.1), ("saturn", 1.04, 0.85, 0.058, 26.7)]
    for i, (name, R, ay, size, tilt) in enumerate(planets):
        top = 1.25 - size - 0.012
        a = Geo().add("brass", washer(0.022, 0.042, ay - 0.016, ay + 0.016, 48))
        a.add("brass", rod((0.03, ay, 0), (R, ay, 0), 0.006))
        a.add("brass", rod((-0.03, ay, 0), (-0.16 - 0.02 * i, ay, 0), 0.006))
        a.add("bronze", tf(sphere(0.022 + 0.003 * i, 32, 16), None, (-0.16 - 0.02 * i, ay, 0)))
        a.add("brass", tf(sphere(0.011, 24, 12), None, (R, ay, 0)))
        a.add("brass", rod((R, ay, 0), (R, top, 0), 0.0045))
        a.add("brass", tf(lathe([(0, 0)] + ogee(0.004, 0, 0.014, 0.012) + [(0, 0.012)], 32), None, (R, top, 0)))
        node = g.node(f"orbit_{name}", g.mesh(f"arm_{name}", a), q=quat((0, 1, 0), 2.4 * i + 0.6),
                      spin=round(0.5 * (0.30 / R) ** 1.5, 4))
        tilt_n = g.node(f"{name}_tilt", t=(R, 1.25, 0), q=quat((0, 0, 1), math.radians(tilt)), parent=node)
        V, F = sphere(size, 96 if size > 0.03 else 64, 48 if size > 0.03 else 32)
        body = Geo().add("planet", (V, F), paint(name, V))
        if name == "saturn":
            prof = [(r, -0.0015) for r in np.linspace(0.075, 0.135, 25)] + [(r, 0.0015) for r in np.linspace(0.135, 0.075, 25)]
            Vr, Fr = lathe(prof, 160, closed=True)
            rr = np.linalg.norm(Vr[:, [0, 2]], axis=1)
            band = 0.5 + 0.5 * np.sin(rr * 420)
            col_r = mix((0.85, 0.76, 0.58), (0.62, 0.52, 0.38), band)
            col_r = np.where(((rr > 0.110) & (rr < 0.116))[:, None], (0.12, 0.1, 0.08), col_r)  # Cassini division
            body.add("planet", (Vr, Fr), col_r)
        g.node(name, g.mesh(name, body), parent=tilt_n, spin=1.2)
        if name == "earth":
            moon_arm = g.node("moon_orbit", parent=tilt_n, spin=1.6)
            Vm, Fm = sphere(0.011, 32, 16)
            g.node("moon", g.mesh("moon", Geo().add("planet", (Vm, Fm), paint("moon", Vm))
                                  .add("steel", rod((-0.072, 0, 0), (-0.04, 0, 0), 0.0015))), (0.072, 0, 0), parent=moon_arm)

    # hand crank on the side of the pedestal (turns around x)
    crank = Geo().add("steel", rod((-0.03, 0, 0), (0.07, 0, 0), 0.008))
    crank.add("brass", tf(box(0.012, 0.13, 0.024), None, (0.07, 0.05, 0)), flat=True)
    crank.add("brass", tf(sphere(0.02, 24, 12), None, (0.07, 0, 0)))
    handle = [(0, 0)] + hard((0.007, 0)) + ogee(0.007, 0.0, 0.016, 0.03) + ogee(0.016, 0.03, 0.009, 0.07)[1:] + [(0, 0.075)]
    crank.add("ivory", tf(lathe(handle, 32), rot((0, 0, 1), -math.pi / 2), (0.076, 0.11, 0)))
    g.node("crank", g.mesh("crank", crank), (0.89, 0.21, 0), spin=1.5, axis="x")

    g.save(out)
    return g


def main():
    out = sys.argv[1] if len(sys.argv) > 1 else os.path.join(os.path.dirname(os.path.abspath(__file__)), "orrery.glb")
    g = build(out)
    print(f"{out}: {len(g.j['nodes'])} nodes, {len(g.j['meshes'])} meshes, {g.tris:,} triangles, "
          f"{os.path.getsize(out) / 1e6:.1f} MB")


if __name__ == "__main__":
    main()