fenxing.zip_dressmrs_fenwei_差分盒维_毯子法_盒维数

function brown(maxr, npath) % BROWN Simulates electrochemical aggregation using Brownian motion % BROWN(MAXR, NPATH) Compute NPATH particles, use maximum radius % of injection circle MAXR. % % Examples: % brown( 50, 500) [ca. 13 seconds on 1.6 GHz Athlon] % brown(100, 2000) [ca. 3 minutes] % brown(200, 8000) [ca. 40 minutes] % See pp. 440-444, "Chaos and Fractals" by Peitgen, Juergens, Saupe % (2004) for further information on this experiment. % Author : Andreas Klimke, IANS, Universitaet Suttgart % Version: 1.0 % Date : Nov. 8, 2004 if nargin < 2 error('Type "help brown" for calling syntax.'); end % Maximum number of random iterations to perform in the inner loop % (100 random variables are computed outside of the loop since the % rand call is expensive and loops involving a call to rand are not % JIT-compiled). maxrand = 100; % Init matrix of sticky particles A = zeros(maxr*2, maxr*2); % Init 1st sticky particle at center of matrix A(maxr,maxr) = 1; % The direction the particle moves dir = 0; % The initial injection radius; this is increased up to maxr. This % improves performance, since the likelyhood of hitting the sticky % particle would be very low in the beginning if one would use the % full injection radius right from the beginning. r = 2; % Counter for total number of iterations niter = 0; % If equal to one, compute new path newval = 1; % Counter for the number of completed particle paths k = 0; while k < npath random = rand(maxrand,1); count = 2; if newval > 0 % Compute new point on injection circle phi = random(1); px = round(r*cos(2*phi*pi)); py = round(r*sin(2*phi*pi)); end newval = 0; while 1 % Arbitrarily select direction to turn to rval = random(count); if rval < 0.25 % EAST px = px + 1; elseif rval < 0.5 % NORTH py = py + 1; elseif rval < 0.75 % WEST px = px - 1; else % SOUTH py = py - 1; end actrsqr = px^2+py^2; if (actrsqr > r^2) % Check if particle has left the area of interest; if yes, % discard and start with new particle path. newval = 1; break; elseif (actrsqr < r^2) % Check if sticky particle is hit; if yes, attach particle % and do the next path. x = px + maxr; y = py + maxr; if (A(x+1,y) > 0 | A(x-1,y) > 0 | A(x,y+1) > 0 | A(x,y-1) > 0) A(x,y) = 1; newval = 1; k = k + 1; break; end end count = count + 1; if count > maxrand break; end niter = niter + 1; end % Compute new injection radius, (the maximum over twice the % radius of the last sticky particle and the current radius % should be enough) actmaxr = 2*sqrt(actrsqr); if actmaxr > r r = min(floor(actmaxr),maxr-1); end end disp(['Number of iterations: ', num2str(niter)]); imagesc(A); colormap(1-gray); axis equal; axis tight;