
clear all
clf
niter = 2000;
global maxX maxY energy destroyed;
energy = 1;
prior_energy = 0;
maxX = 100;
maxY = 100;
destroyed = 0;
f = 0;
nNewParticlesPerFrame = 5;
map = [1.0, 1.0, 1.0
        0, 0.5, 0
        0, 1.0, 0
        1, 0, 0
        1, 0, 0
        0, 0, 1.0];

% Initialize matrix containing all 2D gird points A(x,y)
% 1 <= x <= maxX
% 1 <= y <= maxY
% A(x,y)=0 ... cite is empty
% A(x,y)>0 ... cite is filled
A = zeros(maxX, maxY);
bact = [];
mult = [];
xOlder = nan;
yOlder = nan;
light_src = 1.5;
yBuffer = floor(maxY*.25);
yStart = 1+yBuffer;
% Introduce a solid wall of phyto particles at y=1 
A(:,1) = .4;
for i = 1:37
    r = ceil(rand()*maxY);
    A(r,i) = .6;
    A(i,r) = 1.0;
    bact= [bact; r,i];
    mult = [mult; i,r];
end



for i = 1:niter;
	% Compute new starting point on line y=yStart
	x  = 1+floor(maxX*rand); 
	y  = 2+floor(yStart*rand()); 

    while 1
        xOld = x; 
        yOld = y;
        
        % 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
        if ((r*light_src)<.25)
            x = x+1;
        end
        
        if (.25<=(r*light_src)) && r<.5
            x = x-1;
        end
        if (.5<=r) && r<.75
            y = y+1;
        end
        if r>=.75 
            y = y-1;
        end
      
        if(x>maxX)
            x = 1;
        end
        if(x<1)
            x = maxX;
        end
        if (y>yStart)   % did the particle diffuse through the ceiling?
            y = yStart; % just set back to max if is gets too high
        end
        if(y==0)
            y = 1;
        end
        
        
        xR = x+1; % xRIGHT = x coordinate of my right neighbor
        xL = x-1; % xLEFT  = x coordinate of my left  neighbor
        yU = y+1; % yUP    = y coordinate of my upper neighor
        yD = y-1; % yDOWN  = y coordinate of my lower neighbor
        
        %%% Now apply periodic boundary conditions to xL and xR
        if(xR>maxX)
            xR = 1;
        end
        
        if(xL<1)
            xL = maxX;
        end
        if (yD == 0)
            yD = 1;
        end
       

        if (A(x,y) ~= 0)

            continue;
        end

        
      
       
       
            
   
        if (A(xR,y) == .4 || A(xL,y) == .4 || A(x,yU) == .4 || A(x,yD) == .4)
            A(x,y) = .4; % mark site as taken
            if (y+yBuffer>yStart && y+yBuffer<maxY) 
                yStart = y+yBuffer;
            end
            if(yStart > maxY)
                yStart = 1+yBuffer
            end
            fprintf('i= %5d x=%4d y=%4d yStart=%4d\n',i,x,y,yStart);
            break;
        end
    
    end
 
    [A,bact,mult] = determine_movement_alt(A,bact,mult);
          if(mod(i,nNewParticlesPerFrame)==0)
            if(prior_energy == energy)
                size1 = size(bact);
                size1 = size1(1); 
                if(size1 ~=0)
                 
                     A(bact(1,1),bact(1,2)) = 0;
                     bact = bact(2:size1,:);
                end
            end
            prior_energy = energy;
          end
   
    if (i==niter || mod(i,nNewParticlesPerFrame)==0 || yStart+1==maxY) 
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(map);
        axis equal;
        axis tight;
        if(light_src<1)
        title('Top Left Light Source');
        
        elseif(light_src>1)
        title('Top Right Light Source');
        else
        title('Central Light Source');
        end
        drawnow
%         saveas(gcf, ['Top_Right',sprintf('%03d',f),'.png']) 
%UNCOMMENT THIS TO SAVE IMAGES
        f =f+1;
    end

end

	