DanielHabib/hopf_fibration.py Secret

## hopf_fibration.py
#!/usr/bin/env python3
"""
hopf_fibration.py
18-second 3D voxel animation of the Hopf fibration - one of the most beautiful structures in mathematics.
The Hopf fibration maps the 3-sphere to the 2-sphere, creating interlocking circles that never touch.
Each point on the 2-sphere corresponds to a circle on the 3-sphere, creating a mesmerizing fiber bundle.

The Hopf fibration demonstrates the deep connection between topology, geometry, and physics,
and is fundamental to understanding fiber bundles, gauge theory, and quantum mechanics.

Run:
  pip install spatialstudio numpy
  python hopf_fibration.py

Outputs:
  hopf_fibration.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 = 18  # seconds
COUNT = 1000  # base particle count for flow
OUTPUT = "../outputs/hopf_fibration.splv"
# -------------------------------------------------

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

# Hopf fibration parameters
SCALE = GRID * 0.28  # scale factor for the fibers
FIBER_SAMPLES = 24  # number of sample fibers to display
POINTS_PER_FIBER = 60  # points along each fiber circle
THICKNESS = 2  # thickness for solid appearance

# Flow parameters
FLOW_SPEED = 1.8  # how fast particles flow around the fibers
FIBER_ROTATION_SPEED = 0.5  # how fast the fiber bundle rotates

# Pre-compute flow parameters for flowing particles
np.random.seed(2024)
particle_fiber_indices = np.random.randint(0, FIBER_SAMPLES, COUNT)
particle_t_offsets = np.random.rand(COUNT) * 2 * math.pi  # t parameter offsets
particle_phases = np.random.rand(COUNT) * 2 * math.pi  # phase offsets for variation

# Wander noise bases for transition phases
base_pos = np.random.rand(COUNT, 3) * GRID
phase_offsets = np.random.rand(COUNT, 3) * 2 * math.pi

# Pre-compute fiber base points on the 2-sphere (Riemann sphere)
fiber_base_points = []
for i in range(FIBER_SAMPLES):
    # Distribute points on the 2-sphere using spherical coordinates
    phi = (i / FIBER_SAMPLES) * 2 * math.pi  # azimuthal angle
    theta = math.acos(1 - 2 * (i + 0.5) / FIBER_SAMPLES)  # polar angle for uniform distribution

    # Convert to Cartesian coordinates on unit sphere
    x = math.sin(theta) * math.cos(phi)
    y = math.sin(theta) * math.sin(phi)
    z = math.cos(theta)

    fiber_base_points.append(np.array([x, y, z]))


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


def butter_smoothstep(edge0: float, edge1: float, x: float) -> float:
    """Extra smooth transition curve for buttery smooth animations"""
    t = max(0.0, min(1.0, (x - edge0) / (edge1 - edge0)))
    # Quintic smoothstep: 6t^5 - 15t^4 + 10t^3 (even smoother than cubic)
    return t * t * t * (t * (t * 6 - 15) + 10)


def lerp(a, b, t):
    return a * (1 - t) + b * t


def hsv_bytes(h, s: float = 1.0, v: float = 1.0):
    r, g, b = hsv_to_rgb(h, s, v)
    return int(r * 255), int(g * 255), int(b * 255)


def quaternion_multiply(q1, q2):
    """Multiply two quaternions represented as [w, x, y, z]"""
    w1, x1, y1, z1 = q1
    w2, x2, y2, z2 = q2
    return np.array([
        w1*w2 - x1*x2 - y1*y2 - z1*z2,
        w1*x2 + x1*w2 + y1*z2 - z1*y2,
        w1*y2 - x1*z2 + y1*w2 + z1*x2,
        w1*z2 + x1*y2 - y1*x2 + z1*w2
    ])


def stereographic_inverse(point_2d):
    """
    Inverse stereographic projection: map point on 2-sphere to quaternion on 3-sphere
    Input: point on unit 2-sphere [x, y, z]
    Output: quaternion [w, x, y, z] on unit 3-sphere
    """
    x, y, z = point_2d
    # Use the inverse stereographic projection from the south pole
    denom = 1 + z
    if abs(denom) < 1e-10:  # Handle singularity at south pole
        return np.array([0.0, 0.0, 0.0, 1.0])

    w = x / denom
    u = y / denom
    v = (1 - z) / denom

    # Normalize to ensure it's on the unit 3-sphere
    norm = math.sqrt(w*w + u*u + v*v + 1)
    return np.array([1/norm, w/norm, u/norm, v/norm])


def hopf_fiber_point(base_point_2d, t_param: float) -> np.ndarray:
    """
    Generate a point on the Hopf fiber corresponding to base_point_2d
    t_param: parameter from 0 to 2π around the fiber circle
    Returns: 3D point projected from the 3-sphere to 3D space
    """
    # Get the base quaternion on the 3-sphere
    base_quat = stereographic_inverse(base_point_2d)

    # Create a quaternion that rotates around the fiber
    # The fiber is parameterized by multiplication with e^(it/2) = cos(t/2) + i*sin(t/2)
    fiber_quat = np.array([math.cos(t_param/2), math.sin(t_param/2), 0.0, 0.0])

    # The point on the fiber is the product of these quaternions
    fiber_point_quat = quaternion_multiply(base_quat, fiber_quat)

    # Project from 3-sphere to 3D space using stereographic projection
    w, x, y, z = fiber_point_quat

    # Stereographic projection from 3-sphere to 3D space
    if abs(w - 1.0) < 1e-10:  # Handle singularity
        return np.array([0.0, 0.0, 0.0])

    denom = 1 - w
    return np.array([x/denom, y/denom, z/denom]) * SCALE


def rotate(vec: np.ndarray, ax: np.ndarray, angle: float) -> np.ndarray:
    """Rotate vec around axis ax by angle (rad) using Rodrigues formula."""
    if np.linalg.norm(ax) == 0:
        return vec
    ax = ax / np.linalg.norm(ax)
    return (
        vec * math.cos(angle)
        + np.cross(ax, vec) * math.sin(angle)
        + ax * np.dot(ax, vec) * (1 - math.cos(angle))
    )


def render_solid_fibers(frame, t, cluster, rot_y, rot_x, rot_z, flow_time):
    """Render the solid Hopf fiber bundle"""
    # Make the structure grow from nothing - control how many fibers appear
    growth_factor = cluster * cluster  # Quadratic growth for smooth emergence
    active_fibers = max(1, int(FIBER_SAMPLES * growth_factor))

    for fiber_idx in range(active_fibers):
        base_point = fiber_base_points[fiber_idx]

        # Add gentle rotation to the base points (smooth from the beginning)
        rotation_strength = cluster * cluster  # Quadratic ramp for extra smoothness
        base_point = rotate(base_point, np.array([0, 1, 0]), flow_time * rotation_strength * 0.3)
        base_point = rotate(base_point, np.array([1, 0, 0]), flow_time * rotation_strength * 0.2)

        # Make each fiber grow along its length - control how much of each fiber is visible
        fiber_growth = cluster * cluster * cluster  # Cubic growth for even smoother emergence
        active_points = max(3, int(POINTS_PER_FIBER * fiber_growth))

        for i in range(active_points):
            t_param = (i / POINTS_PER_FIBER) * 2 * math.pi

            # Add flow motion along the fiber
            flow_t = t_param + flow_time * cluster * 0.8

            # Get point on Hopf fiber
            fiber_pos = hopf_fiber_point(base_point, flow_t)

            # Create thickness layers for solid appearance
            for thickness_layer in range(-THICKNESS//2, THICKNESS//2 + 1):
                # Calculate fiber direction for thickness offset
                epsilon = 0.1
                fiber_pos_plus = hopf_fiber_point(base_point, flow_t + epsilon)
                fiber_pos_minus = hopf_fiber_point(base_point, flow_t - epsilon)

                tangent = fiber_pos_plus - fiber_pos_minus
                if np.linalg.norm(tangent) > 0:
                    tangent = tangent / np.linalg.norm(tangent)
                    # Create perpendicular vector for thickness
                    perp = np.array([-tangent[1], tangent[0], 0])
                    if np.linalg.norm(perp) < 0.1:
                        perp = np.array([0, -tangent[2], tangent[1]])
                    if np.linalg.norm(perp) > 0:
                        perp = perp / np.linalg.norm(perp)
                else:
                    perp = np.array([1, 0, 0])

                # Apply thickness offset
                thick_pos = fiber_pos + perp * thickness_layer * 1.2

                # Apply rotations
                rot_pos = rotate(thick_pos, np.array([0, 1, 0]), rot_y)
                rot_pos = rotate(rot_pos, np.array([1, 0, 0]), rot_x)
                rot_pos = rotate(rot_pos, np.array([0, 0, 1]), rot_z)
                world_pos = CENTER + rot_pos

                x, y, z = world_pos.astype(int)

                if 0 <= x < GRID and 0 <= y < GRID and 0 <= z < GRID:
                    # Beautiful gradient coloring based on fiber and position
                    fiber_hue = (fiber_idx / FIBER_SAMPLES) % 1.0
                    position_hue = (t_param / (2 * math.pi)) * 0.4
                    time_hue = (t * 0.4) % 1.0

                    hue = (fiber_hue + position_hue + time_hue) % 1.0

                    # Enhanced saturation for fibers
                    saturation = 0.8 + 0.2 * math.sin(flow_t * 2 + fiber_idx)
                    saturation = min(1.0, saturation)

                    # Gradient brightness with proper lighting and smooth fade-in
                    base_brightness = 0.9 - abs(thickness_layer) / max(1, THICKNESS//2) * 0.2
                    flow_brightness = 0.9 + 0.1 * math.sin(flow_t * 3 + t * 4 * math.pi)

                    # Smooth fade-in based on cluster value (buttery smooth appearance)
                    fade_in = cluster * cluster * cluster  # Cubic fade for extra smoothness

                    brightness = base_brightness * flow_brightness * fade_in
                    # Never let it get too dark - maintain beautiful colors throughout
                    brightness = max(0.4 + 0.3 * fade_in, min(1.0, brightness))

                    frame.set_voxel(x, y, z, hsv_bytes(hue, saturation, brightness))


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

for f in range(TOTAL_FRAMES):
    t = f / TOTAL_FRAMES  # 0-1 progress along video

    # -------- Buttery smooth phase blend: unordered → ordered → unordered --------
    # Start transition much earlier and make it longer for smoothness
    if t < 0.05:
        cluster = 0.0
    elif t < 0.35:
        cluster = butter_smoothstep(0.05, 0.35, t)
    elif t < 0.75:
        cluster = 1.0
    else:
        cluster = 1.0 - butter_smoothstep(0.75, 0.95, t)

    # Use the same smooth curve for ordered_t but with slight offset for natural feel
    if t < 0.08:
        ordered_t = 0.0
    elif t < 0.38:
        ordered_t = butter_smoothstep(0.08, 0.38, t)
    elif t < 0.72:
        ordered_t = 1.0
    else:
        ordered_t = 1.0 - butter_smoothstep(0.72, 0.92, t)

    # Shape rotation (multiple axes for complex motion)
    rot_y = ordered_t * 2 * math.pi * FIBER_ROTATION_SPEED  # rotation around Y
    rot_x = ordered_t * math.pi * 0.4  # tilt for better visibility
    rot_z = ordered_t * math.pi * 0.2  # slight roll

    # Flow time (continuous throughout animation)
    flow_time = t * FLOW_SPEED * 2 * math.pi

    frame = splv.Frame(GRID, GRID, GRID)

    # Render solid fibers with smooth fade-in (start appearing much earlier)
    if cluster > 0.05:
        render_solid_fibers(frame, t, cluster, rot_y, rot_x, rot_z, flow_time)

    # Render flowing particles - control how many appear as structure grows
    particle_growth = cluster * cluster  # Quadratic growth for particles too
    active_particles = max(10, int(COUNT * particle_growth))

    for i in range(active_particles):
        # -------- Flowing particle on Hopf fiber --------

        fiber_idx = particle_fiber_indices[i]
        base_point = fiber_base_points[fiber_idx]

        # Add gentle rotation to base point (smooth from the beginning)
        rotation_strength = cluster * cluster  # Quadratic ramp for extra smoothness
        base_point = rotate(base_point, np.array([0, 1, 0]), flow_time * rotation_strength * 0.3)
        base_point = rotate(base_point, np.array([1, 0, 0]), flow_time * rotation_strength * 0.2)

        # t parameter flows around the fiber circle
        base_t = (particle_t_offsets[i] + flow_time) % (2 * math.pi)

        # Add gentle wave motion (reduced during organized phase)
        wave_intensity = 1.0 - cluster * 0.7
        wave_offset = 0.2 * wave_intensity * math.sin(3 * base_t + particle_phases[i])
        t_param = base_t + wave_offset

        # Get position on Hopf fiber
        fiber_pos = hopf_fiber_point(base_point, t_param)

        # Apply rotations
        rot_pos = rotate(fiber_pos, np.array([0, 1, 0]), rot_y)
        rot_pos = rotate(rot_pos, np.array([1, 0, 0]), rot_x)
        rot_pos = rotate(rot_pos, np.array([0, 0, 1]), rot_z)
        ordered_pos = CENTER + rot_pos

        # -------- Wander noise for transition phases --------
        wander_amp = 5

        # Add gentle organized motion even during chaotic phase for smoother start
        gentle_organized_offset = (ordered_pos - CENTER) * (cluster * 0.08)  # Very subtle pull towards organization

        npos = base_pos[i] + np.array(
            [
                math.sin(t * 2 * math.pi + phase_offsets[i, 0]) * wander_amp,
                math.cos(t * 2 * math.pi + phase_offsets[i, 1]) * wander_amp,
                math.sin(t * 1.5 * math.pi + phase_offsets[i, 2]) * wander_amp,
            ]
        ) + gentle_organized_offset

        pos = lerp(npos, ordered_pos, cluster)
        x, y, z = pos.astype(int)

        if 0 <= x < GRID and 0 <= y < GRID and 0 <= z < GRID:
            # Particle coloring - brighter and more dynamic than fibers
            fiber_hue = (fiber_idx / FIBER_SAMPLES) % 1.0
            time_hue = (t * 0.6) % 1.0
            particle_hue = 0.3 * math.sin(base_t * 2 + t * 5 * math.pi)

            hue = (fiber_hue + time_hue + particle_hue) % 1.0

            # Higher saturation for particles to make them stand out
            saturation = 0.95 + 0.05 * math.sin(t_param * 4 + t * 6 * math.pi)

            # Brighter particles, especially during transitions - never too dark
            base_particle_brightness = 0.85 + 0.15 * math.sin(wave_offset * 8)
            if cluster < 0.7:  # brighter during transitions
                base_particle_brightness += (0.7 - cluster) * 0.25

            # Ensure particles are always bright and beautiful, even at the beginning
            brightness = max(0.6, min(1.0, base_particle_brightness))

            frame.set_voxel(x, y, z, hsv_bytes(hue, saturation, brightness))

    enc.encode(frame)

    if f % FPS == 0:
        print(f"  second {f // FPS + 1} / {DURATION} - Hopf fibration flowing")

enc.finish()
print("Done. Saved", OUTPUT)
	#!/usr/bin/env python3
	"""
	hopf_fibration.py
	18-second 3D voxel animation of the Hopf fibration - one of the most beautiful structures in mathematics.
	The Hopf fibration maps the 3-sphere to the 2-sphere, creating interlocking circles that never touch.
	Each point on the 2-sphere corresponds to a circle on the 3-sphere, creating a mesmerizing fiber bundle.

	The Hopf fibration demonstrates the deep connection between topology, geometry, and physics,
	and is fundamental to understanding fiber bundles, gauge theory, and quantum mechanics.

	Run:
	pip install spatialstudio numpy
	python hopf_fibration.py

	Outputs:
	hopf_fibration.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 = 18 # seconds
	COUNT = 1000 # base particle count for flow
	OUTPUT = "../outputs/hopf_fibration.splv"
	# -------------------------------------------------

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

	# Hopf fibration parameters
	SCALE = GRID * 0.28 # scale factor for the fibers
	FIBER_SAMPLES = 24 # number of sample fibers to display
	POINTS_PER_FIBER = 60 # points along each fiber circle
	THICKNESS = 2 # thickness for solid appearance

	# Flow parameters
	FLOW_SPEED = 1.8 # how fast particles flow around the fibers
	FIBER_ROTATION_SPEED = 0.5 # how fast the fiber bundle rotates

	# Pre-compute flow parameters for flowing particles
	np.random.seed(2024)
	particle_fiber_indices = np.random.randint(0, FIBER_SAMPLES, COUNT)
	particle_t_offsets = np.random.rand(COUNT) * 2 * math.pi # t parameter offsets
	particle_phases = np.random.rand(COUNT) * 2 * math.pi # phase offsets for variation

	# Wander noise bases for transition phases
	base_pos = np.random.rand(COUNT, 3) * GRID
	phase_offsets = np.random.rand(COUNT, 3) * 2 * math.pi

	# Pre-compute fiber base points on the 2-sphere (Riemann sphere)
	fiber_base_points = []
	for i in range(FIBER_SAMPLES):
	# Distribute points on the 2-sphere using spherical coordinates
	phi = (i / FIBER_SAMPLES) * 2 * math.pi # azimuthal angle
	theta = math.acos(1 - 2 * (i + 0.5) / FIBER_SAMPLES) # polar angle for uniform distribution

	# Convert to Cartesian coordinates on unit sphere
	x = math.sin(theta) * math.cos(phi)
	y = math.sin(theta) * math.sin(phi)
	z = math.cos(theta)

	fiber_base_points.append(np.array([x, y, z]))


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


	def butter_smoothstep(edge0: float, edge1: float, x: float) -> float:
	"""Extra smooth transition curve for buttery smooth animations"""
	t = max(0.0, min(1.0, (x - edge0) / (edge1 - edge0)))
	# Quintic smoothstep: 6t^5 - 15t^4 + 10t^3 (even smoother than cubic)
	return t * t * t * (t * (t * 6 - 15) + 10)


	def lerp(a, b, t):
	return a * (1 - t) + b * t


	def hsv_bytes(h, s: float = 1.0, v: float = 1.0):
	r, g, b = hsv_to_rgb(h, s, v)
	return int(r * 255), int(g * 255), int(b * 255)


	def quaternion_multiply(q1, q2):
	"""Multiply two quaternions represented as [w, x, y, z]"""
	w1, x1, y1, z1 = q1
	w2, x2, y2, z2 = q2
	return np.array([
	w1w2 - x1x2 - y1y2 - z1z2,
	w1x2 + x1w2 + y1z2 - z1y2,
	w1y2 - x1z2 + y1w2 + z1x2,
	w1z2 + x1y2 - y1x2 + z1w2
	])


	def stereographic_inverse(point_2d):
	"""
	Inverse stereographic projection: map point on 2-sphere to quaternion on 3-sphere
	Input: point on unit 2-sphere [x, y, z]
	Output: quaternion [w, x, y, z] on unit 3-sphere
	"""
	x, y, z = point_2d
	# Use the inverse stereographic projection from the south pole
	denom = 1 + z
	if abs(denom) < 1e-10: # Handle singularity at south pole
	return np.array([0.0, 0.0, 0.0, 1.0])

	w = x / denom
	u = y / denom
	v = (1 - z) / denom

	# Normalize to ensure it's on the unit 3-sphere
	norm = math.sqrt(ww + uu + v*v + 1)
	return np.array([1/norm, w/norm, u/norm, v/norm])


	def hopf_fiber_point(base_point_2d, t_param: float) -> np.ndarray:
	"""
	Generate a point on the Hopf fiber corresponding to base_point_2d
	t_param: parameter from 0 to 2π around the fiber circle
	Returns: 3D point projected from the 3-sphere to 3D space
	"""
	# Get the base quaternion on the 3-sphere
	base_quat = stereographic_inverse(base_point_2d)

	# Create a quaternion that rotates around the fiber
	# The fiber is parameterized by multiplication with e^(it/2) = cos(t/2) + i*sin(t/2)
	fiber_quat = np.array([math.cos(t_param/2), math.sin(t_param/2), 0.0, 0.0])

	# The point on the fiber is the product of these quaternions
	fiber_point_quat = quaternion_multiply(base_quat, fiber_quat)

	# Project from 3-sphere to 3D space using stereographic projection
	w, x, y, z = fiber_point_quat

	# Stereographic projection from 3-sphere to 3D space
	if abs(w - 1.0) < 1e-10: # Handle singularity
	return np.array([0.0, 0.0, 0.0])

	denom = 1 - w
	return np.array([x/denom, y/denom, z/denom]) * SCALE


	def rotate(vec: np.ndarray, ax: np.ndarray, angle: float) -> np.ndarray:
	"""Rotate vec around axis ax by angle (rad) using Rodrigues formula."""
	if np.linalg.norm(ax) == 0:
	return vec
	ax = ax / np.linalg.norm(ax)
	return (
	vec * math.cos(angle)
	+ np.cross(ax, vec) * math.sin(angle)
	+ ax * np.dot(ax, vec) * (1 - math.cos(angle))
	)


	def render_solid_fibers(frame, t, cluster, rot_y, rot_x, rot_z, flow_time):
	"""Render the solid Hopf fiber bundle"""
	# Make the structure grow from nothing - control how many fibers appear
	growth_factor = cluster * cluster # Quadratic growth for smooth emergence
	active_fibers = max(1, int(FIBER_SAMPLES * growth_factor))

	for fiber_idx in range(active_fibers):
	base_point = fiber_base_points[fiber_idx]

	# Add gentle rotation to the base points (smooth from the beginning)
	rotation_strength = cluster * cluster # Quadratic ramp for extra smoothness
	base_point = rotate(base_point, np.array([0, 1, 0]), flow_time * rotation_strength * 0.3)
	base_point = rotate(base_point, np.array([1, 0, 0]), flow_time * rotation_strength * 0.2)

	# Make each fiber grow along its length - control how much of each fiber is visible
	fiber_growth = cluster * cluster * cluster # Cubic growth for even smoother emergence
	active_points = max(3, int(POINTS_PER_FIBER * fiber_growth))

	for i in range(active_points):
	t_param = (i / POINTS_PER_FIBER) * 2 * math.pi

	# Add flow motion along the fiber
	flow_t = t_param + flow_time * cluster * 0.8

	# Get point on Hopf fiber
	fiber_pos = hopf_fiber_point(base_point, flow_t)

	# Create thickness layers for solid appearance
	for thickness_layer in range(-THICKNESS//2, THICKNESS//2 + 1):
	# Calculate fiber direction for thickness offset
	epsilon = 0.1
	fiber_pos_plus = hopf_fiber_point(base_point, flow_t + epsilon)
	fiber_pos_minus = hopf_fiber_point(base_point, flow_t - epsilon)

	tangent = fiber_pos_plus - fiber_pos_minus
	if np.linalg.norm(tangent) > 0:
	tangent = tangent / np.linalg.norm(tangent)
	# Create perpendicular vector for thickness
	perp = np.array([-tangent[1], tangent[0], 0])
	if np.linalg.norm(perp) < 0.1:
	perp = np.array([0, -tangent[2], tangent[1]])
	if np.linalg.norm(perp) > 0:
	perp = perp / np.linalg.norm(perp)
	else:
	perp = np.array([1, 0, 0])

	# Apply thickness offset
	thick_pos = fiber_pos + perp * thickness_layer * 1.2

	# Apply rotations
	rot_pos = rotate(thick_pos, np.array([0, 1, 0]), rot_y)
	rot_pos = rotate(rot_pos, np.array([1, 0, 0]), rot_x)
	rot_pos = rotate(rot_pos, np.array([0, 0, 1]), rot_z)
	world_pos = CENTER + rot_pos

	x, y, z = world_pos.astype(int)

	if 0 <= x < GRID and 0 <= y < GRID and 0 <= z < GRID:
	# Beautiful gradient coloring based on fiber and position
	fiber_hue = (fiber_idx / FIBER_SAMPLES) % 1.0
	position_hue = (t_param / (2 * math.pi)) * 0.4
	time_hue = (t * 0.4) % 1.0

	hue = (fiber_hue + position_hue + time_hue) % 1.0

	# Enhanced saturation for fibers
	saturation = 0.8 + 0.2 * math.sin(flow_t * 2 + fiber_idx)
	saturation = min(1.0, saturation)

	# Gradient brightness with proper lighting and smooth fade-in
	base_brightness = 0.9 - abs(thickness_layer) / max(1, THICKNESS//2) * 0.2
	flow_brightness = 0.9 + 0.1 * math.sin(flow_t * 3 + t * 4 * math.pi)

	# Smooth fade-in based on cluster value (buttery smooth appearance)
	fade_in = cluster * cluster * cluster # Cubic fade for extra smoothness

	brightness = base_brightness * flow_brightness * fade_in
	# Never let it get too dark - maintain beautiful colors throughout
	brightness = max(0.4 + 0.3 * fade_in, min(1.0, brightness))

	frame.set_voxel(x, y, z, hsv_bytes(hue, saturation, brightness))


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

	for f in range(TOTAL_FRAMES):
	t = f / TOTAL_FRAMES # 0-1 progress along video

	# -------- Buttery smooth phase blend: unordered → ordered → unordered --------
	# Start transition much earlier and make it longer for smoothness
	if t < 0.05:
	cluster = 0.0
	elif t < 0.35:
	cluster = butter_smoothstep(0.05, 0.35, t)
	elif t < 0.75:
	cluster = 1.0
	else:
	cluster = 1.0 - butter_smoothstep(0.75, 0.95, t)

	# Use the same smooth curve for ordered_t but with slight offset for natural feel
	if t < 0.08:
	ordered_t = 0.0
	elif t < 0.38:
	ordered_t = butter_smoothstep(0.08, 0.38, t)
	elif t < 0.72:
	ordered_t = 1.0
	else:
	ordered_t = 1.0 - butter_smoothstep(0.72, 0.92, t)

	# Shape rotation (multiple axes for complex motion)
	rot_y = ordered_t * 2 * math.pi * FIBER_ROTATION_SPEED # rotation around Y
	rot_x = ordered_t * math.pi * 0.4 # tilt for better visibility
	rot_z = ordered_t * math.pi * 0.2 # slight roll

	# Flow time (continuous throughout animation)
	flow_time = t * FLOW_SPEED * 2 * math.pi

	frame = splv.Frame(GRID, GRID, GRID)

	# Render solid fibers with smooth fade-in (start appearing much earlier)
	if cluster > 0.05:
	render_solid_fibers(frame, t, cluster, rot_y, rot_x, rot_z, flow_time)

	# Render flowing particles - control how many appear as structure grows
	particle_growth = cluster * cluster # Quadratic growth for particles too
	active_particles = max(10, int(COUNT * particle_growth))

	for i in range(active_particles):
	# -------- Flowing particle on Hopf fiber --------

	fiber_idx = particle_fiber_indices[i]
	base_point = fiber_base_points[fiber_idx]

	# Add gentle rotation to base point (smooth from the beginning)
	rotation_strength = cluster * cluster # Quadratic ramp for extra smoothness
	base_point = rotate(base_point, np.array([0, 1, 0]), flow_time * rotation_strength * 0.3)
	base_point = rotate(base_point, np.array([1, 0, 0]), flow_time * rotation_strength * 0.2)

	# t parameter flows around the fiber circle
	base_t = (particle_t_offsets[i] + flow_time) % (2 * math.pi)

	# Add gentle wave motion (reduced during organized phase)
	wave_intensity = 1.0 - cluster * 0.7
	wave_offset = 0.2 * wave_intensity * math.sin(3 * base_t + particle_phases[i])
	t_param = base_t + wave_offset

	# Get position on Hopf fiber
	fiber_pos = hopf_fiber_point(base_point, t_param)

	# Apply rotations
	rot_pos = rotate(fiber_pos, np.array([0, 1, 0]), rot_y)
	rot_pos = rotate(rot_pos, np.array([1, 0, 0]), rot_x)
	rot_pos = rotate(rot_pos, np.array([0, 0, 1]), rot_z)
	ordered_pos = CENTER + rot_pos

	# -------- Wander noise for transition phases --------
	wander_amp = 5

	# Add gentle organized motion even during chaotic phase for smoother start
	gentle_organized_offset = (ordered_pos - CENTER) * (cluster * 0.08) # Very subtle pull towards organization

	npos = base_pos[i] + np.array(
	[
	math.sin(t * 2 * math.pi + phase_offsets[i, 0]) * wander_amp,
	math.cos(t * 2 * math.pi + phase_offsets[i, 1]) * wander_amp,
	math.sin(t * 1.5 * math.pi + phase_offsets[i, 2]) * wander_amp,
	]
	) + gentle_organized_offset

	pos = lerp(npos, ordered_pos, cluster)
	x, y, z = pos.astype(int)

	if 0 <= x < GRID and 0 <= y < GRID and 0 <= z < GRID:
	# Particle coloring - brighter and more dynamic than fibers
	fiber_hue = (fiber_idx / FIBER_SAMPLES) % 1.0
	time_hue = (t * 0.6) % 1.0
	particle_hue = 0.3 * math.sin(base_t * 2 + t * 5 * math.pi)

	hue = (fiber_hue + time_hue + particle_hue) % 1.0

	# Higher saturation for particles to make them stand out
	saturation = 0.95 + 0.05 * math.sin(t_param * 4 + t * 6 * math.pi)

	# Brighter particles, especially during transitions - never too dark
	base_particle_brightness = 0.85 + 0.15 * math.sin(wave_offset * 8)
	if cluster < 0.7: # brighter during transitions
	base_particle_brightness += (0.7 - cluster) * 0.25

	# Ensure particles are always bright and beautiful, even at the beginning
	brightness = max(0.6, min(1.0, base_particle_brightness))

	frame.set_voxel(x, y, z, hsv_bytes(hue, saturation, brightness))

	enc.encode(frame)

	if f % FPS == 0:
	print(f" second {f // FPS + 1} / {DURATION} - Hopf fibration flowing")

	enc.finish()
	print("Done. Saved", OUTPUT)
No results found