st: Mandelbrot sets

[Scott Merryman <scott.merryman@gmail.com>]

In the recent issue of The Stata News
(http://www.stata.com/stata-news/statanews.27.4.pdf) there is an
article by Hua Peng on generating and graphing Mandelbrot sets.

There does appear to be an error in the second block of code.  The
line res= J(201*201, 3, .) ,  I believe should be mandelbrot_set =
J(201*201, 3, .)

Warning - it can take a while for the graphs to be produced.

Scott


clear*
mata:
real scalar escape(complex scalar Z, C,
                   real scalar    R, maxiter)
{
real scalar i
for (i=0; i<maxiter; i++) {
   if (norm(Z) <= R) {
    Z = Z*Z+C
   }
else  return(i)
}
return(i)
}
end

mata:
mandelbrot_set = J(201*201, 3, .)
cnt = 1
C   = -0.8+0.156i
for(i=0; i<=200; i++) {
for(j=0; j<=200; j++) {
Z  = C(-2+i*4/200, -2+j*4/200)
n = escape(Z, C, 100, 400)
mandelbrot_set[cnt, .] = (i, j, n)
cnt++
}
}
end

getmata (x y escape)= mandelbrot_set
qui summ escape
local min = r(min)-1
local max = r(max)+1
tw contour escape y x, ccuts(`min'(1)`max') clegend(off)
graphr(m(zero)) xscale(off) xlab(,nogrid) yscale(off) ylab(,nogrid)
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/