%
%  diffusion limited aggregation
%

N = 80;  L = 100;  m = 2;  St = 1000;

patt = [ones(1,N) ; zeros(L-1,N)];  map = patt;

aggr = [ones(1,N) ; ones(1,N) ; zeros(L-2,N)];

mol0 = [(L-1) + i*(1:N)'];  mols = [];

for j = 1:m
  mols = [mols ; mol0];
end

go = [];
while (isempty(go))

for j = 1:St

%  random step

  step = exp(i*floor(4*rand(size(mols)))*pi/2);
  mols = mols + step;

  for k = 1:size(mols,1)

%  periodic & reflection

    rw = real(mols(k));  iw = imag(mols(k));
    if (rw == L+1)
      mols(k) = (L-1) + i*iw;  rw = L-1;
    end
    
    if (iw == N+1)
      mols(k) = rw + i*1;  iw = 1;
    end
    
    if (iw ==   0)
      mols(k) = rw + i*N;  iw = N;
    end

%  aggregate & replenish

    if (aggr(rw,iw) == 1)
      if (rw == L-10)
        k = size(mols,1)-1;  go = 1
      end
      
      patt(rw,iw) = 1;
   	  aggr(rw+1,iw) = 1;  aggr(rw-1,iw) = 1;

	  if (iw == N)
%   	 aggr(rw,   1) = 1;
   	    aggr(rw,   1) = 1;  aggr(rw+1,   1) = 1;  aggr(rw-1,   1) = 1;
	  else
%	    aggr(rw,iw+1) = 1;
	    aggr(rw,iw+1) = 1;  aggr(rw+1,iw+1) = 1;  aggr(rw-1,iw+1) = 1;
	  end

	  if (iw == 1)
%   	  aggr(rw,   N) = 1;
   	    aggr(rw,   N) = 1;  aggr(rw+1,   N) = 1;  aggr(rw-1,   N) = 1;
	  else
%	    aggr(rw,iw-1) = 1;
	    aggr(rw,iw-1) = 1;  aggr(rw+1,iw-1) = 1;  aggr(rw-1,iw-1) = 1;
	  end

  	  mols(k) = (L-1) + i*floor(N*rand(1)+1);
    end
  end
end

%  plot

map = patt;
for k = 1:size(mols,1)

  rw = real(mols(k));  iw = imag(mols(k));
  map(rw,iw) = -1;
end

image(map,'cdatamapping','scaled')
axis equal;  axis image
pause(0.1)

end

%  final plot

image([map map],'cdatamapping','scaled')
axis equal;  axis image