# -*- coding: utf-8 -*- """CAOS - three double pendulums 1 mm apart at the start. The physics is RK4 on the exact double-pendulum equations, in SI units (L1 = L2 = 1 m), integrated before rendering and sampled per frame. The script's millimeter is literal: 0.001 rad on a 1 m rod. """ import math, os, sys sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import bpy from base import * NAME_KEY = "caos" G, L1, L2, M1, M2 = 9.81, 1.0, 1.0, 1.0, 1.0 TH0 = 2.4 # 138 degrees: chaotic regime, visible divergence at ~2.5 s DELTA = 0.001 # radians = 1 mm at the tip of the first rod SCALE = 0.80 # drawing scale (physics stays in meters) FRAME_H = 6.0 # meters of height visible in the frame PZ_A, PZ_B = 0.47, 1.11 # pivot height before / after beat 8 K_B = 0.66 # the rig shrinks to make room for the graph # usable band measured on the overlay: the chip reaches z=+2.27 and the # subtitle box starts at z=-1.09 (three lines) / -1.25 (two lines) GX0, GZ0, GW, GH = -1.45, -0.92, 2.90, 0.80 COLS = ["ambar", "rosa", "celeste"] ESTELA = 34 # tail frames BEAT_RELEASE = 4 # released when this beat starts # --- integrator --------------------------------------------------------------- def deriv(s): t1, t2, w1, w2 = s d = t1 - t2 den = 2 * M1 + M2 - M2 * math.cos(2 * d) a1 = (-G * (2 * M1 + M2) * math.sin(t1) - M2 * G * math.sin(t1 - 2 * t2) - 2 * math.sin(d) * M2 * (w2 * w2 * L2 + w1 * w1 * L1 * math.cos(d))) / (L1 * den) a2 = (2 * math.sin(d) * (w1 * w1 * L1 * (M1 + M2) + G * (M1 + M2) * math.cos(t1) + w2 * w2 * L2 * M2 * math.cos(d))) / (L2 * den) return (w1, w2, a1, a2) def rk4(s, h): k1 = deriv(s) k2 = deriv(tuple(s[i] + h / 2 * k1[i] for i in range(4))) k3 = deriv(tuple(s[i] + h / 2 * k2[i] for i in range(4))) k4 = deriv(tuple(s[i] + h * k3[i] for i in range(4))) return tuple(s[i] + h / 6 * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) for i in range(4)) def joints(s): """(elbow, tip) in meters, relative to the pivot.""" t1, t2, _, _ = s cx, cz = L1 * math.sin(t1), -L1 * math.cos(t1) return (cx, cz), (cx + L2 * math.sin(t2), cz - L2 * math.cos(t2)) def integrate(n_frames, sub=20): """States of the three pendulums, sampled per frame.""" estimate = [(TH0 + k * DELTA, TH0, 0.0, 0.0) for k in range(3)] h = 1.0 / (FPS * sub) tray = [[joints(s) for s in estimate]] for _ in range(n_frames): for _ in range(sub): estimate = [rk4(s, h) for s in estimate] tray.append([joints(s) for s in estimate]) return tray # --- scene ------------------------------------------------------------------- def build_scene(T): sc = scene_setup() lens = 70.0 dist = lens / 36.0 * FRAME_H # distance to see FRAME_H meters of height # the camera is centred on z=0: that way coordinates match the frame and # everything can be placed relative to the chip and the subtitles cam = camera_obj((0.0, -dist, 0.0), (0.0, 0.0, 0.0), lens=lens) cam.data.sensor_fit = 'VERTICAL' cam.data.sensor_height = 36.0 light_obj("key", 'AREA', (-3.6, -5.2, 4.6), 900, "blanco", size_u=5.0, sight=(0, 0, PZ_A)) light_obj("fill", 'AREA', (4.2, -4.6, -0.6), 300, "celeste", size_u=5.0, sight=(0, 0, PZ_A)) light_obj("rim", 'AREA', (0.8, 4.6, 2.6), 520, "blanco", size_u=4.0, sight=(0, 0, PZ_A)) # pivot hub: a metal disc facing the camera, with no support blocking the view m_hub = material("hub", "riel", rough=0.35, metal=0.9) hub = cylinder("hub", 0.105, 0.10, m_hub) hub.rotation_euler = (math.pi / 2, 0, 0) hoop = cylinder("aro", 0.145, 0.05, material("aro", "gris", rough=0.25, metal=1.0, emit=0.5)) hoop.rotation_euler = (math.pi / 2, 0, 0) pend = [] for k, c in enumerate(COLS): m = material(f"p{k}", c, rough=0.28, metal=0.55, emit=0.25) m_bola = material(f"b{k}", c, rough=0.18, metal=0.2, emit=0.9) pend.append({ "y": -0.035 * (k - 1), "col": c, "v1": cylinder(f"v1_{k}", 0.034, 1.0, m), "v2": cylinder(f"v2_{k}", 0.029, 1.0, m), "codo": sphere(f"codo_{k}", 0.055, m_bola), "punta": sphere(f"punta_{k}", 0.090, m_bola), "estela": curve_poly(f"estela_{k}", [[(0, 0, 0)] * ESTELA], thickness_px=0.032, radios=[[0.0] * ESTELA], mat=material(f"e{k}", c, rough=0.5, emit=2.4)), }) legend = [] for k, (c, txt) in enumerate(zip(COLS, ("arranca +0 mm", "arranca +1 mm", "arranca +2 mm"))): t = txt_m(txt, size_u=0.23, color=c, align='LEFT') t.location = (-1.25, -0.6, -0.28 - 0.38 * k) p = sphere(f"pt_{k}", 0.075, material(f"pm{k}", c, emit=2.2)) p.location = (-1.42, -0.6, -0.28 - 0.38 * k) legend.append((t, p)) m_axes = material("ejes", "gris", rough=0.6, emit=0.9) graph = { "x": cylinder("gx", 0.014, GW, m_axes), "y": cylinder("gy", 0.014, GH, m_axes), "curva": curve_poly("gcurva", [[(0, 0, 0)] * 2], thickness_px=0.030, radios=[[1.0] * 2], mat=material("gc", "rosa", emit=2.8)), "arr": txt_m("2 metros", size_u=0.20, color="gris", align='LEFT'), "aba": txt_m("1 mm", size_u=0.20, color="gris", align='RIGHT'), } graph["x"].rotation_euler = (0, math.pi / 2, 0) graph["x"].location = (GX0 + GW / 2, 0.25, GZ0) graph["y"].location = (GX0, 0.25, GZ0 + GH / 2) # the labels go where the curve does not pass: top left and bottom right # (the curve starts bottom-left and ends top-right) graph["arr"].location = (GX0 + 0.12, 0.25, GZ0 + GH + 0.02) graph["aba"].location = (GX0 + GW - 0.10, 0.25, GZ0 + 0.13) return dict(cam=cam, pend=pend, legend=legend, graph=graph, hub=hub, hoop=hoop) def main(): T = Timeline(NAME_KEY) f_release = T.span(BEAT_RELEASE)[0] n_sim = T.n_frames - f_release + 2 print(f"[{NAME_KEY}] integrating {n_sim} physics frames...", flush=True) tray = integrate(n_sim) obj = build_scene(T) pend, graph = obj["pend"], obj["graph"] # separation between tip 0 and tip 2, in meters, per frame sep = [math.dist(p[0][1], p[2][1]) for p in tray] LMIN, LMAX = math.log10(0.001), math.log10(2.0) # maximum separation reached: monotonic, it reads at a glance sep, mx = [], 0.0 for p in tray: mx = max(mx, math.dist(p[0][1], p[2][1])) sep.append(mx) LMIN, LMAX = math.log10(0.001), math.log10(2.0) def refresh(f): idx = max(0, min(f - f_release, len(tray) - 1)) # the rig shrinks and rises in beat 8 to make room for the graph m = suave(T.p(f, 8)) pz = PZ_A + (PZ_B - PZ_A) * m k_scale = SCALE * (1.0 + (K_B - 1.0) * m) obj["hub"].location = (0, 0.06, pz) obj["hoop"].location = (0, 0.02, pz) obj["hub"].scale = obj["hoop"].scale = (1 - 0.3 * m,) * 3 fan = (1.0 - suave(T.p(f, 0))) * 0.34 if f < f_release else 0.0 for k, g in enumerate(pend): (cx, cz), (px, pz2) = tray[idx][k] if fan: ang = TH0 + (k - 1) * fan cx, cz = L1 * math.sin(ang), -L1 * math.cos(ang) px, pz2 = cx + L2 * math.sin(ang), cz - L2 * math.cos(ang) y = g["y"] o = (0.0, y, pz) elbow = (cx * k_scale, y, pz + cz * k_scale) tip_pt = (px * k_scale, y, pz + pz2 * k_scale) orient(g["v1"], o, elbow) orient(g["v2"], elbow, tip_pt) g["v1"].scale = (1, 1, L1 * k_scale) g["v2"].scale = (1, 1, L2 * k_scale) g["codo"].location = elbow g["punta"].location = tip_pt g["codo"].scale = g["punta"].scale = (1 - 0.3 * m,) * 3 # the trail comes from the trajectory, it does not accumulate: works with skipped frames pts, rad = [], [] for j in range(ESTELA): i2 = idx - (ESTELA - 1 - j) q = tray[max(0, i2)][k][1] pts.append((q[0] * k_scale, y, pz + q[1] * k_scale)) alive = 1.0 if (i2 > 0 and f >= f_release) else 0.0 rad.append(alive * (j / (ESTELA - 1.0)) ** 2.2) rebuild_curve(g["estela"], [pts], [rad]) visible = suave(T.p(f, 1)) * (1.0 - suave(T.p(f, 3))) for t, p in obj["legend"]: t.scale = p.scale = (visible, visible, visible) gv = suave(T.p(f, 8)) for k in ("x", "y", "arr", "aba"): graph[k].scale = (gv, gv, gv) n_g = min(len(sep), int(8.0 * FPS)) # window: the first 8 s advance = suave((T.p(f, 8) - 0.12) / 0.6) * n_g pts, rad = [], [] for i in range(n_g): x = GX0 + GW * i / (n_g - 1.0) v = (math.log10(max(sep[i], 1e-3)) - LMIN) / (LMAX - LMIN) pts.append((x, 0.25, GZ0 + GH * max(0.0, min(1.0, v)))) rad.append(1.0 if i <= advance else 0.0) if len(graph["curva"].data.splines[0].points) != len(pts): graph["curva"].data.splines.clear() sp = graph["curva"].data.splines.new('POLY') sp.points.add(len(pts) - 1) rebuild_curve(graph["curva"], [pts], [[r * gv for r in rad]]) if env("MODE", "MODO") == "sim": for s in (0, 1, 2, 2.5, 3, 4, 5, 6, 8): i = min(int(s * FPS), len(sep) - 1) print(f" t={s:4.1f}s separacion = {sep[i]*100:8.2f} cm") return render_sequence(NAME_KEY, T, refresh) main()