Skip to content

Instantly share code, notes, and snippets.

@baptiste
Created September 3, 2024 07:37
Show Gist options
  • Save baptiste/24d59b7d2a02d9411a4365801c55e3e1 to your computer and use it in GitHub Desktop.
Save baptiste/24d59b7d2a02d9411a4365801c55e3e1 to your computer and use it in GitHub Desktop.
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"))
#
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment