Orbit trapped Julia fractal

[This article was first published on Saturn Elephant, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Given a complex number \(z_c\), called the Julia point, the corresponding Julia fractal is obtained by iterating \(z \mapsto z^2 + z_c\) for each complex number \(z\) until the modulus of \(z\) exceeds a certain threshold or the maximal number of iterations is attained. Then a color is assigned to \(z\).

An orbit trapped Julia fractal is obtained in the same way, but the iteration is stopped whenever \(z\) is close enough to a given set such as a square or a circle. In the example shown below, we take the two axes as this trapping set.

I also add something: instead of looking at the distance between \(z\) and the two axes, I look at the distance between \(z\) and the axes after having rotated \(z\) by an angle \(\alpha\). Then I’ll vary \(\alpha\) to get an animation.

So here is the code of this algorithm; the color assigned to the final \(z\) is defined in function of the value of the trapping function (the distance):

# trapping function:
#   distance (up to factor 1/0.03) between alpha-rotated z and axes
f <- function(z, alpha) {
  z <- z * exp(1i*alpha)
  min(abs(Re(z)), abs(Im(z))) / 0.03

# choose the Julia point
juliaPoint <- complex(real = -0.687, imaginary = 0.299015)

# main function
Julia <- Vectorize(function(x, y, juliaPoint, alpha) {
  # counter
  i <- 0L
  # current point, to be iterated
  z <- complex(real = x, imaginary = y)
  # iterations
  while(i < 100L && Mod(z) < 100 && (i < 2L || f(z, alpha) > 1)) {
    z <- z^2 + juliaPoint
    i <- i + 1L
  # now assign a color to the resulting z
  fz <- 2 * f(z, alpha)
  hsv( # h, s, v must be in (0, 1)
    h = (Arg(z) + pi) / (2 * pi), 
    s = min(1, fz), 
    v = max(0, min(1, 2 - fz))

The condition i < 2L ensures that the iteration is not stopped at the beginning. Let’s plot a first image:

# run the orbit trapping of Julia
n <- 2048L
x_ <- seq(-2,   2,   length.out = n)
y_ <- seq(-1.5, 1.5, length.out = n)
img <- t(outer(x_, y_, Julia, juliaPoint = juliaPoint, alpha = 0))
# plot
opar <- par(mar = c(0, 0, 0, 0), bg = "black")
plot(c(-100, 100), c(-100, 100), type = "n", asp = 3/4, 
     xlab = NA, ylab = NA, axes = FALSE, xaxs = "i", yaxs = "i")
rasterImage(img, -100, -100, 100, 100)

And here is the animation obtained by varying the angle \(\alpha\):

Such a fractal is easily and efficiently rendered as a shader. Click here to play with the shader, in which the cursor of the mouse is used to take the Julia point. I also modified the trapping function (\(z^3\) instead of \(z\)).

To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)