// 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 plotFractal(R,cmax) // Plots the fractal associated with the data in R. // R(i,j)=k if the point (i,j) has left the Mandelbrot set at iteration k, // R(i,j)=-1 if the point is in the set. // cmax: the number of colors in the jet color map. f=scf(); 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 function R = computeMandelbrotNaive(xsize,ysize,nmax,xmin,xmax,ymin,ymax) // Compute the fractal associated with the data in R. // R(i,j)=k if the point (i,j) has left the Mandelbrot set at iteration k, // R(i,j)=-1 if the point is in the set. xvect = linspace( xmin, xmax, xsize ); yvect = linspace( ymin, ymax, ysize ); R = 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 k0); 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'; 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); PPS = floor(xsize*ysize/t); mprintf("PPS = %d\n",PPS); plotFractal(R,1000);