// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

// Reference:
//
// Jonathan Blanchard's blog
// Sat, 02/27/2010
// http://jblopen.com/node/18
// Writing external functions for Scilab
// Fractal generation example
//
// Mandelbrot set
// From Wikipedia, the free encyclopedia
// http://en.wikipedia.org/wiki/Mandelbrot_set#Computer_drawings

function plotMandelbrotNaive(xsize,ysize,nmax,xmin,xmax,ymin,ymax)
    xvect = linspace( xmin, xmax, xsize );
    yvect = linspace( ymin, ymax, ysize );
    res = zeros(xsize,ysize);
    for i = 1:xsize
        for j = 1:ysize
            x = xvect(i);
            y = yvect(j);
            x0 = x;
            y0 = y;
            k = 0;
            while( x*x + y*y < 4 & k < nmax )
                xtemp = x*x - y*y + x0;
                y = 2*x*y + y0;
                x = xtemp;
                k=k+1;
            end
            if k<nmax then
                res(i,j) = k;
            else
                res(i,j) = -1;
            end
        end
    end
    // The maximum number of colors
    CMAX = 1000;
    f=gcf();
    f.color_map = jetcolormap(CMAX);
    D = res;
    k = find(res<>-1);
    D(k) = floor(res(k)/max(res(k))*CMAX);
    k = find(res==-1);
    D(k) = CMAX;
    Matplot(D);
    f.children.isoview="on";
    f.children.axes_visible=["off" "off" "off"];
endfunction

xsize = 50;
ysize = 50;
nmax = 1000;
xmin = 0.2675;
xmax = 0.2685;
ymin = 0.591;
ymax = 0.592;

tic();
plotMandelbrotNaive(xsize,ysize,nmax,xmin,xmax,ymin,ymax);
t = toc();
mprintf("Time = %f (s)\n",t);


////////////////////////////////////////////////////
//
// Vectorized
//

function plotMandelbrotVect(xsize,ysize,nmax,xmin,xmax,ymin,ymax)
    xvect = linspace( xmin, xmax, xsize );
    yvect = linspace( ymin, ymax, ysize );
    [X,Y]=meshgrid(xvect,yvect);

    Z = zeros(xsize,ysize);
    R = -ones(xsize,ysize);
    W = zeros(xsize,ysize);
    C=X+%i*Y;
    J = 1:xsize*ysize;
    for k=0:nmax
        L = J(J>0);
        Z(L)=Z(L).^2+C(L);
        W(L) = abs(Z(L));
        M = find(W(L)>2);
        R(L(M))=k;
        J(L(M)) = 0;
    end
    R = R';
    // The maximum number of colors
    CMAX = 1000;
    f=gcf();
    f.color_map = jetcolormap(CMAX);
    D = R;
    k = find(R<>-1);
    D(k) = floor(R(k)/max(R(k))*CMAX);
    k = find(R==-1);
    D(k) = CMAX;
    Matplot(D);
    f.children.isoview="on";
    f.children.axes_visible=["off" "off" "off"];
endfunction

stacksize("max");
xsize = 200;
ysize = 200;
nmax = 1000;
xmin = 0.2675;
xmax = 0.2685;
ymin = 0.591;
ymax = 0.592;

tic();
plotMandelbrotVect(xsize,ysize,nmax,xmin,xmax,ymin,ymax);
t = toc();
mprintf("Time = %f (s)\n",t);