[R] Fractals in R and having fun! (and more persp and color)
ucgamdo@ucl.ac.uk
ucgamdo at ucl.ac.uk
Wed Sep 17 13:45:51 CEST 2003
Well, I started playing with fractals in R, and wrote a function to
generate de Mandelbrot set, which might
be of interest to some people
###########################################################################
# Mandelbrot set
###########################################################################
mandelbrot <- function(x = c(-3.0, 1.0), # x coordinates
y = c(-1.8, 1.8), # y coordinates
b = .05, # by 'steps'
iter = 20) # maximum number of iterations
{
m = NULL # will store the results, this is the 'image' matrix
for(i in seq(x[1], x[2], by = b)) {
r = NULL # stores part of the iteration results
for(j in seq(y[1], y[2], by = b)) {
it = iter # will hold iteration at which point (i, j) breaks off
c = c(i, j) # initial point
z = c(i, j) # i: real part; j: imaginary part
for(k in 1:iter) {
# the Mandelbrot iteration formulae: z -> z*z + c
z = c(z[1]^2 - z[2]^2, 2 * z[1]*z[2]) + c
# tests if point breaks off
if((z[1] + z[2])^2 > 4) { it = k; break }
}
r = c(r, it) # stores iteration results
}
# constructs the 'image' matrix
m = rbind(m, r)
}
# the output fractal object
fractal = list(x = seq(x[1], x[2], by = b), # x coordinates
y = seq(y[1], y[2], by = b), # y coordinates
z = m) # it matrix
}
######################################################################
# here goes how it works
######################################################################
frac <- mandelbrot()
image(frac) # perhaps not very nice!
# more resolution, beware, this might take some time to calculate
frac <- mandelbrot(b = .01)
image(frac)
# now here comes the fun, lets do a persp plot of the set!
persp(log(frac$z), border = NA, theta = 225, phi = 45, col = "gray", shade
= .4)
###END###
Well, here comes what would have been my question. How to put some nice
colors to the above perspective plot? I wanted to post this question some
weeks ago, but it was answered yesterday after all the discussion on persp
and colors. If you see my previous post, you'll the definition of
surf.colors used below:
persp(log(frac$z), border = NA, theta = 225, phi = 45, col =
surf.colors(frac$z, col = c(gray(seq(0.5, 1, len = 39)), "black")), shade =
.4)
that was what I always wanted to do, and this was the source of my
confusion when checking the code in persp demo.
But now some new question arises, which I'll leave to the interested hacker:
Any C programmer who volunteers to write and link c code to R to calculate de
Mandelbrot set faster?
Does anyone have any idea on how to 'smooth' the valleys and ridges that
rise towards the set?, i.e. to avoid the stairs like appearence of the plot?
I hope you can have some fun with this!
More information about the R-help
mailing list