first formulation

#!/usr/bin/env python3

###############################################################################
# Generative animal-like animation in Python
#
# Optimised MP4 version of the compact JavaScript / p5.js formula:
#
# a=(m,d=mag(k=9*cos(i/81),e=i/765-13)/4)=>point(
#   (q=79-2*sin(k*3)+
#      sin(k*k<19?t*3+d*4:d/2+4)/2*k*
#      (9+5*sin(d*d-e/6-t+m))
#   )*sin(c=d*d/9-t/16+m)+200,
#   (q+50)*cos(c)+200
# )
#
# t=0,draw=$=>{
#   t||createCanvas(w=400,w);
#   background(9).stroke(w,96);
#   for(t+=PI/45,i=2e4;i--;)a(i%2*9)
# }
#
# Main optimisation:
#
#   - save MP4 directly with FFMpegWriter instead of GIF with PillowWriter
#   - precompute immutable arrays
#   - update one existing scatter object frame by frame
###############################################################################

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation, FFMpegWriter

# =============================================================================
# 1. Tunable parameters
# =============================================================================

n_points = 20_000

# 540 * 4 frames at 25 fps gives:
#
#   2160 / 25 = 86.4 seconds
n_frames = 540

fps = 25

canvas_width = 400
canvas_height = 400

# Rendered output size.
#
# Final pixel size is:
#
#   figsize_inches * dpi
#
# Since figsize is computed as pixels / dpi below, this gives 800 x 800.
render_width_px = 600
render_height_px = 600
dpi = 150

# Visual tuning.
#
# Matplotlib scatter "s" is marker area in points squared.
point_size = 1.5
point_alpha = 0.45

background_colour = "#090909"

output_file = "living_animal_second_formula_python_thicker.mp4"

# FFmpeg encoding settings.
#
# crf:
#   lower = better quality / larger file
#   18 is visually high quality
#
# preset:
#   faster presets encode faster but may produce larger files
ffmpeg_crf = 18
ffmpeg_preset = "veryfast"
ffmpeg_threads = 4

# =============================================================================
# 2. Static point data
# =============================================================================

i = np.arange(n_points)

# JavaScript:
#
#   k = 9 * cos(i / 81)
k = 9 * np.cos(i / 81)

# JavaScript:
#
#   e = i / 765 - 13
e = i / 765 - 13

# JavaScript:
#
#   d = mag(k, e) / 4
#
# In p5.js:
#
#   mag(k, e) = sqrt(k^2 + e^2)
d = np.sqrt(k**2 + e**2) / 4

# JavaScript:
#
#   m = i % 2 * 9
m = (i % 2) * 9

# =============================================================================
# 3. Precomputed static terms
# =============================================================================

# These quantities do not change with frame/time.
# Precomputing them avoids repeated work inside compute_frame().
k2 = k**2
d2 = d**2

phase_mask = k2 < 19
phase_else = d / 2 + 4
d_times_4 = d * 4

q_static = 79 - 2 * np.sin(k * 3)
q_multiplier = k / 2

wave_static = d2 - e / 6 + m
c_static = d2 / 9 + m

# Preallocated coordinate array for matplotlib scatter offsets.
#
# Matplotlib expects an N x 2 array:
#
#   [[x1, y1],
#    [x2, y2],
#    ...]
offsets = np.empty((n_points, 2), dtype=np.float64)

# Preallocated phase array.
phase_switch = np.empty(n_points, dtype=np.float64)

# =============================================================================
# 4. Frame computation
# =============================================================================

def compute_frame(frame_id):
    """
    Compute one frame of the animation.

    Parameters
    ----------
    frame_id : int
        Integer frame number: 1, 2, 3, ...

        The frame number is converted to continuous time as:

            t = frame_id * pi / 45

        matching the original JavaScript update:

            t += PI / 45

    Returns
    -------
    offsets : numpy.ndarray
        N x 2 array of x/y coordinates.
    """

    t = frame_id * np.pi / 45

    # JavaScript ternary:
    #
    #   k*k < 19 ? t*3 + d*4 : d/2 + 4
    #
    # Instead of np.where(), fill a preallocated array.
    phase_switch[:] = phase_else
    phase_switch[phase_mask] = t * 3 + d_times_4[phase_mask]

    # Main deformation term:
    #
    #   q = 79 - 2*sin(k*3) +
    #       sin(phase_switch)/2 * k *
    #       (9 + 5*sin(d*d - e/6 - t + m))
    q = (
        q_static
        + np.sin(phase_switch)
        * q_multiplier
        * (9 + 5 * np.sin(wave_static - t))
    )

    # Angular term:
    #
    #   c = d*d/9 - t/16 + m
    c = c_static - t / 16

    # Final coordinates:
    #
    #   x = q * sin(c) + 200
    #   y = (q + 50) * cos(c) + 200
    offsets[:, 0] = q * np.sin(c) + 200
    offsets[:, 1] = (q + 50) * np.cos(c) + 200

    return offsets

# =============================================================================
# 5. Matplotlib figure setup
# =============================================================================

fig_width_inches = render_width_px / dpi
fig_height_inches = render_height_px / dpi

fig, ax = plt.subplots(
    figsize=(fig_width_inches, fig_height_inches),
    dpi=dpi
)

fig.patch.set_facecolor(background_colour)
ax.set_facecolor(background_colour)

ax.set_xlim(0, canvas_width)
ax.set_ylim(0, canvas_height)
ax.set_aspect("equal")

ax.set_xticks([])
ax.set_yticks([])

for spine in ax.spines.values():
    spine.set_visible(False)

ax.margins(0)
plt.subplots_adjust(left=0, right=1, top=1, bottom=0)

# =============================================================================
# 6. Initial frame and scatter object
# =============================================================================

initial_offsets = compute_frame(1).copy()

points = ax.scatter(
    initial_offsets[:, 0],
    initial_offsets[:, 1],
    s=point_size,
    c="white",
    alpha=point_alpha,
    linewidths=0
)

# =============================================================================
# 7. Animation update function
# =============================================================================

def update(frame_id):
    current_offsets = compute_frame(frame_id)

    # Update the existing scatter object.
    # This avoids rebuilding the plot.
    points.set_offsets(current_offsets)

    return (points,)

# =============================================================================
# 8. Build and save animation directly as MP4
# =============================================================================

animation = FuncAnimation(
    fig,
    update,
    frames=np.arange(1, n_frames + 1),
    interval=1000 / fps,
    blit=True
)

writer = FFMpegWriter(
    fps=fps,
    codec="libx264",
    bitrate=-1,
    extra_args=[
        "-preset", ffmpeg_preset,
        "-crf", str(ffmpeg_crf),
        "-threads", str(ffmpeg_threads),
        "-pix_fmt", "yuv420p",
    ],
)

animation.save(
    output_file,
    writer=writer
)

plt.close(fig)

print(f"Saved: {output_file}")