"""Dirac belt (the belt trick, or plate trick): a cube tied to fixed anchor pads by elastic ribbons, posed at any moment of a full 720 degree untangle, with the cube's ribbons never passing through each other or the cube. Every point p in space is turned rigidly about the centre by a rotation that depends only on its distance |p| from the centre: points inside the ball of radius R0 around the cube turn with the cube, points outside R1 stay still, and each sphere in between turns as one piece. That map is smooth, one-to-one and volume-preserving at every moment, so ribbons that start apart stay apart and keep their volume. The motion follows Andrew Hanson's 360 and 720 degree quaternion paths (Visualizing Quaternions, 2006), then a six-ribbon finale. Ribbons are B-spline solids fitted through the exact mapped points (top and bottom surfaces, ruled edges and ends). The cube sits at the origin; ribbon +z points up. Its +x face is marked red so its turn is visible. """ import math import cadquery as cq from steve import Model time = 24.0 # s Moment in the motion ribbons = "as in the video" # Ribbons cube_size = 40 # mm Cube size ribbon_width = 18 # mm Ribbon width ribbon_thickness = 0.8 # mm Ribbon thickness ribbon_length = 108 # mm Ribbon length # ---- geometry, in the cube's half-size `a` (the video's units: a = 0.5) ---- a = cube_size / 2 W = ribbon_width TH = ribbon_thickness L = a + ribbon_length # anchor face distance from the centre R0 = 1.76 * a # rigid ball: turns with the cube (cube corner is 1.732 a) R1 = L - 0.3 * a # outside this nothing moves DELTA = 0.15 # ramp width of the log-radius profile PAD_HALF = 0.9 * a PAD_T = 0.2 * a # Proportions the fitted ribbon surfaces are checked for (the video: width 0.45, thickness 0.02, length 2.7 x cube) if not W <= 0.5 * cube_size: raise ValueError("ribbon_width must be at most 0.5 x cube_size.") if not TH <= 0.02 * cube_size: raise ValueError("ribbon_thickness must be at most 0.02 x cube_size.") if not TH <= W / 4: raise ValueError("ribbon_thickness must be at most a quarter of ribbon_width.") if not ribbon_length >= 2.5 * cube_size: raise ValueError("ribbon_length must be at least 2.5 x cube_size.") # Timeline in seconds, as in the 37.5 s video: rest, 360 degree turn (A1), Hanson 360 loop, second turn (A2), # Hanson 720 untangle, four side ribbons grow, six-ribbon finale loop. PHASES = [("rest", 0.0, 3.0), ("A1", 3.0, 9.0), ("H360", 9.0, 15.0), ("A2", 15.0, 21.0), ("H720", 21.0, 28.0), ("grow", 28.0, 29.5), ("coda", 29.5, 37.5)] # (name, axis e, width direction u); the ribbon's front (gold) face looks along e x u RIBBONS = [("pz", (0, 0, 1), (0, -1, 0)), ("nz", (0, 0, -1), (0, 1, 0)), ("px", (1, 0, 0), (0, 1, 0)), ("nx", (-1, 0, 0), (0, -1, 0)), ("py", (0, 1, 0), (-1, 0, 0)), ("ny", (0, -1, 0), (1, 0, 0))] SHOWN = {"as in the video": None, "1": 1, "2": 2, "6": 6}[ribbons] def ease(x): x = min(max(x, 0.0), 1.0) return x * x * x * (10 - 15 * x + 6 * x * x) def _E(x): return 2.5 * x ** 4 - 3 * x ** 5 + x ** 6 # integral of smootherstep def h(r): """distance from the centre -> path parameter s: 1 inside R0, 0 outside R1, smooth (C3) and monotone""" rho = min(max(math.log(max(r, 1e-12) / R0) / math.log(R1 / R0), 0.0), 1.0) d = DELTA if rho < d: w = d * _E(rho / d) elif rho > 1 - d: w = d / 2 + (1 - 2 * d) + d * (0.5 - _E((1 - rho) / d)) else: w = d / 2 + (rho - d) return 1.0 - w / (1.0 - d) T = min(max(float(time), 0.0), 37.5) PHASE, U = "coda", 1.0 for _name, _t0, _t1 in PHASES: if T <= _t1 + 1e-12: PHASE, U = _name, ease((T - _t0) / (_t1 - _t0)) break def qmul(p, q): aw, ax, ay, az = p bw, bx, by, bz = q return (aw * bw - ax * bx - ay * by - az * bz, aw * bx + ax * bw + ay * bz - az * by, aw * by - ax * bz + ay * bw + az * bx, aw * bz + ax * by - ay * bx + az * bw) def Q(s): """unit quaternion (w, x, y, z) for path parameter s (0 = anchors, 1 = cube) at this moment""" if PHASE in ("rest", "grow"): return (1.0, 0.0, 0.0, 0.0) if PHASE in ("A1", "A2"): tau = 0.5 * U if PHASE == "A1" else 0.5 + 0.5 * U ang = 2 * math.pi * s * tau return (math.cos(ang), 0.0, 0.0, math.sin(ang)) if PHASE == "H360": al = 2 * math.pi * U c, sn = math.cos(math.pi * s), math.sin(math.pi * s) return (c, 0.0, sn * math.sin(al), sn * math.cos(al)) if PHASE == "H720": b = 0.5 * math.pi * U c2 = math.cos(2 * math.pi * s) return (math.sin(b) ** 2 + math.cos(b) ** 2 * c2, 0.0, math.cos(b) * math.sin(b) * (1 - c2), math.cos(b) * math.sin(2 * math.pi * s)) hs = 0.5 * math.pi * s # coda: loop of the six-ribbon finale al = 2 * math.pi * U q1 = (math.cos(hs), math.sin(hs) * math.cos(al), math.sin(hs) * math.sin(al), 0.0) q0c = (math.cos(hs), -math.sin(hs), 0.0, 0.0) return qmul(q1, q0c) def rotate(q, p): w, vx, vy, vz = q px, py, pz = p tx, ty, tz = 2 * (vy * pz - vz * py), 2 * (vz * px - vx * pz), 2 * (vx * py - vy * px) return (px + w * tx + (vy * tz - vz * ty), py + w * ty + (vz * tx - vx * tz), pz + w * tz + (vx * ty - vy * tx)) def phi(p): return rotate(Q(h(math.sqrt(p[0] ** 2 + p[1] ** 2 + p[2] ** 2))), p) def grown(index): """(pad scale, ribbon length fraction) of ribbon `index` at this moment""" if SHOWN is not None: return (1.0, 1.0) if index < SHOWN else (0.0, 0.0) if index < 2: return 1.0, 1.0 if T < 28.0: return 0.0, 0.0 g = (T - 28.0) / 1.5 return ease(g / 0.5), ease((g - 0.3) / 0.7) def box_at(centre, e, u, size_e, size_u, size_n): """box of the given sizes along e, u and n = e x u, centred on `centre`""" n = (e[1] * u[2] - e[2] * u[1], e[2] * u[0] - e[0] * u[2], e[0] * u[1] - e[1] * u[0]) o = [centre[i] - (e[i] * size_e + u[i] * size_u + n[i] * size_n) / 2 for i in range(3)] return cq.Solid.makeBox(size_e, size_u, size_n).moved(cq.Location(cq.Plane(o, e, n))) def ribbon_solid(e, u, frac, nl=240, nw=12): n = (e[1] * u[2] - e[2] * u[1], e[2] * u[0] - e[0] * u[2], e[0] * u[1] - e[1] * u[0]) tip = a + frac * (L - a) def pt(l, w, k): return tuple(l * e[i] + w * u[i] + k * n[i] for i in range(3)) if frac < 1.0: # growing: only in the grow phase, where nothing has turned yet, so it is a plain slab return box_at(pt((a + tip) / 2, 0, 0), e, u, tip - a, W, TH) # length samples are denser in the moving shell R0..R1 (where ribbons stretch) than in the still parts ls = [] for i in range(nl + 1): x = i / nl ls.append(a + (L - a) * (x - 0.10 * math.sin(2 * math.pi * x) / (2 * math.pi))) ws = [-W / 2 + W * j / nw for j in range(nw + 1)] layers = [] for k in (TH / 2, -TH / 2): layers.append(cq.Face.makeSplineApprox( [[cq.Vector(*phi(pt(l, w, k))) for w in ws] for l in ls], tol=5e-5 * a, minDeg=3, maxDeg=5)) top, bot = layers faces = [top, bot] for et in top.Edges(): # sides and ends: ruled between matching top and bottom edges mt = et.positionAt(0.5) eb = min(bot.Edges(), key=lambda x: (x.positionAt(0.5) - mt).Length) if (et.startPoint() - eb.startPoint()).Length > (et.startPoint() - eb.endPoint()).Length: eb = cq.Edge(eb.wrapped.Reversed()) faces.append(cq.Face.makeRuledSurface(et, eb)) return cq.Solid.makeSolid(cq.Shell.makeShell(faces)).fix() # cube: turns rigidly with Q(1); edges rounded; +x face marked qc = Q(1.0) mark = box_at((a - 0.02 * a, 0, 0), (1, 0, 0), (0, 1, 0), 0.04 * a, 1.6 * a, 1.6 * a) # inlaid, flush cube = cq.Workplane("XY").box(2 * a, 2 * a, 2 * a).edges().fillet(0.07 * a).val().cut(mark) ang = 2 * math.degrees(math.acos(min(1.0, abs(qc[0])))) if ang > 1e-9: axis = (qc[1], qc[2], qc[3]) if qc[0] >= 0 else (-qc[1], -qc[2], -qc[3]) cube = cube.rotate(cq.Vector(0, 0, 0), cq.Vector(*axis), ang) mark = mark.rotate(cq.Vector(0, 0, 0), cq.Vector(*axis), ang) model = Model() model.add("cube", cube, op="new", name="Cube") model.add("mark", mark, op="new", name="Marked face") pads, ribs = [], [] for index, (key, e, u) in enumerate(RIBBONS): gp, gr = grown(index) if gp > 1e-3: c = tuple((L + PAD_T / 2) * x for x in e) model.add(f"pad_{key}", box_at(c, e, u, gp * PAD_T, gp * 2 * PAD_HALF, gp * 2 * PAD_HALF), op="new", name=f"Anchor pad {key}") pads.append(f"pad_{key}") if gr > 1e-3: model.add(f"ribbon_{key}", ribbon_solid(e, u, gr), op="new", name=f"Ribbon {key}") ribs.append(f"ribbon_{key}") model.material("aluminum-6061-t6", bodies=["cube"], finish="brushed", color="#cfd1d4") model.material("pa6", bodies=["mark"], color="#d91f0a") # red inlay on the +x face if pads: model.material("aluminum-6061-t6", bodies=pads, finish="anodized", color="#2a2c31") if ribs: model.material("tpu-60a", bodies=ribs, finish="satin", color="#cc7a1a") # elastic: ribbons stretch result = model.result()