superhelix.py

#!/usr/bin/env python3
"""
superhelix.py
20-second 3D voxel animation of a morphing superhelix with DNA-like structure.
Features dynamic radius oscillation, color gradients, particle trails, and smooth rotations.

The superhelix equation: x = (r + s sin n t) cos t, y = (r + s sin n t) sin t, z = s cos n t
Where r is base radius, s is oscillation amplitude, n controls frequency of bulges.

Run:
  pip install spatialstudio numpy
  python superhelix.py

Outputs:
  superhelix.splv
"""

import math
import numpy as np
from colorsys import hsv_to_rgb
from spatialstudio import splv

# -------------------------------------------------
GRID = 128            # cubic voxel grid size
FPS = 30              # frames per second
DURATION = 20         # seconds
COUNT = 600           # number of particles along helix
OUTPUT = "../outputs/superhelix.splv"
# -------------------------------------------------

TOTAL_FRAMES = FPS * DURATION
CENTER = np.array([GRID // 2] * 3, dtype=float)

# Superhelix parameters
BASE_RADIUS = GRID * 0.15      # base radius of helix
OSC_AMPLITUDE = GRID * 0.08    # amplitude of radius oscillation
HARMONIC_FREQ = 4              # frequency of radius oscillation (n parameter)
HELIX_HEIGHT = GRID * 0.6      # total height of helix
TURNS = 3                      # number of complete turns along height

# Animation parameters
ROTATION_SPEED = 1.5           # rotations per duration around Y axis
MORPH_CYCLES = 2               # how many times parameters morph during animation
TRAIL_LENGTH = 8               # length of particle trails

def smoothstep(edge0: float, edge1: float, x: float) -> float:
    """Smooth interpolation function"""
    t = max(0.0, min(1.0, (x - edge0) / (edge1 - edge0)))
    return t * t * (3 - 2 * t)

def hsv_bytes(h: float, s: float = 1.0, v: float = 1.0) -> tuple:
    """Convert HSV to RGB bytes"""
    r, g, b = hsv_to_rgb(h % 1.0, s, v)
    return int(r * 255), int(g * 255), int(b * 255)

def superhelix_pos(t_param: float, time_factor: float = 0.0) -> np.ndarray:
    """
    Calculate position on superhelix with time-varying parameters
    t_param: parameter along helix (0 to 2π * TURNS)
    time_factor: global animation time (0 to 1) for morphing effects
    """
    # Morph the harmonic frequency over time for visual interest
    n = HARMONIC_FREQ + 2 * math.sin(time_factor * MORPH_CYCLES * 2 * math.pi)
    
    # Morph oscillation amplitude
    s = OSC_AMPLITUDE * (0.5 + 0.5 * math.cos(time_factor * MORPH_CYCLES * 2 * math.pi))
    
    # Base radius with gentle variation
    r = BASE_RADIUS * (1 + 0.2 * math.sin(time_factor * math.pi))
    
    # Superhelix equations
    radius_at_t = r + s * math.sin(n * t_param)
    x = radius_at_t * math.cos(t_param)
    y = radius_at_t * math.sin(t_param)
    z = s * math.cos(n * t_param) + (t_param / (2 * math.pi * TURNS)) * HELIX_HEIGHT - HELIX_HEIGHT/2
    
    return np.array([x, y, z])

def rotate_y(vec: np.ndarray, angle: float) -> np.ndarray:
    """Rotate vector around Y axis"""
    cos_a, sin_a = math.cos(angle), math.sin(angle)
    return np.array([
        vec[0] * cos_a + vec[2] * sin_a,
        vec[1],
        -vec[0] * sin_a + vec[2] * cos_a
    ])

def rotate_x(vec: np.ndarray, angle: float) -> np.ndarray:
    """Rotate vector around X axis"""
    cos_a, sin_a = math.cos(angle), math.sin(angle)
    return np.array([
        vec[0],
        vec[1] * cos_a - vec[2] * sin_a,
        vec[1] * sin_a + vec[2] * cos_a
    ])

# Pre-compute parameter values along helix
t_vals = np.linspace(0, 2 * math.pi * TURNS, COUNT, endpoint=False)

# Initialize particle trail history
particle_trails = [[] for _ in range(COUNT)]

enc = splv.Encoder(GRID, GRID, GRID, framerate=FPS, outputPath=OUTPUT)
print(f"Encoding {TOTAL_FRAMES} frames for superhelix animation...")

for frame_idx in range(TOTAL_FRAMES):
    global_time = frame_idx / TOTAL_FRAMES  # 0 to 1
    
    # Rotation angles
    rot_y = global_time * ROTATION_SPEED * 2 * math.pi
    rot_x = math.sin(global_time * math.pi) * 0.3  # gentle nodding motion
    
    frame = splv.Frame(GRID, GRID, GRID)
    
    for i in range(COUNT):
        # Calculate position on superhelix
        helix_pos = superhelix_pos(t_vals[i], global_time)
        
        # Apply rotations
        rotated_pos = rotate_y(helix_pos, rot_y)
        rotated_pos = rotate_x(rotated_pos, rot_x)
        
        # Translate to center
        world_pos = CENTER + rotated_pos
        
        # Add to trail history
        particle_trails[i].append(world_pos.copy())
        if len(particle_trails[i]) > TRAIL_LENGTH:
            particle_trails[i].pop(0)
        
        # Render particle trail
        for trail_idx, trail_pos in enumerate(particle_trails[i]):
            x, y, z = trail_pos.astype(int)
            
            if 0 <= x < GRID and 0 <= y < GRID and 0 <= z < GRID:
                # Color based on position along helix and trail age
                hue_base = (i / COUNT) * 0.8  # spread across most of hue spectrum
                hue_offset = global_time * 0.2  # slowly shift colors over time
                hue = (hue_base + hue_offset) % 1.0
                
                # Trail fading effect
                trail_alpha = (trail_idx + 1) / len(particle_trails[i])
                saturation = 0.8 + 0.2 * trail_alpha
                brightness = 0.4 + 0.6 * trail_alpha
                
                # Add some sparkle to the newest trail points
                if trail_idx == len(particle_trails[i]) - 1:
                    brightness = min(1.0, brightness + 0.3 * math.sin(global_time * 20 + i * 0.1))
                
                color = hsv_bytes(hue, saturation, brightness)
                frame.set_voxel(x, y, z, color)
                
                # Add extra brightness for main particle (trail head)
                if trail_idx == len(particle_trails[i]) - 1:
                    # Add neighboring voxels for larger particle effect
                    for dx in [-1, 0, 1]:
                        for dy in [-1, 0, 1]:
                            for dz in [-1, 0, 1]:
                                nx, ny, nz = x + dx, y + dy, z + dz
                                if (0 <= nx < GRID and 0 <= ny < GRID and 0 <= nz < GRID and 
                                    abs(dx) + abs(dy) + abs(dz) <= 1):
                                    bright_color = hsv_bytes(hue, saturation, min(1.0, brightness * 0.7))
                                    frame.set_voxel(nx, ny, nz, bright_color)
    
    enc.encode(frame)
    
    # Progress indicator
    if frame_idx % FPS == 0:
        seconds_done = frame_idx // FPS + 1
        print(f"  second {seconds_done} / {DURATION}")

enc.finish()
print("Done. Saved", OUTPUT)