spiral.R

library(purrr)
library(ggplot2)
library(tibble)

# https://community.wolfram.com/groups/-/m/t/2575939
x <- function(t, m, gamma, epsilon0, q, epsilon) {
  numerator <- m * epsilon0 * exp(-gamma * t / m) * (
    (epsilon * q * exp(gamma * t / m)) - (gamma * sin((epsilon * q * t) / m)) - (epsilon * q * cos((epsilon * q * t) / m))
  )
  denominator <- gamma^2 + (epsilon^2 * q^2)
  return(numerator / denominator)
}

y <- function(t, m, gamma, epsilon0, q, epsilon) {
  numerator <- epsilon0 * exp(-gamma * t / m) * (
    (epsilon^2 * q^2 * exp(gamma * t / m)) - (gamma * m * cos((epsilon * q * t) / m)) + 
      (epsilon * m * q * sin((epsilon * q * t) / m)) + (gamma^2 * exp(gamma * t / m)) + 
      (gamma * m * exp(gamma * t / m))
  )
  denominator <- gamma^2 + (epsilon^2 * q^2)
  return(numerator / denominator)
}

# Example parameters
t <- seq(0,100, length=1e3)
m <- 1
gamma <- 0.5
epsilon0 <- 1
q <- 1
epsilon <- 1

# Calculating x(t) and y(t)
x_val <- x(t, m, gamma, epsilon0, q, epsilon)
y_val <- y(t, m, gamma, epsilon0, q, epsilon)

set.seed(123)
spiral <- function(m = 1,
                   gamma = 0.5,
                   epsilon0 = 1,
                   q = 1,
                   epsilon = 1, x0=0,y0=0){
  t0 <- runif(1, -1, 1)
  tf <- runif(1, 50, 300)
  t <- seq(t0, tf, length=1e3)
  data.frame(x = x0+x(t, m, gamma, epsilon0, q, epsilon),
             y = y0+y(t, m, gamma, epsilon0, q, epsilon))
  
}
# 
# plot(spiral(), t='l', xlim=c(-1,1), ylim=c(-1,5))
# lines(spiral(m=2), t='l')
# lines(spiral(q=2), t='l')
# lines(spiral(epsilon0=1.2), t='l')
# lines(spiral(q=-1.2), t='l')
# lines(spiral(m=-1.2), t='l')
# plot(spiral(t = seq(0,100,length=2e3),m=2.2, q=4.1), t='l', asp=1)
# plot(spiral(t = seq(0,100,length=2e3),m=2.2, q=-4.1), t='l', asp=1)
# lines(spiral(t = seq(0,100,length=2e3),m=2.2, q=2.1), t='l', asp=1)

params <- tribble(
  ~m, ~q, ~gamma, ~y0,~x0,
  2.2, 4.1, 0.5, 1.5, 0,
  2.2, -4.1,0.5, 1.5, 0,
  2.2, 2.1,0.5,1.5, 0,
  2.2, 0.1, 0.5,  3.5, 0,
  2.2, 0.1, 0.05, -1, 0,
  2.2, -2.1, 0.3, 5, 0,
  2.2, -2.1, 0.08, 12, 0,
  2.2, -2.1, 0.10, 3, -0.5,
  100.2, -0.1, 4.10, -1.5, -0.2,
  100.2, -0.1, 4.10, -1.5, -0.31,
  50.2, -0.12, 5.10, -2.12, 0.31,
  50.2, -0.12, 2.10, -1.12, 0.21,
  50.2, -6.12, 2.10, -2.12, 0.21,
  50.2, -12.12, 5.10, 1.12, 0.21,
  50.2, 23.12, 5.10, 1.12, 0.21,
  500.2, 43.12, 5.10, -2.12, 0.21,
  600.2, -4.12, 5.10, -2.12, -0.21,
  
  600.2, 4.12, 5.10, -2.12, 2.21,
  400.2, 2.12, 5.10, -2.12, 2.811,
  
  15.2, -6.12, 2.10, 6.12, 2.21,
  25.2, -12.12, 5.10, 8.12, 2.21,
  25.2, 23.12, 5.10, 9.12, 3.21,
  50.2, 43.12, 5.10, 12.12, 3.21,
  
  1.2, 0.1, 0.05, -2, 2,
  1.2, -2.1, 0.3, 3, 1.7,
  1.2, 2.1, 0.10, 6, -1,
  1.2, -2.1, 0.10, 2, -0.5
)

d <- pmap_df(params, spiral, .id='id')

p <- ggplot(d, aes(x,y , group=id)) +
  geom_path(col='white',lwd=0.2) +
  coord_equal(xlim=c(-5,5),ylim=c(0,20)) +
  theme_void() +
  theme(legend.position = 'none', plot.background = element_rect(fill='black'))
# 
# grid.rect(gp=gpar(fill='black'))
# print(p, vp=viewport(gp=gpar(fill='black')))

p
ggsave(p, file='bubble.svg',width=5,height=10)

# library(vwline)
# # x <- c(.2, .8, .8, .2)
# # y <- c(.9, .3, .9, .3)
# w <- c(0, .2)
# 
# d <- pmap_df(params, spiral, .id='id', t=seq(1,100, length=100))
# 
# for(l in split(d, d$id))
# grid.brushXspline(verticalBrush, l$x, l$y, w, shape=1, gp=gpar(col="black"))
#