# Orbit trapped Julia fractal

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
function(z, alpha) {
f <- z * exp(1i*alpha)
z <-min(abs(Re(z)), abs(Im(z))) / 0.03
}
# choose the Julia point
complex(real = -0.687, imaginary = 0.299015)
juliaPoint <-
# main function
Vectorize(function(x, y, juliaPoint, alpha) {
Julia <-# counter
0L
i <-# current point, to be iterated
complex(real = x, imaginary = y)
z <-# iterations
while(i < 100L && Mod(z) < 100 && (i < 2L || f(z, alpha) > 1)) {
z^2 + juliaPoint
z <- i + 1L
i <-
}# now assign a color to the resulting z
2 * f(z, alpha)
fz <-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
2048L
n <- seq(-2, 2, length.out = n)
x_ <- seq(-1.5, 1.5, length.out = n)
y_ <- t(outer(x_, y_, Julia, juliaPoint = juliaPoint, alpha = 0))
img <-# plot
par(mar = c(0, 0, 0, 0), bg = "black")
opar <-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\)).