clear;

nParticles = 300;
maxD = 30;
global maxD;

nNewParticlesPerFrame = 10;

% Initialize matrix containing all 2D gird points A(x,y)
% 1 <= x <= maxDX
% 1 <= y <= maxDY
% A(x,y)=0 ... cite is empty
% A(x,y)>0 ... cite is filled
A = zeros(maxD, maxD, maxD);


%A(maxD/2,maxD/2,maxD/2) = 1;
%A(:,:,1) = 1;
%A(:,:,maxD) = 1;
%A(:,1,:) = 1;
%A(:, maxD, :) = 1;
%A(1, :, :) = 1;
%A(maxD, :, :) = 1;
  center = maxD/2;
  radius = 14;

  for x=1:maxD
      for y=1:maxD
          for z=1:maxD
              dist = norm([center,center,center]-[x,y,z]);
                 if (dist >= radius && dist < radius+2)
                     A(x,y,z) = 1;
                 end
          end
      end
  end

for i = 1:nParticles

    x = 1+ floor(maxD*rand);
    y = 1+ floor(maxD*rand);
    z = 1+ floor(maxD*rand);
        
    x = center;
    y = center;
    z = center;

    while 1
        xOld = x; %store current x and y values
        yOld = y;
        zOld = z;
        
        % Arbitrarily select direction to turn to
		r = rand();
        %%% Determine the direction of motion using the random number 'r'
        %%% Update 'x' and 'y' with the following probabilities
        %%% 25% go LEFT
        %%% 25% go RIGHT
        %%% 25% go UP
        %%% 25% go DOWN
        %%% --> need to enter about 10 lines of code
        if r < 1/6
            x = xOld-1;
        elseif r < 2/6
            x = xOld+1;
        elseif r < 3/6
            y = yOld+1;
        elseif r < 4/6;
            y = yOld-1;
        elseif r < 5/6;
            z = zOld +1;
        else
            z = zOld -1;
        end
       
        x = boundaries(x);
        y = boundaries(y);
        z = boundaries(z);
                   
          
       if (A(x,y,z)>0) 
          break;
       end

        
        xR = boundaries(x+1);
        xL = boundaries(x-1);
        yR = boundaries(y+1);
        yL = boundaries(y-1);
        zU = boundaries(z+1);
        zD = boundaries(z-1);
        
        
      
        
 		% Any of my 4 neighor filled already? If yes, get stuck here
        if (A(xR,y,z) > 0 | A(xL,y,z) > 0 | A(x,yR,z) > 0 | A(x,yL,z) > 0 | A(x,y,zU)>0 | A(x,y,zD) > 0)
                fprintf('x=%d,y=%d,z=%d\n', x,y,z);
            A(x,y,z) = 1; % mark site as taken
            break;
        end

    end
    
    
    % create B, a matrix like A except with borders removed
    B = zeros(maxD, maxD, maxD);

      for x=1:maxD
        for y=1:maxD
          for z=1:maxD;
              dist = norm([center,center,center]-[x,y,z]);
                 if (dist > radius && dist < radius+2)
                     B(x,y,z) = 0;
                 else
                     B(x,y,z) = A(x,y,z);
                 end
          end
         end
      end

    if (i==nParticles | mod(i,nNewParticlesPerFrame)==0) 
       plotBetter(B);
      % fprintf('plotted');
       pause(0.001);
    end
end

	