code for animating strange attractors in blender

attractors.py

import bpy
import bmesh
import random
from mathutils import Vector, noise, Matrix
from math import sin, cos, tau, pi, radians
from utils.interpolation import *
from easing_functions import *

frame_start = 1
total_frames = 120
bpy.context.scene.render.fps = 30
bpy.context.scene.frame_start = frame_start
bpy.context.scene.frame_end = total_frames
random.seed(1)

dt = 0.02
num_paths = 700
points_per_path = 100
start_pos_scale = 2.0
initial_path_steps = 50
scale = 1.0
bounds_limit = 6.0
path_offsets = [random.uniform(-0.09, 0.09) for _ in range(num_paths)]
# colors: 2D24D0,faf0ca,f4d35e,ee964b,f95738
# colors: 8c224b,c11425,ed7d51,df633d,465664
# colors: 4059ad,78e3f6,3cff7a,33d040,d8fa2d

# https://www.dynamicmath.xyz/strange-attractors/

def thomas(pos):
    b = 0.208186
    x, y, z = pos
    x1 = x + (sin(y) - b * x) * dt
    y1 = y + (sin(z) - b * y) * dt
    z1 = z + (sin(x) - b * z) * dt
    return x1, y1, z1


def dradas(pos):
    x, y, z = pos
    a = 3.0
    b = 2.7
    c = 1.7
    d = 2.0
    e = 9.0
    x1 = x + (y - a * x + b * y * z) * dt
    y1 = y + (c * y - x * z + z) * dt
    z1 = z + (d*x*y - e * z) * dt
    return x1, y1, z1


def rabinovich_fabrikant(pos):
    x, y, z = pos
    a = 0.14
    g = 0.10
    x1 = x + (y * (z - 1 + x*x) + g * x) * dt
    y1 = y + (x * (3 * z + 1 - x * x) + g * y) * dt
    z1 = z + (-2 * z * (a + x * y)) * dt
    return x1, y1, z1


def compute_path(pos):
    path = []
    for _ in range(initial_path_steps):
        pos = rabinovich_fabrikant(pos)
        # if Vector(pos).length > 6.0:
        if max(pos) > bounds_limit:
            return None
    for _ in range(points_per_path):
        pos = rabinovich_fabrikant(pos)
        # if Vector(pos).length > 6.0:
        if max(pos) > bounds_limit:
            return None
        path.append(pos)
    return path


def compute_paths():
    paths = []
    while len(paths) < num_paths:
        path = []
        start_pos = (
            random.uniform(-start_pos_scale, start_pos_scale),
            random.uniform(-start_pos_scale, start_pos_scale),
            random.uniform(-start_pos_scale, start_pos_scale),
        )
        path = compute_path(start_pos)
        if path:
            paths.append(path)
    return paths


paths = compute_paths()


def create_spline(curve):
    spline = curve.splines.new(type='NURBS')
    spline.use_cyclic_u = False
    spline.points.add(points_per_path - 1)
    spline.use_endpoint_u = True


def setup():
    col = bpy.data.collections.get('generated')
    if not col:
        col = bpy.data.collections.new('generated')
        bpy.context.scene.collection.children.link(col)
        material = bpy.data.materials.get('wire')
        parent_empty = bpy.data.objects.get('parent')
        if not parent_empty:
            parent_empty = bpy.data.objects.new('parent', None)
            bpy.data.collections['Collection'].objects.link(parent_empty)

        for i in range(num_paths):
            curve = bpy.data.curves.new('curve', type='CURVE')
            curve.dimensions = '3D'
            curve.bevel_depth = 0.033
            curve.materials.append(material)
            curve.use_fill_caps = False
            curve.taper_object = bpy.data.objects["BezierCurve"]

            # for _ in range(num_paths):
            create_spline(curve)

            # for i, spline in enumerate(curve.splines):
            path = paths[i]
            for j, point in enumerate(curve.splines[0].points):
                x, y, z = path[j]
                x *= scale
                y *= scale
                z *= scale
                point.co = (x, y, z, 1.0)

            obj = bpy.data.objects.new('curves', curve)
            col.objects.link(obj)
            obj.parent = parent_empty


def frame_update(scene, degp):
    frame = scene.frame_current
    t = (((frame - 1) % total_frames) + 1) / float(total_frames)

    objects = bpy.data.collections['generated'].all_objects
    for i, obj in enumerate(objects):
        if not obj.data: continue
        # time_offset = path_offsets[i]
        # tt = ((t + time_offset) % 1.0) * 3
        # t_end = linearstep(1.0, 2.0, tt)
        # t_start = linearstep(1.0, 2.0, tt - 1.0)

        tt = linearstep(0.1, 0.9, t + path_offsets[i])
        t_end = QuadEaseOut().ease(linearstep(0.0, 0.65, tt))
        t_start = QuadEaseIn().ease(linearstep(0.35, 1.0, tt))

        curve = obj.data
        curve.bevel_factor_end = t_end
        curve.bevel_factor_start = t_start

    bpy.data.objects.get('parent').rotation_euler.z = mix(radians(7.0), radians(-7.0), t)
    # camera = bpy.data.objects.get('Camera')
    # camera.location.y = 1.0 * sin(tau * t)


setup()
bpy.app.handlers.frame_change_pre.clear()
bpy.app.handlers.frame_change_post.clear()
bpy.app.handlers.frame_change_post.append(frame_update)
# bpy.app.handlers.frame_change_pre.append(frame_update)