###############################################################################
# Living Animal Animation in R
#
# Translation/adaptation of the Wolfram Mathematica animation from:
# https://community.wolfram.com/groups/-/m/t/3713047
#
# This script:
#   1. Creates a cloud of 10,000 points
#   2. Applies nonlinear trigonometric transformations
#   3. Recomputes the point positions over time
#   4. Separates points into base/body points and highlighted tip points
#   5. Changes the colour of the highlighted tips globally over time
#   6. Exports the result as an animated GIF
#
# Main R packages:
#   - ggplot2   : plotting
#   - gganimate : animation
#   - gifski    : GIF rendering backend
#   - dplyr     : readable data manipulation
#   - purrr     : readable iteration over frames
###############################################################################

# -----------------------------------------------------------------------------
# Load required packages
# -----------------------------------------------------------------------------

library(ggplot2)
library(gganimate)
library(gifski)
library(dplyr)
library(purrr)

# -----------------------------------------------------------------------------
# Animation parameters
#
# n_frames:
#   Number of frames in the animation.
#
# fps:
#   Frames per second.
#
# Higher fps makes motion smoother but sometimes less organic.
# You found that fps = 25 looks better than fps = 30.
#
# duration in seconds is approximately:
#
#   n_frames / fps
#
# With 180 frames and 25 fps:
#
#   180 / 25 = 7.2 seconds
# -----------------------------------------------------------------------------

n_frames <- 180
fps <- 25

output_file <- "living_animal.gif"

# -----------------------------------------------------------------------------
# Generate the base point index
#
# Mathematica equivalent:
#
#   x = N @ Range[0, 9999]
#
# These are not yet the final x-coordinates.
# They are an index used to generate the point cloud.
# -----------------------------------------------------------------------------

x <- seq(0, 9999)

# -----------------------------------------------------------------------------
# Precompute quantities that do not depend on time
#
# These values are static, so we compute them once and reuse them for every frame.
#
# Mathematica logic:
#
#   y = x / 235
#   k = 4 Cos[x / 21]
#   e = y / 8 - 20
#   d = Sqrt[k^2 + e^2]
#
# Meaning of the variables:
#
#   y_index:
#     auxiliary scaled index
#
#   k:
#     oscillatory component
#
#   e:
#     slow linear component
#
#   d:
#     distance-like quantity combining k and e
#
#   group:
#     separates highlighted "tip" points from the rest of the body
# -----------------------------------------------------------------------------

base_data <- tibble(
  x_index = x,

  # Auxiliary variable used in the original formula
  y_index = x / 235,

  # Oscillatory component
  k = 4 * cos(x / 21)
) |>
  mutate(
    # Slow linear component
    e = y_index / 8 - 20,

    # Distance-like quantity
    d = sqrt(k^2 + e^2),

    # Equivalent to the Mathematica selection of highlighted points:
    #
    #   k^2 > 15
    #
    # These points are drawn larger and with a colour that changes over time.
    group = if_else(k^2 > 15, "tips", "base")
  )

# -----------------------------------------------------------------------------
# Function generating one animation frame
#
# Inputs:
#
#   t:
#     current animation time
#
#   frame_id:
#     integer frame number used by gganimate
#
#   data:
#     the static base data created above
#
# Output:
#
#   A tibble containing final x/y coordinates, group, colour, and frame number.
#
# The main time-dependent formula is q.
# It controls much of the organic deformation of the shape.
# -----------------------------------------------------------------------------

make_frame <- function(t, frame_id, data = base_data) {

  # ---------------------------------------------------------------------------
  # Global colour oscillation
  #
  # Original Mathematica colour expression:
  #
  #   Blend[{White, Red}, Sin[t]^2]
  #
  # This means:
  #
  #   sin(t)^2 = 0  -> white
  #   sin(t)^2 = 1  -> red
  #
  # In RGB terms:
  #
  #   white = rgb(1, 1, 1)
  #   red   = rgb(1, 0, 0)
  #
  # So the blend can be represented as:
  #
  #   rgb(1, 1 - intensity, 1 - intensity)
  #
  # where:
  #
  #   intensity = sin(t)^2
  # ---------------------------------------------------------------------------

  red_intensity <- sin(t)^2

  tips_colour <- rgb(
    red = 1,
    green = 1 - red_intensity,
    blue = 1 - red_intensity
  )

  # ---------------------------------------------------------------------------
  # Compute the frame-specific coordinates
  # ---------------------------------------------------------------------------

  data |>
    mutate(
      # Equivalent to Mathematica:
      #
      #   c = d - t
      #
      # This creates a rotating phase shift over time.
      c = d - t,

      # Main nonlinear deformation term.
      #
      # Mathematica equivalent:
      #
      #   q = 3 Sin[2 k] + 0.3/k +
      #       Sin[y/19] k (9 + 2 Sin[14 e - 3 d + 2 t])
      #
      # This combines:
      #
      #   - local oscillations through sin(2 * k)
      #   - a small inverse-k term
      #   - a larger time-varying deformation
      q =
        3 * sin(2 * k) +
        0.3 / k +
        sin(y_index / 19) * k *
        (9 + 2 * sin(14 * e - 3 * d + 2 * t)),

      # Final x coordinate.
      #
      # Mathematica equivalent:
      #
      #   q + 50 Cos[c] + 200
      #
      # The +200 recentres the image horizontally.
      x = q + 50 * cos(c) + 200,

      # Final y coordinate.
      #
      # Mathematica equivalent:
      #
      #   400 - (q Sin[c] + 39 d - 475)
      #
      # This simplifies to:
      #
      #   875 - q Sin[c] - 39 d
      #
      # The expression also flips/repositions the vertical direction.
      y = 875 - q * sin(c) - 39 * d,

      # Frame identifier used by gganimate.
      frame = frame_id,

      # Colour for each point.
      #
      # Base/body points remain white.
      # Tip points use the globally changing colour.
      colour = if_else(group == "tips", tips_colour, "white")
    ) |>
    select(x, y, group, colour, frame)
}

# -----------------------------------------------------------------------------
# Define the animation timeline
#
# The original Mathematica animation is usually parameterised over 0 to 2*pi.
#
# Here we use 0 to 6*pi, which gives a longer animation and makes the movement
# easier to appreciate.
#
# To make it closer to the original, use:
#
#   times <- seq(0, 2 * pi, length.out = n_frames)
# -----------------------------------------------------------------------------

times <- seq(
  from = 0,
  to = 6 * pi,
  length.out = n_frames
)

# -----------------------------------------------------------------------------
# Generate the complete animation data set
#
# map2_dfr() iterates over two inputs at once:
#
#   1. times
#   2. frame numbers
#
# For each pair, it calls:
#
#   make_frame(time, frame_number)
#
# The "_dfr" suffix means that all resulting data frames are row-bound into one
# large data frame.
#
# The final data frame has:
#
#   number of rows = number of points * number of frames
#
# In this script:
#
#   10,000 * 180 = 1,800,000 rows
# -----------------------------------------------------------------------------

df <- map2_dfr(
  times,
  seq_along(times),
  make_frame
)

# -----------------------------------------------------------------------------
# Build the animated plot
#
# We draw two layers:
#
#   1. highlighted tips
#      - larger
#      - semi-transparent
#      - colour oscillates between white and red
#
#   2. base/body points
#      - smaller
#      - white
#
# The colour column already contains literal colour names or RGB hex values.
# Therefore, we use scale_colour_identity().
# -----------------------------------------------------------------------------

p <- ggplot(df, aes(x, y)) +

  # ---------------------------------------------------------------------------
  # Highlighted tip points
  # ---------------------------------------------------------------------------

  geom_point(
    data = filter(df, group == "tips"),
    aes(colour = colour),
    alpha = 0.5,
    size = 2.2
  ) +

  # ---------------------------------------------------------------------------
  # Base/body points
  # ---------------------------------------------------------------------------

  geom_point(
    data = filter(df, group == "base"),
    aes(colour = colour),
    alpha = 0.75,
    size = 0.45
  ) +

  # Use colour values exactly as stored in the data
  scale_colour_identity() +

  # ---------------------------------------------------------------------------
  # Visible plotting window
  #
  # These limits reproduce the cropped view of the original animation.
  #
  # expand = FALSE removes extra padding around the plot.
  # ---------------------------------------------------------------------------

  coord_cartesian(
    xlim = c(100, 300),
    ylim = c(75, 320),
    expand = FALSE
  ) +

  # Remove axes, ticks, labels and grid lines
  theme_void() +

  # Black background
  theme(
    plot.background = element_rect(fill = "black", colour = NA),
    panel.background = element_rect(fill = "black", colour = NA)
  ) +

  # ---------------------------------------------------------------------------
  # Animation instruction
  #
  # transition_manual(frame) tells gganimate:
  #
  #   use each distinct value of "frame" as one animation frame
  # ---------------------------------------------------------------------------

  transition_manual(frame)

# -----------------------------------------------------------------------------
# Render the animation
#
# This writes the animated GIF to the current working directory.
#
# You can check your working directory with:
#
#   getwd()
#
# You can change it with:
#
#   setwd("path/to/folder")
# -----------------------------------------------------------------------------

animate(
  p,
  nframes = n_frames,
  fps = fps,
  width = 550,
  height = 550,
  renderer = gifski_renderer(output_file)
)

###############################################################################
# Output
#
# The file:
#
#   living_animal.gif
#
# is created in the current working directory.
#
# On Linux, open it with:
#
#   xdg-open living_animal.gif
#
# or:
#
#   firefox living_animal.gif
#
###############################################################################