diff options
Diffstat (limited to 'blender/caos.py')
| -rw-r--r-- | blender/caos.py | 219 |
1 files changed, 219 insertions, 0 deletions
diff --git a/blender/caos.py b/blender/caos.py new file mode 100644 index 0000000..99ef9db --- /dev/null +++ b/blender/caos.py @@ -0,0 +1,219 @@ +# -*- coding: utf-8 -*- +"""CAOS - tres pendulos dobles con 1 mm de diferencia inicial. + +La fisica es un RK4 sobre las ecuaciones exactas del pendulo doble, en unidades +SI (L1 = L2 = 1 m), integrada antes de renderizar y muestreada por frame. El +milimetro del guion es literal: 0,001 rad sobre una varilla de 1 m. +""" +import math, os, sys +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) +import bpy +from base import * + +NOMBRE = "caos" +G, L1, L2, M1, M2 = 9.81, 1.0, 1.0, 1.0, 1.0 +TH0 = 2.4 # 138 grados: regimen caotico, separacion visible a los ~2,5 s +DELTA = 0.001 # radianes = 1 mm en la punta de la primera varilla +ESC = 0.80 # escala de dibujo (la fisica sigue en metros) +ALTO = 6.0 # metros visibles de alto en el encuadre +PZ_A, PZ_B = 0.47, 1.11 # altura del pivote antes / despues del beat 8 +K_B = 0.66 # el rig se achica para dejar lugar al grafico +# banda util medida sobre el overlay: el chip llega a z=+2,27 y la caja de +# subtitulos arranca en z=-1,09 (tres lineas) / -1,25 (dos lineas) +GX0, GZ0, GW, GH = -1.45, -0.92, 2.90, 0.80 +COLS = ["ambar", "rosa", "celeste"] +ESTELA = 34 # frames de cola +BEAT_SUELTA = 4 # se sueltan al empezar este beat + + +# --- integrador --------------------------------------------------------------- +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 articulaciones(s): + """(codo, punta) en metros, relativo al pivote.""" + 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 integrar(n_frames, sub=20): + """Estados de los tres pendulos, muestreados por frame.""" + est = [(TH0 + k * DELTA, TH0, 0.0, 0.0) for k in range(3)] + h = 1.0 / (FPS * sub) + tray = [[articulaciones(s) for s in est]] + for _ in range(n_frames): + for _ in range(sub): + est = [rk4(s, h) for s in est] + tray.append([articulaciones(s) for s in est]) + return tray + + +# --- escena ------------------------------------------------------------------- +def construir(T): + sc = escena() + lente = 70.0 + dist = lente / 36.0 * ALTO # distancia para ver ALTO metros de alto + # la camara se centra en z=0: asi las coordenadas coinciden con el cuadro + # y se puede ubicar todo respecto del chip y de los subtitulos + cam = camara((0.0, -dist, 0.0), (0.0, 0.0, 0.0), lente=lente) + cam.data.sensor_fit = 'VERTICAL' + cam.data.sensor_height = 36.0 + + luz("key", 'AREA', (-3.6, -5.2, 4.6), 900, "blanco", tam=5.0, mira=(0, 0, PZ_A)) + luz("fill", 'AREA', (4.2, -4.6, -0.6), 300, "celeste", tam=5.0, mira=(0, 0, PZ_A)) + luz("rim", 'AREA', (0.8, 4.6, 2.6), 520, "blanco", tam=4.0, mira=(0, 0, PZ_A)) + + # cubo del pivote: un disco metalico mirando a camara, sin soporte que tape + m_hub = material("hub", "riel", rug=0.35, metal=0.9) + hub = cilindro("hub", 0.105, 0.10, m_hub) + hub.rotation_euler = (math.pi / 2, 0, 0) + aro = cilindro("aro", 0.145, 0.05, material("aro", "gris", rug=0.25, metal=1.0, emis=0.5)) + aro.rotation_euler = (math.pi / 2, 0, 0) + + pend = [] + for k, c in enumerate(COLS): + m = material(f"p{k}", c, rug=0.28, metal=0.55, emis=0.25) + m_bola = material(f"b{k}", c, rug=0.18, metal=0.2, emis=0.9) + pend.append({ + "y": -0.035 * (k - 1), "col": c, + "v1": cilindro(f"v1_{k}", 0.034, 1.0, m), + "v2": cilindro(f"v2_{k}", 0.029, 1.0, m), + "codo": esfera(f"codo_{k}", 0.055, m_bola), + "punta": esfera(f"punta_{k}", 0.090, m_bola), + "estela": curva_poly(f"estela_{k}", [[(0, 0, 0)] * ESTELA], grosor=0.032, + radios=[[0.0] * ESTELA], + mat=material(f"e{k}", c, rug=0.5, emis=2.4)), + }) + + leyenda = [] + for k, (c, txt) in enumerate(zip(COLS, ("arranca +0 mm", "arranca +1 mm", "arranca +2 mm"))): + t = texto(txt, tam=0.23, color=c, align='LEFT') + t.location = (-1.25, -0.6, -0.28 - 0.38 * k) + p = esfera(f"pt_{k}", 0.075, material(f"pm{k}", c, emis=2.2)) + p.location = (-1.42, -0.6, -0.28 - 0.38 * k) + leyenda.append((t, p)) + + m_ejes = material("ejes", "gris", rug=0.6, emis=0.9) + graf = { + "x": cilindro("gx", 0.014, GW, m_ejes), + "y": cilindro("gy", 0.014, GH, m_ejes), + "curva": curva_poly("gcurva", [[(0, 0, 0)] * 2], grosor=0.030, + radios=[[1.0] * 2], mat=material("gc", "rosa", emis=2.8)), + "arr": texto("2 metros", tam=0.20, color="gris", align='LEFT'), + "aba": texto("1 mm", tam=0.20, color="gris", align='RIGHT'), + } + graf["x"].rotation_euler = (0, math.pi / 2, 0) + graf["x"].location = (GX0 + GW / 2, 0.25, GZ0) + graf["y"].location = (GX0, 0.25, GZ0 + GH / 2) + # las etiquetas van donde la curva no pasa: arriba a la izquierda y abajo a + # la derecha (la curva sale de abajo-izquierda y termina arriba-derecha) + graf["arr"].location = (GX0 + 0.12, 0.25, GZ0 + GH + 0.02) + graf["aba"].location = (GX0 + GW - 0.10, 0.25, GZ0 + 0.13) + return dict(cam=cam, pend=pend, leyenda=leyenda, graf=graf, hub=hub, aro=aro) + + +def main(): + T = Tiempo(NOMBRE) + f_suelta = T.rango(BEAT_SUELTA)[0] + n_sim = T.n_frames - f_suelta + 2 + print(f"[{NOMBRE}] integrando {n_sim} frames de fisica...", flush=True) + tray = integrar(n_sim) + obj = construir(T) + pend, graf = obj["pend"], obj["graf"] + + # separacion entre la punta 0 y la punta 2, en metros, por frame + sep = [math.dist(p[0][1], p[2][1]) for p in tray] + LMIN, LMAX = math.log10(0.001), math.log10(2.0) + + # separacion maxima alcanzada: monotona, se lee de un vistazo + 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 actualizar(f): + idx = max(0, min(f - f_suelta, len(tray) - 1)) + # el rig se achica y sube en el beat 8 para dejarle lugar al grafico + m = suave(T.p(f, 8)) + pz = PZ_A + (PZ_B - PZ_A) * m + k_esc = ESC * (1.0 + (K_B - 1.0) * m) + obj["hub"].location = (0, 0.06, pz) + obj["aro"].location = (0, 0.02, pz) + obj["hub"].scale = obj["aro"].scale = (1 - 0.3 * m,) * 3 + + abanico = (1.0 - suave(T.p(f, 0))) * 0.34 if f < f_suelta else 0.0 + for k, g in enumerate(pend): + (cx, cz), (px, pz2) = tray[idx][k] + if abanico: + ang = TH0 + (k - 1) * abanico + 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) + codo = (cx * k_esc, y, pz + cz * k_esc) + punta = (px * k_esc, y, pz + pz2 * k_esc) + orientar(g["v1"], o, codo) + orientar(g["v2"], codo, punta) + g["v1"].scale = (1, 1, L1 * k_esc) + g["v2"].scale = (1, 1, L2 * k_esc) + g["codo"].location = codo + g["punta"].location = punta + g["codo"].scale = g["punta"].scale = (1 - 0.3 * m,) * 3 + # la estela sale de la trayectoria, no se acumula: vale con frames salteados + pts, rad = [], [] + for j in range(ESTELA): + i2 = idx - (ESTELA - 1 - j) + q = tray[max(0, i2)][k][1] + pts.append((q[0] * k_esc, y, pz + q[1] * k_esc)) + vivo = 1.0 if (i2 > 0 and f >= f_suelta) else 0.0 + rad.append(vivo * (j / (ESTELA - 1.0)) ** 2.2) + rehacer_curva(g["estela"], [pts], [rad]) + + vis = suave(T.p(f, 1)) * (1.0 - suave(T.p(f, 3))) + for t, p in obj["leyenda"]: + t.scale = p.scale = (vis, vis, vis) + + gv = suave(T.p(f, 8)) + for k in ("x", "y", "arr", "aba"): + graf[k].scale = (gv, gv, gv) + n_g = min(len(sep), int(8.0 * FPS)) # ventana: los primeros 8 s + avance = 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 <= avance else 0.0) + if len(graf["curva"].data.splines[0].points) != len(pts): + graf["curva"].data.splines.clear() + sp = graf["curva"].data.splines.new('POLY') + sp.points.add(len(pts) - 1) + rehacer_curva(graf["curva"], [pts], [[r * gv for r in rad]]) + + if os.environ.get("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_secuencia(NOMBRE, T, actualizar) + + +main() |