### Copyright (C) 2002 Martin Maechler
###
### This is free software; you can redistribute it and/or modify it under
### the terms of the GNU General Public License as published by the Free
### Software Foundation; either version 2 of the License, or (at your
### option) any later version.
###
### The software is distributed in the hope that it will be useful, but
### WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
### or FITNESS FOR A PARTICULAR PURPOSE.
### See the GNU General Public License for more details.
###
### You should have received a copy of the GNU General Public License
### along with R; if not, you can obtain it via the World Wide Web at
### 'http://www.gnu.org/copyleft/gpl.html', or by writing to the Free
### Software Foundation, 59 Temple Place -- Suite 330, Boston, MA
### 02111-3307, USA.

mandelbrot01 <- function(xlim = c(-2, 0.5), ylim = c(-1,1),
                         nx, ny = nx * diff(ylim)/diff(xlim), iter = 120, r = 2)
{
  ## Purpose: Mandelbrot set for nx x ny  equispaced c in  xlim + i*ylim
  ## ----------------------------------------------------------------------
  ## Arguments: xlim, ylim: specify range of `c'
  ##            nx, ny    : number of points in x- and y-direction
  ##            iter      : number of iterations
  ##            r         : radius; if Mod(z) > r, we assume divergence
  ## ----------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 22 Feb 2002, 16:50
    x <- seq(xlim[1], xlim[2], length = nx)
    y <- seq(ylim[1], ylim[2], length = ny)
    cc <- outer(x, y, function(x,y) complex(real=x, imag = y))
    z <- cc # first z = 0 ==> z <- z^2 + c = c
    for(i in 1:iter)
        z <- z*z + cc # z*z is almost a factor 3 faster than z^2
    return(x = x, y = y, is.mandelbrot = is.finite(z) & Mod(z) < r)
}

system.time(m1 <- mandelbrot01(nx = 200, ny = 200))
## 7.00 {Pentium III, 700 MHz}
image(m1$x, m1$y, m1$is.man, asp = 1, col = c("blue", "yellow"))

system.time(m2 <- mandelbrot01(nx = 250))
image(m2$x, m2$y, m2$is.man, asp = 1, col = c("blue", "yellow"))


mandelbrot <- function(xlim = c(-2, 0.5), ylim = c(-1,1),
                       nx, ny = nx * diff(ylim)/diff(xlim), iter = 200, r = 2,
                       trace = FALSE)
{
    ## Purpose: Mandelbrot iteration count for nx x ny
    ##          equispaced c in  xlim + i*ylim;  #{iter; |z| > r}
    ## ----------------------------------------------------------------------
    ## Arguments: xlim, ylim: specify range of `c'
    ##            nx, ny    : number of points in x- and y-direction
    ##            iter      : number of iterations
    ##            r         : radius; if Mod(z) > r, we assume divergence
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 22 Feb 2002, 17:30
    if((nx <- as.integer(nx)) <= 0) stop("`nx' must be in 1,2,..")
    if((ny <- as.integer(ny)) <= 0) stop("`ny' must be in 1,2,..")
    if((iter <- as.integer(iter)) <= 0) stop("`iter' must be in 1,2,..")
    x <- seq(xlim[1], xlim[2], length = nx)
    if(ylim[1] > ylim[2]) stop("`ylim' is not increasing")
    Ny <- ny
    if(prod(ylim) >= 0)
        y <- seq(ylim[1], ylim[2], length = ny)
    else { ## make sure to have y = 0, and make use of symmetry around y = 0:
        symm.neg <- abs(ylim[1]) > ylim[2]
        Ny <- ceiling(ny * abs(ylim[1])/diff(ylim))
        y <- {
            if(symm.neg) seq(ylim[1], 0, length = Ny)
            else         seq(0, ylim[2], length = Ny)
        }
    }
    acc <- cc <- outer(x, y, function(x,y) complex(real=x, imag = y))
    az <- z <- cc ## first z = 0 ==> z <- z^2 + c = c
    sml <- rep(TRUE, nx * Ny)
    isml <- 1:(nx*Ny)                   # need this really?
    res <- matrix(iter+ 1:1, nx, Ny)    #integer
    ## acc and az are the `active' values of c (= cc) and z, i.e.,
    ## acc == c(cc[sml]),  az == c(z[sml])
    for(i in 1:iter) {
        asmall <- Mod(az <- az*az + acc) <= r
        ## NOT(asmall)  are those that are newly > r
        res[isml[!asmall]] <- i
        az    <- az[asmall]
        acc  <- acc[asmall]
        isml <-isml[asmall]
        sml[sml] <- asmall
        if(trace) cat(i,": no{sml}= ", length(az),"\n")
    }
    if(prod(ylim) < 0) { ## fill `symmetry' result
        if(symm.neg) {
            i <- y >= -ylim[2] ## (-y) <= ylim[2]
            y <- c(y, -rev(y[i])[-1])
            res <- res[, c(1:Ny, rev(which(i))[-1])]
        } else {
            i <- y <= -ylim[1]
            y <- c(-rev(y[i][-1]), y)
            res <- res[, c(rev(which(i))[-1], 1:Ny)]
        }
    }
    return(x = x, y = y, n.iter = res)
}

system.time(mm <- mandelbrot(nx = 250, trace =TRUE))
## 6.19 {Pentium III, 700 MHz}
image(mm$x, mm$y, mm$n.iter, asp = 1)
n.it <- mm$n.iter; i.max <- max(n.it)
image(mm$x, mm$y, n.it, asp = 1, frame.plot = FALSE,
      xlab = "Re(c)", ylab = "Im(z)", main = "Mandelbrot set + #{iter.}",
      breaks = c(quantile(n.it[n.it < i.max], probs=seq(0,1, len=16)), i.max),
      col = c(rev(heat.colors(15)), "black"))


system.time(mm2 <- mandelbrot(nx = 500, iter = 200))
## 18.47

iterImage <- function(xyi, n.col = 16,
                      col.type = c("heat", "rainbow", "terrain", "topo", "cm",
                      "gray", "zebra", "b&w"), col = NULL,
                      breaks = c(quantile(n.it[n.it < i.mx],
                                 probs=seq(0,1, len=n.col),names=FALSE), i.mx),
                      asp = 1, frame.plot = FALSE,
                      xlab = "Re(c)", ylab = "Im(z)",
                      tit = "", main = "Set & #{iter.}")
{
  ## Purpose: image() plot of mandelbrot like result
  ## ----------------------------------------------------------------------
  ## Arguments:
  ## ----------------------------------------------------------------------
  ## Author: Martin Maechler, Date: 22 Feb 2002, 19:02
    if(dim(n.it <- xyi[[3]]) !=
       c(nx <- length(x <- xyi$x), ny <- length(y <- xyi$y)))
        stop("invalid `xyi' argument")
    if(is.null(col)) {
        col.type <- match.arg(col.type)
        n1 <- n.col - 1
        col <- c(switch(col.type,
                        heat   = rev(heat.colors(n1)),
                        rainbow= {k <- n1%/%6 ; rainbow(n1)[c(k:n1,1:(k-1))]},
                        terrain= rev(terrain.colors(n1)),
                        topo   = rev(topo.colors(n1)),
                        cm     = rev(  cm.colors(n1)),
                        gray   = gray((1:(n1))/n.col),
                        "b&w"=rep(c("white","black"), len = n1),
                        zebra=rep(c("gray85","gray30"),len = n1),
                        ), "black")
    } else col.type <- "user specified"
    i.mx <- max(n.it)
    image(x, y, n.it, breaks = breaks, col = col,
          asp = asp, frame.plot = frame.plot,
          xlab = xlab, ylab = ylab, main = paste(tit, main))
    mtext(paste(nx," x ",ny,"  points  --- color type: \"", col.type,"\",  ",
                ## should do the following only when i.mx > 1
                col[length(col)]," means `> ",i.mx - 1,"iter'", sep=""))
}


iterImage(mm2, tit = "Mandelbrot")
iterImage(mm2, 32, tit = "Mandelbrot")

iterImage(mm2, 64, tit = "Mandelbrot", col.type = "topo")
iterImage(mm2, 77, tit = "Mandelbrot", col.type = "rainbow")
iterImage(mm2, 99, tit = "Mandelbrot", col.type = "terrain")
iterImage(mm2, 30, tit = "Mandelbrot", col.type = "cm")
iterImage(mm2, 30, tit = "Mandelbrot", col.type = "gray") # really beautiful
iterImage(mm2, 30, tit = "Mandelbrot", col.type = "zebra")
iterImage(mm2, 30, tit = "Mandelbrot", col.type = "b&w")# Black & White
