orrery.py
Автор: User 254f7fd2
23 KB · text/x-pythonОткрыть в JupyterСкачать
"""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()