[R] Fractals in R and having fun! (and more persp and color)

[ucgamdo@ucl.ac.uk]

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!