1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
|
# -*- 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 +
|
