par(mai = c(0, 0, 0, 0), bg="black")

d <- 1:1100
theta <- 2 * pi/((1 + sqrt(5)) / 2)^2 # golden ratio
df <- data.frame(x = d * cos(theta * d), y = d * sin(theta * d))
plot(df,
     cex = d / 500,
     pch = 19,
     col = rainbow(6000),
     asp = 1,
     axes=F
     )