Packages

install.packages(c(
  "tidyverse",
  "ggplot2",
  "ggforce",
  "ambient",
  "colorspace",
  "gifski",
  "tuneR",
  "png",
  "scales"
))
install.packages(c("ggplot2", "ggforce", "ambient", "colorspace"))
install.packages(c("ggplot2", "ggforce", "ambient", "colorspace"))
install.packages("deldir")
Kubrickplot7

n <- 1000  # finer resolution
x <- seq(-4, 4, length.out = n)
y <- seq(-4, 4, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Central void (isolation)
void <- -1.2 * exp(-r^2)

# Multi-layered grief waves
wave1 <- sin(6 * r + theta) * exp(-0.8 * r)
wave2 <- 0.6 * sin(10 * r - 2*theta) * exp(-1.2 * r)
wave3 <- 0.4 * sin(14 * r + 4*theta) * exp(-1.6 * r)

#  Stronger asymmetry (loss)
asym <- 0.5 * sin(4 * theta) * exp(-0.6 * r)

#  Finer fractal noise (chaotic emotion)
noise <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 6,
  x = grid$x,
  y = grid$y,
  frequency = 1.5
) * 0.4

# Combine layers
grid$z <- void + wave1 + wave2 + wave3 + asym + noise

# Color palette: dark blues, purples, and subtle red hints
colors <- diverging_hcl(150, palette = "Blue-Red 3")
colors <- rev(colors)  # make center darker

# Base plot
p <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))

# Optional: ghost Bezier curves for fading memories

spirals <- data.frame(
  x = c(0, 0.5, 1, 1.5, 2),
  y = c(0, 0.3, 0.6, 0.9, 1),
  xend = c(0.2, 0.7, 1.2, 1.7, 2.2),
  yend = c(0.2, 0.5, 0.8, 1.1, 1.3)
)

p + geom_bezier(aes(x = x, y = y, xend = xend, yend = yend),
                data = spirals, color = "#ffffff20", size = 0.3, inherit.aes = FALSE)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning in geom_bezier(aes(x = x, y = y, xend = xend, yend = yend), data =
## spirals, : Ignoring unknown aesthetics: xend and yend
## Warning: Computation failed in `stat_bezier()`.
## Caused by error in `compute_panel()`:
## ! Only support for quadratic and cubic beziers
## ℹ Make sure each group consists of 3 or 4 rows

p

library(ggplot2)
library(ggforce)


# Grid
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Stronger angular distortion (clearly visible)
set.seed(42)
theta_distort <- theta + 0.25 * sin(3 * r)

# More dynamic spiral layering
wave1 <- sin(6 * r + 4 * theta_distort) * exp(-0.4 * r)
wave2 <- 0.8 * sin(10 * r - 5 * theta_distort) * exp(-0.7 * r)
wave3 <- 0.6 * sin(18 * r + 8 * theta_distort) * exp(-1.1 * r)

# Noise weighted toward center
noise <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 3,
  x = grid$x,
  y = grid$y,
  frequency = 1
)

noise <- noise * exp(-0.3 * r)

# Much stronger central sphere
sphere <- 2.5 * exp(-r^2 * 2)

# Radial darkening toward edges (adds depth)
vignette <- -0.8 * (r / max(r))^2

# Combine
z <- wave1 + wave2 + wave3 + noise + sphere + vignette

# Normalize values for stronger contrast
z <- (z - min(z)) / (max(z) - min(z))

grid$z <- z

# More dramatic gradient
colors <- c(
  "#0f1c2e",   # deep blue
  "#2e4057",   # muted blue
  "#5d6d7e",   # cool grey-blue
  "#9b59b6",   # muted purple
  "#f39c12"    # orange core
)

Plot2 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)

# Memory-efficient grid
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Add subtle random perturbation to angles for organic feel
set.seed(42)
theta_distort <- theta + runif(length(theta), -0.1, 0.1)

# Spiral layers
wave1 <- sin(8 * r + 5 * theta_distort) * exp(-0.5 * r)
wave2 <- 0.6 * sin(12 * r - 3 * theta_distort) * exp(-0.8 * r)
wave3 <- 0.4 * sin(16 * r + 7 * theta_distort) * exp(-1.0 * r)

# Fractal noise stronger near the center
noise_base <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 4,
  x = grid$x,
  y = grid$y,
  frequency = 1.2
)
# weight noise by radial decay
noise <- noise_base * exp(-0.5 * r)

# Central orange sphere (larger and softer glow)
sphere <- exp(-r^2 * 3)  # adjust multiplier for glow size

# Combine layers numerically
grid$z <- wave1 + wave2 + wave3 + noise + sphere

# Gradient palette (blues/purples fading into orange near center)
colors <- c("#2e4057", "#5d6d7e", "#7f8c8d", "#f39c12")

# Plot
p2 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
p2

Plot2

library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)

# Grid (memory-safe)
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Angular distortion (structured, not random)
theta_distort <- theta + 0.3 * sin(2 * r)

# Spiral equation
spiral_field <- 10 * r + 6 * theta_distort

# Create thin ridge lines using cosine sharpening
ridge <- cos(spiral_field)

# Sharpen ridges (higher power = thinner lines)
ridge_sharp <- abs(ridge)^20

# Fade outward slightly
ridge_sharp <- ridge_sharp * exp(-0.25 * r)

# Add faint secondary counter-rotation for tension
counter <- abs(cos(8 * r - 5 * theta))^18
counter <- 0.4 * counter * exp(-0.4 * r)

# Strong central sphere
sphere <- 2.5 * exp(-r^2 * 2)

# Darken edges for depth
vignette <- -0.6 * (r / max(r))^2

# Combine
z <- ridge_sharp + counter + sphere + vignette

# Normalize (important for contrast)
z <- (z - min(z)) / (max(z) - min(z))

grid$z <- z

# Controlled palette (ink-like)
colors <- c(
  "#0b132b",  # deep navy
  "#1c2541",  # muted blue
  "#3a506b",  # steel blue
  "#5d6d7e",  # cool grey-blue
  "#9b59b6",  # muted purple
  "#f39c12"   # orange core
)

P3 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
P3

library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)

# Grid (memory-safe)
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Angular distortion (structured, not random)
theta_distort <- theta + 0.3 * sin(2 * r)

# Spiral equation
spiral_field <- 10 * r + 6 * theta_distort

# Create thin ridge lines using cosine sharpening
ridge <- cos(spiral_field)

# Sharpen ridges (higher power = thinner lines)
ridge_sharp <- abs(ridge)^15

# Fade outward slightly
ridge_sharp <- ridge_sharp * exp(-0.25 * r)

# Add faint secondary counter-rotation for tension
counter <- abs(cos(8 * r - 5 * theta))^18
counter <- 0.4 * counter * exp(-0.4 * r)

# Strong central sphere
sphere <- 2.5 * exp(-r^2 * 2)

# Darken edges for depth
vignette <- -0.6 * (r / max(r))^2

# Combine
z <- ridge_sharp + counter + sphere + vignette

# Normalize (important for contrast)
z <- (z - min(z)) / (max(z) - min(z))

grid$z <- z

# Controlled palette (ink-like)
colors <- c(
  "#0b132b",  # deep navy
  "#1c2541",  # muted blue
  "#3a506b",  # steel blue
  "#5d6d7e",  # cool grey-blue
  "#f39c12"   # orange core
)

P4 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
P4

library(ggplot2)
library(deldir)
## deldir 2.0-4      Nickname: "Idol Comparison"
## 
##      The syntax of deldir() has changed since version 
##      0.0-10.  Read the help!!!.
set.seed(42)

# Number of colony seeds (increase for denser mould)
n_points <- 80

# Random seed locations
x <- runif(n_points, -5, 5)
y <- runif(n_points, -5, 5)

# Compute Voronoi tessellation
vor <- deldir(x, y)
tiles <- tile.list(vor)

# Convert to dataframe for ggplot
mould_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
  data.frame(
    x = tiles[[i]]$x,
    y = tiles[[i]]$y,
    id = i
  )
}))

# Add organic distortion
mould_df$x <- mould_df$x + 0.2 * sin(3 * mould_df$y)
mould_df$y <- mould_df$y + 0.2 * cos(3 * mould_df$x)

mould_df$x <- mould_df$x + runif(nrow(mould_df), -0.1, 0.1)
mould_df$y <- mould_df$y + runif(nrow(mould_df), -0.1, 0.1)

n_points <- 150
# Plot
Mould1 <- ggplot(mould_df, aes(x, y, group = id)) +
  geom_polygon(
    fill = "#4b5d3a", 
    color = "#1f2a1f",
    alpha = 0.85
  ) +
  coord_fixed(xlim = c(-5,5), ylim = c(-5,5)) +
  theme_void() +
  theme(
    plot.background = element_rect(fill = "#0d140d", color = NA)
  )
library(ggplot2)

n <- 250

Du <- 0.16
Dv <- 0.08
F  <- 0.060
k  <- 0.062
dt <- 1

# Initialize
U <- matrix(1, n, n)
V <- matrix(0, n, n)

# Add random noise across grid (break symmetry)
set.seed(1)
V <- V + matrix(runif(n*n, 0, 0.05), n, n)

# Seed larger irregular region in center
center <- (n/2 - 20):(n/2 + 20)
U[center, center] <- 0.50
V[center, center] <- 0.25

laplace <- function(M) {
  -4*M +
    rbind(M[-1,], M[n,]) +
    rbind(M[1,], M[-n,]) +
    cbind(M[,-1], M[,n]) +
    cbind(M[,1], M[,-n])
}

steps <- 2000  # important increase

for (i in 1:steps) {
  Lu <- laplace(U)
  Lv <- laplace(V)
  
  uvv <- U * V * V
  
  U <- U + (Du * Lu - uvv + F * (1 - U)) * dt
  V <- V + (Dv * Lv + uvv - (F + k) * V) * dt
}

df <- expand.grid(x = 1:n, y = 1:n)
df$z <- as.vector(V)

df$z <- (df$z - min(df$z)) / (max(df$z) - min(df$z))

colors <- c(
  "#050a05",
  "#0f1f0f",
  "#1f3b1f",
  "#3f6f3f",
  "#7fbf7f",
  "#e6ffe6"
)

ggplot(df, aes(x, y, fill = z)) +
  geom_raster() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))

library(ggplot2)

n <- 250

Du <- 0.16
Dv <- 0.08
F  <- 0.048
k  <- 0.060
dt <- 1

U <- matrix(1, n, n)
V <- matrix(0, n, n)

set.seed(2)

# Stronger noise everywhere
V <- V + matrix(runif(n*n, 0, 0.08), n, n)

# Seed infection in top-left corner (not center)
seed <- 1:40
U[seed, seed] <- 0.50
V[seed, seed] <- 0.25
U[200:230, 150:180] <- 0.50
V[200:230, 150:180] <- 0.25

laplace <- function(M) {
  -4*M +
    rbind(M[-1,], M[n,]) +
    rbind(M[1,], M[-n,]) +
    cbind(M[,-1], M[,n]) +
    cbind(M[,1], M[,-n])
}

steps <- 3500

for (i in 1:steps) {
  Lu <- laplace(U)
  Lv <- laplace(V)
  
  uvv <- U * V * V
  
  U <- U + (Du * Lu - uvv + F * (1 - U)) * dt
  V <- V + (Dv * Lv + uvv - (F + k) * V) * dt
}

df <- expand.grid(x = 1:n, y = 1:n)
df$z <- as.vector(V)

df$z <- (df$z - min(df$z)) / (max(df$z) - min(df$z))

colors <- c(
  "#020402",
  "#0a1a0a",
  "#1b3a1b",
  "#356b35",
  "#6bbf6b",
  "#d9ffd9"
)

P6 <- ggplot(df, aes(x, y, fill = z)) +
  geom_raster() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
P6

library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)
library(gifski)
library(tuneR)
library(png)

Kubrick Palette - Static Stills

The final Kubrick-esque colour system applied to the gravitational wave spiral. Palette: Void → HAL’s Iris → Corridor Fog → Asylum White → Period Amber → Overlook Blood.

# Grid
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r     <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Asymmetric distortion - breaks perfect rotational symmetry like a real merger
set.seed(42)
theta_distort <- theta + 0.15 * sin(3 * theta) + runif(length(theta), -0.1, 0.1)

# Chirp effect - frequency increases toward centre (bodies accelerating inward)
wave1 <- sin((8  + 2.0 / (r + 0.1)) * r + 5 * theta_distort) * exp(-0.5 * r)
wave2 <- 0.6 * sin((12 + 1.5 / (r + 0.1)) * r - 3 * theta_distort) * exp(-0.8 * r)
wave3 <- 0.4 * sin((16 + 1.0 / (r + 0.1)) * r + 7 * theta_distort) * exp(-1.0 * r)

# Fractal noise weighted by radial decay
noise_base <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 4,
  x = grid$x,
  y = grid$y,
  frequency = 1.2
)
noise <- noise_base * exp(-0.5 * r)

# Tighter singularity + accretion disc ring
sphere <- exp(-r^2 * 8) * 5 + 0.3 * exp(-((r - 0.5)^2) * 20)

# Radial vignette - outer edge fades to void
vignette <- 1 - pmin(1, (r / 5)^2)

# Combine
grid$z <- (wave1 + wave2 + wave3 + noise + sphere) * vignette

# Kubrick palette
colors_kubrick <- c("#0A0A0F", "#0D0D14", "#1C2B35", "#4A7C8E",
                    "#C8BFA8", "#F5F0E8", "#C8922A", "#8B1A1A")
kubrick_still <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(
    colors = colors_kubrick,
    limits = c(-1.5, 4.5),
    oob = scales::squish,
    guide = "none"
  ) +
  theme_void() +
  theme(plot.background = element_rect(fill = "#0A0A0F", color = NA))

kubrick_still

# Save high-res still
ggsave("kubrick_still.png", kubrick_still,
       width = 8, height = 8, dpi = 300, bg = "#0A0A0F")
kubrick_still

Animated GIF - Seamless Loop

Audio-synced animation: waves propagate outward, singularity flares on beat.

# Audio sync parameters
fps             <- 24
duration_gif    <- 5        # seconds for GIF loop
n_frames_gif    <- fps * duration_gif   # 120 frames
bpm             <- 120
beat_freq       <- bpm / 60

beat_envelope <- function(i, fps, beat_freq, decay_secs = 0.10) {
  beat_phase <- (i * beat_freq / fps) %% 1
  exp(-beat_phase / (decay_secs * beat_freq))
}

hum_envelope <- function(i, fps) {
  0.06 * sin(2 * pi * 0.5 * (i / fps))
}
frame_dir <- tempdir()

for (i in seq_len(n_frames_gif)) {

  t          <- (2 * pi * (i - 1)) / n_frames_gif
  beat       <- beat_envelope(i, fps, beat_freq)
  hum        <- hum_envelope(i, fps)

  sphere_intensity <- 3 + beat * 4
  sphere_frame <- exp(-r^2 * 8) * sphere_intensity +
                  0.3 * exp(-((r - 0.5)^2) * 20)

  phase_jolt <- beat * 0.4

  wave1 <- sin((8  + 2.0 / (r + 0.1)) * r + 5 * theta_distort - 2.5 * t - phase_jolt) *
            exp(-0.5 * r)
  wave2 <- 0.6 * sin((12 + 1.5 / (r + 0.1)) * r - 3 * theta_distort - 3.5 * t -
            phase_jolt * 1.2) * exp(-0.8 * r)
  wave3 <- 0.4 * sin((16 + 1.0 / (r + 0.1)) * r + 7 * theta_distort - 4.5 * t -
            phase_jolt * 1.5) * exp(-1.0 * r)

  grid$z <- (wave1 + wave2 + wave3 + noise + sphere_frame) * vignette + hum

  p <- ggplot(grid, aes(x, y, fill = z)) +
    geom_raster(interpolate = TRUE) +
    coord_fixed() +
    scale_fill_gradientn(
      colors = colors_kubrick,
      limits = c(-2.5, 4.5),
      oob = scales::squish,
      guide = "none"
    ) +
    theme_void() +
    theme(plot.background = element_rect(fill = "#0A0A0F", color = NA))

  ggsave(
    filename = file.path(frame_dir, sprintf("frame_%03d.png", i)),
    plot = p, width = 6, height = 6, dpi = 150, bg = "#0A0A0F"
  )

  cat(sprintf("Frame %d / %d\n", i, n_frames_gif))
}

# Stitch - drop duplicate boundary frame, loop = 0 for no pause
frames <- sort(list.files(frame_dir, pattern = "frame_.*\\.png", full.names = TRUE))
frames <- frames[-length(frames)]
gifski(frames, gif_file = "kubrick_spiral.gif",
       width = 900, height = 900, delay = 1/24, loop = 0)

cat("Done: kubrick_spiral.gif\n")

MRI Audio

Synthesised beat at 120 BPM with chest-wall low-pass filter.

sample_rate <- 44100
duration    <- 30
bpm         <- 120
beat_freq   <- bpm / 60

# Lub: lower frequency, longer - closure of mitral/tricuspid valves
make_lub <- function(sr) {
  t        <- seq(0, 0.12, by = 1/sr)
  envelope <- exp(-t / 0.04) * (1 - exp(-t / 0.003))
  tone     <- sin(2 * pi * 55  * t) * 0.7 +
              sin(2 * pi * 90  * t) * 0.3 +
              sin(2 * pi * 140 * t) * 0.15
  tone * envelope
}

# Dub: higher frequency, shorter - closure of aortic/pulmonary valves
make_dub <- function(sr) {
  t        <- seq(0, 0.07, by = 1/sr)
  envelope <- exp(-t / 0.02) * (1 - exp(-t / 0.001))
  tone     <- sin(2 * pi * 75  * t) * 0.5 +
              sin(2 * pi * 120 * t) * 0.4 +
              sin(2 * pi * 200 * t) * 0.2
  tone * envelope
}

# Chest-wall filter - muffles everything above ~250 Hz
apply_stethoscope_filter <- function(signal, sr) {
  smoothed <- signal
  for (k in 1:8) {
    smoothed <- stats::filter(smoothed, rep(1/5, 5), sides = 2)
  }
  smoothed[is.na(smoothed)] <- 0
  breath <- rnorm(length(signal), 0, 0.008)
  as.numeric(smoothed) + breath
}

lub <- make_lub(sample_rate)
dub <- make_dub(sample_rate)

total_samples <- duration * sample_rate
signal        <- numeric(total_samples)
beat_interval <- round(sample_rate / beat_freq)

for (i in seq(0, duration * beat_freq - 1)) {
  beat_start <- round(i * beat_interval) + 1

  # Lub at beat onset
  e1  <- min(beat_start + length(lub) - 1, total_samples)
  len <- e1 - beat_start + 1
  signal[beat_start:e1] <- signal[beat_start:e1] + lub[1:len] * 0.9

  # Dub 120ms after lub
  dub_start <- beat_start + round(0.12 * sample_rate)
  if (dub_start <= total_samples) {
    e2   <- min(dub_start + length(dub) - 1, total_samples)
    len2 <- e2 - dub_start + 1
    signal[dub_start:e2] <- signal[dub_start:e2] + dub[1:len2] * 0.6
  }
}

signal <- apply_stethoscope_filter(signal, sample_rate)
signal <- signal / max(abs(signal), na.rm = TRUE)

# Stereo - fixes AAC mono encoder error
signal_int <- as.integer(signal * 32767)
mri_wave   <- Wave(left  = signal_int,
                   right = signal_int,
                   samp.rate = sample_rate, bit = 16)
writeWave(mri_wave, "mri_soundscape.wav")
cat("Audio saved: mri_soundscape.wav\n")

Full 30-Second MP4 - Pipe Direct to ffmpeg

Renders 720 frames piped directly to ffmpeg with audio muxed in. No bulk PNG storage — only one frame on disk at a time.

# Requires ffmpeg installed:
# conda install -c conda-forge ffmpeg -y

Sys.setenv(PATH = paste("/opt/conda/bin", Sys.getenv("PATH"), sep = ":"))

fps      <- 24
duration <- 30
n_frames <- fps * duration   # 720 frames
beat_freq <- 2.0

lub_envelope <- function(i, fps, beat_freq, decay_secs = 0.08) {
  beat_phase <- (i * beat_freq / fps) %% 1
  exp(-beat_phase / (decay_secs * beat_freq))
}

dub_envelope <- function(i, fps, beat_freq, decay_secs = 0.04) {
  dub_offset <- 0.12 * beat_freq
  beat_phase  <- (i * beat_freq / fps) %% 1
  dub_phase   <- (beat_phase - dub_offset) %% 1
  ifelse(dub_phase < 0.15,
         exp(-dub_phase / (decay_secs * beat_freq)),
         0)
}

breath_envelope <- function(i, fps) {
  0.04 * sin(2 * pi * 0.25 * (i / fps))
}

ffmpeg_cmd <- paste(
  "ffmpeg -y",
  "-framerate 24",
  "-f rawvideo -pixel_format rgb24",
  "-video_size 900x900",
  "-i pipe:0",
  "-i mri_soundscape.wav",
  "-c:v libx264 -pix_fmt yuv420p",
  "-c:a aac -ac 2",
  "-shortest",
  "kubrick_mri_final.mp4"
)

pipe_con <- pipe(ffmpeg_cmd, open = "wb")
tmp_png  <- tempfile(fileext = ".png")

for (i in seq_len(n_frames)) {

  t      <- (2 * pi * (i - 1)) / n_frames
  lub    <- lub_envelope(i, fps, beat_freq)
  dub    <- dub_envelope(i, fps, beat_freq)
  breath <- breath_envelope(i, fps)

  sphere_intensity <- 3 + lub * 4 + dub * 1.5
  sphere_frame <- exp(-r^2 * 8) * sphere_intensity +
                  0.3 * exp(-((r - 0.5)^2) * 20)

  phase_jolt <- lub * 0.4 + dub * 0.15

  wave1 <- sin((8  + 2.0 / (r + 0.1)) * r + 5 * theta_distort - 2.5 * t - phase_jolt) *
            exp(-0.5 * r)
  wave2 <- 0.6 * sin((12 + 1.5 / (r + 0.1)) * r - 3 * theta_distort - 3.5 * t -
            phase_jolt * 1.2) * exp(-0.8 * r)
  wave3 <- 0.4 * sin((16 + 1.0 / (r + 0.1)) * r + 7 * theta_distort - 4.5 * t -
            phase_jolt * 1.5) * exp(-1.0 * r)

  grid$z <- (wave1 + wave2 + wave3 + noise + sphere_frame) * vignette + breath

  p <- ggplot(grid, aes(x, y, fill = z)) +
    geom_raster(interpolate = TRUE) +
    coord_fixed() +
    scale_fill_gradientn(
      colors = colors_kubrick,
      limits = c(-2.5, 4.5),
      oob = scales::squish,
      guide = "none"
    ) +
    theme_void() +
    theme(plot.background = element_rect(fill = "#0A0A0F", color = NA))

  ggsave(tmp_png, p, width = 900/96, height = 900/96, dpi = 96, bg = "#0A0A0F")

  img     <- png::readPNG(tmp_png)
  raw_rgb <- as.raw(as.integer(img[,,1:3] * 255))
  writeBin(raw_rgb, pipe_con)
  file.remove(tmp_png)

  if (i %% 24 == 0) cat(sprintf("  %d / %d sec rendered\n", i %/% 24, duration))
}

close(pipe_con)
cat("Done: kubrick_mri_final.mp4\n")
knitr::opts_chunk$set(echo = TRUE)
install.packages(c(
  "tidyverse",
  "ggplot2",
  "ggforce",
  "ambient",
  "colorspace",
  "gifski",
  "tuneR",
  "png",
  "scales"
))
install.packages(c("ggplot2", "ggforce", "ambient", "colorspace"))
install.packages(c("ggplot2", "ggforce", "ambient", "colorspace"))
install.packages("deldir")
# Packages
library(tidyverse)


library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)


# Grid for plotting
n <- 800
x <- seq(-4, 4, length.out = n)
y <- seq(-4, 4, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

#  Base void (isolation)
void <- -exp(-r^2)

#  Grief waves (damped oscillations)
waves <- sin(6 * r + theta) * exp(-0.8 * r)

#  Slight asymmetry (loss)
asym <- 0.3 * sin(3 * theta) * exp(-r)

#  Textural noise (ambient)
noise <- fracture(
  noise = gen_perlin, 
  fractal = fbm,
  octaves = 4,
  x = grid$x,
  y = grid$y,
  frequency = 1
) * 0.5

# Combining everything
grid$z <- void + waves + asym + noise

# Color palette for emotional depth
colors <- diverging_hcl(100, palette = "Blue-Red 3")

# Plot
Kubrickplot7 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))

Kubrickplot7
n <- 1000  # finer resolution
x <- seq(-4, 4, length.out = n)
y <- seq(-4, 4, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Central void (isolation)
void <- -1.2 * exp(-r^2)

# Multi-layered grief waves
wave1 <- sin(6 * r + theta) * exp(-0.8 * r)
wave2 <- 0.6 * sin(10 * r - 2*theta) * exp(-1.2 * r)
wave3 <- 0.4 * sin(14 * r + 4*theta) * exp(-1.6 * r)

#  Stronger asymmetry (loss)
asym <- 0.5 * sin(4 * theta) * exp(-0.6 * r)

#  Finer fractal noise (chaotic emotion)
noise <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 6,
  x = grid$x,
  y = grid$y,
  frequency = 1.5
) * 0.4

# Combine layers
grid$z <- void + wave1 + wave2 + wave3 + asym + noise

# Color palette: dark blues, purples, and subtle red hints
colors <- diverging_hcl(150, palette = "Blue-Red 3")
colors <- rev(colors)  # make center darker

# Base plot
p <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))

# Optional: ghost Bezier curves for fading memories

spirals <- data.frame(
  x = c(0, 0.5, 1, 1.5, 2),
  y = c(0, 0.3, 0.6, 0.9, 1),
  xend = c(0.2, 0.7, 1.2, 1.7, 2.2),
  yend = c(0.2, 0.5, 0.8, 1.1, 1.3)
)

p + geom_bezier(aes(x = x, y = y, xend = xend, yend = yend),
                data = spirals, color = "#ffffff20", size = 0.3, inherit.aes = FALSE)
p
library(ggplot2)
library(ggforce)


# Grid
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Stronger angular distortion (clearly visible)
set.seed(42)
theta_distort <- theta + 0.25 * sin(3 * r)

# More dynamic spiral layering
wave1 <- sin(6 * r + 4 * theta_distort) * exp(-0.4 * r)
wave2 <- 0.8 * sin(10 * r - 5 * theta_distort) * exp(-0.7 * r)
wave3 <- 0.6 * sin(18 * r + 8 * theta_distort) * exp(-1.1 * r)

# Noise weighted toward center
noise <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 3,
  x = grid$x,
  y = grid$y,
  frequency = 1
)

noise <- noise * exp(-0.3 * r)

# Much stronger central sphere
sphere <- 2.5 * exp(-r^2 * 2)

# Radial darkening toward edges (adds depth)
vignette <- -0.8 * (r / max(r))^2

# Combine
z <- wave1 + wave2 + wave3 + noise + sphere + vignette

# Normalize values for stronger contrast
z <- (z - min(z)) / (max(z) - min(z))

grid$z <- z

# More dramatic gradient
colors <- c(
  "#0f1c2e",   # deep blue
  "#2e4057",   # muted blue
  "#5d6d7e",   # cool grey-blue
  "#9b59b6",   # muted purple
  "#f39c12"    # orange core
)

Plot2 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))

library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)

# Memory-efficient grid
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Add subtle random perturbation to angles for organic feel
set.seed(42)
theta_distort <- theta + runif(length(theta), -0.1, 0.1)

# Spiral layers
wave1 <- sin(8 * r + 5 * theta_distort) * exp(-0.5 * r)
wave2 <- 0.6 * sin(12 * r - 3 * theta_distort) * exp(-0.8 * r)
wave3 <- 0.4 * sin(16 * r + 7 * theta_distort) * exp(-1.0 * r)

# Fractal noise stronger near the center
noise_base <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 4,
  x = grid$x,
  y = grid$y,
  frequency = 1.2
)
# weight noise by radial decay
noise <- noise_base * exp(-0.5 * r)

# Central orange sphere (larger and softer glow)
sphere <- exp(-r^2 * 3)  # adjust multiplier for glow size

# Combine layers numerically
grid$z <- wave1 + wave2 + wave3 + noise + sphere

# Gradient palette (blues/purples fading into orange near center)
colors <- c("#2e4057", "#5d6d7e", "#7f8c8d", "#f39c12")

# Plot
p2 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
p2
Plot2

library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)

# Grid (memory-safe)
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Angular distortion (structured, not random)
theta_distort <- theta + 0.3 * sin(2 * r)

# Spiral equation
spiral_field <- 10 * r + 6 * theta_distort

# Create thin ridge lines using cosine sharpening
ridge <- cos(spiral_field)

# Sharpen ridges (higher power = thinner lines)
ridge_sharp <- abs(ridge)^20

# Fade outward slightly
ridge_sharp <- ridge_sharp * exp(-0.25 * r)

# Add faint secondary counter-rotation for tension
counter <- abs(cos(8 * r - 5 * theta))^18
counter <- 0.4 * counter * exp(-0.4 * r)

# Strong central sphere
sphere <- 2.5 * exp(-r^2 * 2)

# Darken edges for depth
vignette <- -0.6 * (r / max(r))^2

# Combine
z <- ridge_sharp + counter + sphere + vignette

# Normalize (important for contrast)
z <- (z - min(z)) / (max(z) - min(z))

grid$z <- z

# Controlled palette (ink-like)
colors <- c(
  "#0b132b",  # deep navy
  "#1c2541",  # muted blue
  "#3a506b",  # steel blue
  "#5d6d7e",  # cool grey-blue
  "#9b59b6",  # muted purple
  "#f39c12"   # orange core
)

P3 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
P3

library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)

# Grid (memory-safe)
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

r <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Angular distortion (structured, not random)
theta_distort <- theta + 0.3 * sin(2 * r)

# Spiral equation
spiral_field <- 10 * r + 6 * theta_distort

# Create thin ridge lines using cosine sharpening
ridge <- cos(spiral_field)

# Sharpen ridges (higher power = thinner lines)
ridge_sharp <- abs(ridge)^15

# Fade outward slightly
ridge_sharp <- ridge_sharp * exp(-0.25 * r)

# Add faint secondary counter-rotation for tension
counter <- abs(cos(8 * r - 5 * theta))^18
counter <- 0.4 * counter * exp(-0.4 * r)

# Strong central sphere
sphere <- 2.5 * exp(-r^2 * 2)

# Darken edges for depth
vignette <- -0.6 * (r / max(r))^2

# Combine
z <- ridge_sharp + counter + sphere + vignette

# Normalize (important for contrast)
z <- (z - min(z)) / (max(z) - min(z))

grid$z <- z

# Controlled palette (ink-like)
colors <- c(
  "#0b132b",  # deep navy
  "#1c2541",  # muted blue
  "#3a506b",  # steel blue
  "#5d6d7e",  # cool grey-blue
  "#f39c12"   # orange core
)

P4 <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
P4
library(ggplot2)
library(deldir)

set.seed(42)

# Number of colony seeds (increase for denser mould)
n_points <- 80

# Random seed locations
x <- runif(n_points, -5, 5)
y <- runif(n_points, -5, 5)

# Compute Voronoi tessellation
vor <- deldir(x, y)
tiles <- tile.list(vor)

# Convert to dataframe for ggplot
mould_df <- do.call(rbind, lapply(seq_along(tiles), function(i) {
  data.frame(
    x = tiles[[i]]$x,
    y = tiles[[i]]$y,
    id = i
  )
}))

# Add organic distortion
mould_df$x <- mould_df$x + 0.2 * sin(3 * mould_df$y)
mould_df$y <- mould_df$y + 0.2 * cos(3 * mould_df$x)

mould_df$x <- mould_df$x + runif(nrow(mould_df), -0.1, 0.1)
mould_df$y <- mould_df$y + runif(nrow(mould_df), -0.1, 0.1)

n_points <- 150
# Plot
Mould1 <- ggplot(mould_df, aes(x, y, group = id)) +
  geom_polygon(
    fill = "#4b5d3a", 
    color = "#1f2a1f",
    alpha = 0.85
  ) +
  coord_fixed(xlim = c(-5,5), ylim = c(-5,5)) +
  theme_void() +
  theme(
    plot.background = element_rect(fill = "#0d140d", color = NA)
  )
library(ggplot2)

n <- 250

Du <- 0.16
Dv <- 0.08
F  <- 0.060
k  <- 0.062
dt <- 1

# Initialize
U <- matrix(1, n, n)
V <- matrix(0, n, n)

# Add random noise across grid (break symmetry)
set.seed(1)
V <- V + matrix(runif(n*n, 0, 0.05), n, n)

# Seed larger irregular region in center
center <- (n/2 - 20):(n/2 + 20)
U[center, center] <- 0.50
V[center, center] <- 0.25

laplace <- function(M) {
  -4*M +
    rbind(M[-1,], M[n,]) +
    rbind(M[1,], M[-n,]) +
    cbind(M[,-1], M[,n]) +
    cbind(M[,1], M[,-n])
}

steps <- 2000  # important increase

for (i in 1:steps) {
  Lu <- laplace(U)
  Lv <- laplace(V)
  
  uvv <- U * V * V
  
  U <- U + (Du * Lu - uvv + F * (1 - U)) * dt
  V <- V + (Dv * Lv + uvv - (F + k) * V) * dt
}

df <- expand.grid(x = 1:n, y = 1:n)
df$z <- as.vector(V)

df$z <- (df$z - min(df$z)) / (max(df$z) - min(df$z))

colors <- c(
  "#050a05",
  "#0f1f0f",
  "#1f3b1f",
  "#3f6f3f",
  "#7fbf7f",
  "#e6ffe6"
)

ggplot(df, aes(x, y, fill = z)) +
  geom_raster() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
library(ggplot2)

n <- 250

Du <- 0.16
Dv <- 0.08
F  <- 0.048
k  <- 0.060
dt <- 1

U <- matrix(1, n, n)
V <- matrix(0, n, n)

set.seed(2)

# Stronger noise everywhere
V <- V + matrix(runif(n*n, 0, 0.08), n, n)

# Seed infection in top-left corner (not center)
seed <- 1:40
U[seed, seed] <- 0.50
V[seed, seed] <- 0.25
U[200:230, 150:180] <- 0.50
V[200:230, 150:180] <- 0.25

laplace <- function(M) {
  -4*M +
    rbind(M[-1,], M[n,]) +
    rbind(M[1,], M[-n,]) +
    cbind(M[,-1], M[,n]) +
    cbind(M[,1], M[,-n])
}

steps <- 3500

for (i in 1:steps) {
  Lu <- laplace(U)
  Lv <- laplace(V)
  
  uvv <- U * V * V
  
  U <- U + (Du * Lu - uvv + F * (1 - U)) * dt
  V <- V + (Dv * Lv + uvv - (F + k) * V) * dt
}

df <- expand.grid(x = 1:n, y = 1:n)
df$z <- as.vector(V)

df$z <- (df$z - min(df$z)) / (max(df$z) - min(df$z))

colors <- c(
  "#020402",
  "#0a1a0a",
  "#1b3a1b",
  "#356b35",
  "#6bbf6b",
  "#d9ffd9"
)

P6 <- ggplot(df, aes(x, y, fill = z)) +
  geom_raster() +
  scale_fill_gradientn(colors = colors, guide = "none") +
  coord_fixed() +
  theme_void() +
  theme(plot.background = element_rect(fill = "black", color = NA))
P6
library(tidyverse)
library(ggplot2)
library(ggforce)
library(ambient)
library(colorspace)
library(gifski)
library(tuneR)
library(png)
# Grid
n <- 600
x <- seq(-5, 5, length.out = n)
y <- seq(-5, 5, length.out = n)
grid <- expand.grid(x = x, y = y)

# Polar coordinates
r     <- sqrt(grid$x^2 + grid$y^2)
theta <- atan2(grid$y, grid$x)

# Asymmetric distortion - breaks perfect rotational symmetry like a real merger
set.seed(42)
theta_distort <- theta + 0.15 * sin(3 * theta) + runif(length(theta), -0.1, 0.1)

# Chirp effect - frequency increases toward centre (bodies accelerating inward)
wave1 <- sin((8  + 2.0 / (r + 0.1)) * r + 5 * theta_distort) * exp(-0.5 * r)
wave2 <- 0.6 * sin((12 + 1.5 / (r + 0.1)) * r - 3 * theta_distort) * exp(-0.8 * r)
wave3 <- 0.4 * sin((16 + 1.0 / (r + 0.1)) * r + 7 * theta_distort) * exp(-1.0 * r)

# Fractal noise weighted by radial decay
noise_base <- fracture(
  noise = gen_perlin,
  fractal = fbm,
  octaves = 4,
  x = grid$x,
  y = grid$y,
  frequency = 1.2
)
noise <- noise_base * exp(-0.5 * r)

# Tighter singularity + accretion disc ring
sphere <- exp(-r^2 * 8) * 5 + 0.3 * exp(-((r - 0.5)^2) * 20)

# Radial vignette - outer edge fades to void
vignette <- 1 - pmin(1, (r / 5)^2)

# Combine
grid$z <- (wave1 + wave2 + wave3 + noise + sphere) * vignette

# Kubrick palette
colors_kubrick <- c("#0A0A0F", "#0D0D14", "#1C2B35", "#4A7C8E",
                    "#C8BFA8", "#F5F0E8", "#C8922A", "#8B1A1A")
kubrick_still <- ggplot(grid, aes(x, y, fill = z)) +
  geom_raster(interpolate = TRUE) +
  coord_fixed() +
  scale_fill_gradientn(
    colors = colors_kubrick,
    limits = c(-1.5, 4.5),
    oob = scales::squish,
    guide = "none"
  ) +
  theme_void() +
  theme(plot.background = element_rect(fill = "#0A0A0F", color = NA))

kubrick_still

# Save high-res still
ggsave("kubrick_still.png", kubrick_still,
       width = 8, height = 8, dpi = 300, bg = "#0A0A0F")
kubrick_still
# Audio sync parameters
fps             <- 24
duration_gif    <- 5        # seconds for GIF loop
n_frames_gif    <- fps * duration_gif   # 120 frames
bpm             <- 120
beat_freq       <- bpm / 60

beat_envelope <- function(i, fps, beat_freq, decay_secs = 0.10) {
  beat_phase <- (i * beat_freq / fps) %% 1
  exp(-beat_phase / (decay_secs * beat_freq))
}

hum_envelope <- function(i, fps) {
  0.06 * sin(2 * pi * 0.5 * (i / fps))
}
frame_dir <- tempdir()

for (i in seq_len(n_frames_gif)) {

  t          <- (2 * pi * (i - 1)) / n_frames_gif
  beat       <- beat_envelope(i, fps, beat_freq)
  hum        <- hum_envelope(i, fps)

  sphere_intensity <- 3 + beat * 4
  sphere_frame <- exp(-r^2 * 8) * sphere_intensity +
                  0.3 * exp(-((r - 0.5)^2) * 20)

  phase_jolt <- beat * 0.4

  wave1 <- sin((8  + 2.0 / (r + 0.1)) * r + 5 * theta_distort - 2.5 * t - phase_jolt) *
            exp(-0.5 * r)
  wave2 <- 0.6 * sin((12 + 1.5 / (r + 0.1)) * r - 3 * theta_distort - 3.5 * t -
            phase_jolt * 1.2) * exp(-0.8 * r)
  wave3 <- 0.4 * sin((16 + 1.0 / (r + 0.1)) * r + 7 * theta_distort - 4.5 * t -
            phase_jolt * 1.5) * exp(-1.0 * r)

  grid$z <- (wave1 + wave2 + wave3 + noise + sphere_frame) * vignette + hum

  p <- ggplot(grid, aes(x, y, fill = z)) +
    geom_raster(interpolate = TRUE) +
    coord_fixed() +
    scale_fill_gradientn(
      colors = colors_kubrick,
      limits = c(-2.5, 4.5),
      oob = scales::squish,
      guide = "none"
    ) +
    theme_void() +
    theme(plot.background = element_rect(fill = "#0A0A0F", color = NA))

  ggsave(
    filename = file.path(frame_dir, sprintf("frame_%03d.png", i)),
    plot = p, width = 6, height = 6, dpi = 150, bg = "#0A0A0F"
  )

  cat(sprintf("Frame %d / %d\n", i, n_frames_gif))
}

# Stitch - drop duplicate boundary frame, loop = 0 for no pause
frames <- sort(list.files(frame_dir, pattern = "frame_.*\\.png", full.names = TRUE))
frames <- frames[-length(frames)]
gifski(frames, gif_file = "kubrick_spiral.gif",
       width = 900, height = 900, delay = 1/24, loop = 0)

cat("Done: kubrick_spiral.gif\n")
sample_rate <- 44100
duration    <- 30
bpm         <- 120
beat_freq   <- bpm / 60

# Lub: lower frequency, longer - closure of mitral/tricuspid valves
make_lub <- function(sr) {
  t        <- seq(0, 0.12, by = 1/sr)
  envelope <- exp(-t / 0.04) * (1 - exp(-t / 0.003))
  tone     <- sin(2 * pi * 55  * t) * 0.7 +
              sin(2 * pi * 90  * t) * 0.3 +
              sin(2 * pi * 140 * t) * 0.15
  tone * envelope
}

# Dub: higher frequency, shorter - closure of aortic/pulmonary valves
make_dub <- function(sr) {
  t        <- seq(0, 0.07, by = 1/sr)
  envelope <- exp(-t / 0.02) * (1 - exp(-t / 0.001))
  tone     <- sin(2 * pi * 75  * t) * 0.5 +
              sin(2 * pi * 120 * t) * 0.4 +
              sin(2 * pi * 200 * t) * 0.2
  tone * envelope
}

# Chest-wall filter - muffles everything above ~250 Hz
apply_stethoscope_filter <- function(signal, sr) {
  smoothed <- signal
  for (k in 1:8) {
    smoothed <- stats::filter(smoothed, rep(1/5, 5), sides = 2)
  }
  smoothed[is.na(smoothed)] <- 0
  breath <- rnorm(length(signal), 0, 0.008)
  as.numeric(smoothed) + breath
}

lub <- make_lub(sample_rate)
dub <- make_dub(sample_rate)

total_samples <- duration * sample_rate
signal        <- numeric(total_samples)
beat_interval <- round(sample_rate / beat_freq)

for (i in seq(0, duration * beat_freq - 1)) {
  beat_start <- round(i * beat_interval) + 1

  # Lub at beat onset
  e1  <- min(beat_start + length(lub) - 1, total_samples)
  len <- e1 - beat_start + 1
  signal[beat_start:e1] <- signal[beat_start:e1] + lub[1:len] * 0.9

  # Dub 120ms after lub
  dub_start <- beat_start + round(0.12 * sample_rate)
  if (dub_start <= total_samples) {
    e2   <- min(dub_start + length(dub) - 1, total_samples)
    len2 <- e2 - dub_start + 1
    signal[dub_start:e2] <- signal[dub_start:e2] + dub[1:len2] * 0.6
  }
}

signal <- apply_stethoscope_filter(signal, sample_rate)
signal <- signal / max(abs(signal), na.rm = TRUE)

# Stereo - fixes AAC mono encoder error
signal_int <- as.integer(signal * 32767)
mri_wave   <- Wave(left  = signal_int,
                   right = signal_int,
                   samp.rate = sample_rate, bit = 16)
writeWave(mri_wave, "mri_soundscape.wav")
cat("Audio saved: mri_soundscape.wav\n")
# Requires ffmpeg installed:
# conda install -c conda-forge ffmpeg -y

Sys.setenv(PATH = paste("/opt/conda/bin", Sys.getenv("PATH"), sep = ":"))

fps      <- 24
duration <- 30
n_frames <- fps * duration   # 720 frames
beat_freq <- 2.0

lub_envelope <- function(i, fps, beat_freq, decay_secs = 0.08) {
  beat_phase <- (i * beat_freq / fps) %% 1
  exp(-beat_phase / (decay_secs * beat_freq))
}

dub_envelope <- function(i, fps, beat_freq, decay_secs = 0.04) {
  dub_offset <- 0.12 * beat_freq
  beat_phase  <- (i * beat_freq / fps) %% 1
  dub_phase   <- (beat_phase - dub_offset) %% 1
  ifelse(dub_phase < 0.15,
         exp(-dub_phase / (decay_secs * beat_freq)),
         0)
}

breath_envelope <- function(i, fps) {
  0.04 * sin(2 * pi * 0.25 * (i / fps))
}

ffmpeg_cmd <- paste(
  "ffmpeg -y",
  "-framerate 24",
  "-f rawvideo -pixel_format rgb24",
  "-video_size 900x900",
  "-i pipe:0",
  "-i mri_soundscape.wav",
  "-c:v libx264 -pix_fmt yuv420p",
  "-c:a aac -ac 2",
  "-shortest",
  "kubrick_mri_final.mp4"
)

pipe_con <- pipe(ffmpeg_cmd, open = "wb")
tmp_png  <- tempfile(fileext = ".png")

for (i in seq_len(n_frames)) {

  t      <- (2 * pi * (i - 1)) / n_frames
  lub    <- lub_envelope(i, fps, beat_freq)
  dub    <- dub_envelope(i, fps, beat_freq)
  breath <- breath_envelope(i, fps)

  sphere_intensity <- 3 + lub * 4 + dub * 1.5
  sphere_frame <- exp(-r^2 * 8) * sphere_intensity +
                  0.3 * exp(-((r - 0.5)^2) * 20)

  phase_jolt <- lub * 0.4 + dub * 0.15

  wave1 <- sin((8  + 2.0 / (r + 0.1)) * r + 5 * theta_distort - 2.5 * t - phase_jolt) *
            exp(-0.5 * r)
  wave2 <- 0.6 * sin((12 + 1.5 / (r + 0.1)) * r - 3 * theta_distort - 3.5 * t -
            phase_jolt * 1.2) * exp(-0.8 * r)
  wave3 <- 0.4 * sin((16 + 1.0 / (r + 0.1)) * r + 7 * theta_distort - 4.5 * t -
            phase_jolt * 1.5) * exp(-1.0 * r)

  grid$z <- (wave1 + wave2 + wave3 + noise + sphere_frame) * vignette + breath

  p <- ggplot(grid, aes(x, y, fill = z)) +
    geom_raster(interpolate = TRUE) +
    coord_fixed() +
    scale_fill_gradientn(
      colors = colors_kubrick,
      limits = c(-2.5, 4.5),
      oob = scales::squish,
      guide = "none"
    ) +
    theme_void() +
    theme(plot.background = element_rect(fill = "#0A0A0F", color = NA))

  ggsave(tmp_png, p, width = 900/96, height = 900/96, dpi = 96, bg = "#0A0A0F")

  img     <- png::readPNG(tmp_png)
  raw_rgb <- as.raw(as.integer(img[,,1:3] * 255))
  writeBin(raw_rgb, pipe_con)
  file.remove(tmp_png)

  if (i %% 24 == 0) cat(sprintf("  %d / %d sec rendered\n", i %/% 24, duration))
}

close(pipe_con)
cat("Done: kubrick_mri_final.mp4\n")