clear;
clf;

f = 0; % keep the images numbered
nParticles = 1000; % CaCO3 particles introduced to fluid
maxX = 60;
maxY = 120;

nNewParticlesPerFrame = 10;

A = zeros(maxX, maxY);

A(10:11,40:80) = 1; % start with a small coral clipping in the flow


for i = 1:nParticles
    % compute particle release along the right boundary
    x  = maxX; % always at maxX
    y  = floor(rand()*maxY); % random entry to grid with respect to y
    
    while 1
        xOld = x; %store current x and y values
        yOld = y;
        
        % Arbitrarily select motion in flow path
        r = rand();
        
        if (r <= 0.75) % 75% flowmotion to the left
            x = x - 1;
        elseif (0.75 < r) && (r <= 0.875) % 12.5% up
            y = y + 1;
        elseif (0.875 < r) && (r <= 1) % 12.5 % down
            y = y - 1;
        end
        
        if (x < 10) % if particle reaches y-axis it flows past accreting precipitate
            break;
        end
        
        % to save computation time, if particle reaches top boundary or
        % x-axis, reflect to continue motion in flow path
        if (y > maxY)
            y = maxY;
        end
        if (y < 1)
            y = 1;
        end
        
        if (A(x,y)>0) % this site has been taken already, reject the move
            x=xOld;
            y=yOld;
            continue;
        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 neighbor
        yD = y-1; % yDOWN  = y coordinate of my lower neighbor
        
        if yU > maxY
            yU = 1;
        end
        if yD < 1;
            yD = maxY;
        end
        
        f = f + 1; % number the snapshots sequentially
        
        if  (A(xL,y) > 0 | A(x,yD) > 0 | A(x,yU) > 0)
            A(x,y) = 1;
            break;
        end
    end
    
    if (i==nParticles | mod(i,nNewParticlesPerFrame)==0 )
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(1-prism);
        axis equal;
        axis tight;
        drawnow
        saveas(gcf, ['DLA aaprecipitation',sprintf('%03d',f),'.png']);
    end
end


jParticles = 500; % acid ions in solution are less prolific than CaCO3 ions
                  % to slow dissolution
jNewParticlesPerFrame = 1;

for j = 1:jParticles
    % compute particle release along the right boundary
    x  = maxX; % always at maxX
    y  = floor(rand()*maxY); % random entry to grid with respect to y
    
    while 1
        xOld = x; %store current x and y values
        yOld = y;
        
        % Arbitrarily select motion in flow path
        r = rand();
        
        if (r <= 0.75) % 75% flowmotion to the left
            x = x - 1;
        elseif (0.75 < r) && (r <= 0.875) % 12.5% up
            y = y + 1;
        elseif (0.875 < r) && (r <= 1) % 12.5 % down
            y = y - 1;
        end
        
        if (x < 2) % if particle reaches x-axis it flows past accreting precipitate
            break;
        end
        
        % to save computation time, if particle reaches top boundary or
        % x-axis, reflect to continue motion in flow path
        if (y > maxY)
            y = maxY;
        end
        if (y < 1)
            y = 1;
        end
        
        if (A(x,y)>0) % this site has been taken already, reject the move
            x=xOld;
            y=yOld;
            continue;
        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 neighbor
        yD = y-1; % yDOWN  = y coordinate of my lower neighbor
        
        if yU > maxY
            yU = 1;
        end
        if yD < 1;
            yD = maxY;
        end
        f = f + 1; %number the snapshots sequentially
        if  ((A(x,yU) > 0) |  (A(x,yD) > 0) | (A(xL,y) > 0))
            A(x,yU) = 0;
            A(xL,y) = 0;
            A(x,yD) = 0;
            break;
        end
        
    end
    
    if (j==jParticles | mod(j,jNewParticlesPerFrame)==0 ) % |yStart+1==maxY)
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(1-prism);
        axis equal;
        axis tight;
        drawnow
        saveas(gcf, ['DLA ablation',sprintf('%03d',f),'.png']);
    end
end


