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)
par(opar)

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$$).