library(viridis)
library(scico)
options(repr.plot.width=10, repr.plot.height=10)

Loading required package: viridisLite

Mandelbrot Set

mandelbrot_vectorized <- function(px = 0, py = 0, d = 2, size=1500, num_iter=100) {
    # variables
    x <- seq(px - d, px + d, length.out=size)
    y <- seq(py - d, py + d, length.out=size)
    c <- outer(x,y*1i,FUN="+")
    z <- matrix(0.0, nrow=length(x), ncol=length(y))
    k <- matrix(0.0, nrow=length(x), ncol=length(y))
    m <- matrix(0.0, nrow=length(x), ncol=length(y))
    
    for (rep in 1:num_iter) {
        index <- which(Mod(z) < 2)
        z[index] <- z[index]^2 + c[index]
        k[index] <- k[index] + exp(-abs(z[index]))
    }
    
    k[index] = num_iter

    return(k)
}

pal2 = c('#421e0f','#19071a','#09012f','#040449','#000764','#0c2c8a','#1852b1','#397dd1',
     '#86b5e5','#d3ecf8','#f1e9bf','#f8c95f','#ffaa00','#cc8000','#995700','#6a3403','#421e0f')

clr = colorRampPalette(pal2)(300)
clr[300] = 'black'

display <- function(z) {
    z[z==max(z)] = NA
    z = z %% (length(pal2)*2)
    z[is.na(z)] = length(pal2)*2 + 1
    par(mar=c(0,0,0,0))
    image(z,useRaster=TRUE, col=clr, axes=FALSE)
}

x = mandelbrot_vectorized(num_iter = 1000)
display(x)

x <- mandelbrot_vectorized(n = 1000, px=-0.759, py=0.125, d=0.03)
display(x)

x <- mandelbrot_vectorized(n = 1000, px=-0.743030, py=0.126433, d=0.009)
display(x)

x <- mandelbrot_vectorized(num_iter = 1000, px=-0.7435669, py=0.1314023, d=0.001)
display(x)

x <- mandelbrot_vectorized(num_iter = 1000, px=-0.7436499, py=0.13188204, d=0.0003)
display(x)

x <- mandelbrot_vectorized(num_iter = 1000, px=-0.74364085, py=0.13182733, d=0.00005)
display(x)

x <- mandelbrot_vectorized(num_iter = 1000, px=-0.743643135, py=0.131825963, d=0.000007)
display(x)

x <- mandelbrot_vectorized(num_iter = 1000, px=-0.743644786, py=0.1318252536, d=0.0000014)
display(x)