// 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-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);