%...............................................................

%INITIAL CONDITIONS OF MATRIX/SEA
maxX = 200;
maxY = 150;
threshold = 5; %Number of nutrient particles needed for growth
A = zeros (maxX+4, maxY); %(A) is the region, with X,Y indices for location and threshold to count the amount of nutrients it has.
N = zeros (maxX+4, maxY); %Keeps track of the nutrient value at this particular (X,Y) in the matrix (A)
nParticles = 30000; %Number of nutrient particles that will diffuse into the system in total.
nNewParticlesPerFrame = 100;
nutrientpower = 2; %How much nutrient does each nutrientparticle have?

%Introduce seafloor conditions;create a slope at start %%%%%%%%%%%%
A(3:maxX+4,3) = 1;
for i = 3:90;
    A(i,floor(1/80*(4250-47*sqrt(165*i-i^2))))=1;
end
    
%Introduce buffer to speed up diffusion
yBuffer = 5 ;
yStart = 40+yBuffer;

%...............................................................

for i = 1:nParticles
    % Starting point of nutrient diffusion
    x  = 3+floor(maxX*rand);
    y  = yStart; %always start at upper limit

%...............................................................

%NUTRIENT DIFFUSION
    while 1
        xOld = x; %store current x and y values
        yOld = y; %main diffusion loop for nutrients movement
        r=rand();
        if r>0 && r<(1/4); %right
            x = x+1;
        else if r>(1/4) && r<=(2/4); %left
                x=x-1;
            else if r>(2/4) && r<=(3/4); %up
                    y=y+1;
                else if r>(3/4) && (r<=1); %down
                        y=y-1;
                    end
                end
            end
        end
        
        %NO PBC, restrict x to a margin of 2 away from matrix ends (2
        %spaces are needed to check. x will check x-1, which will check
        %x-1-1.
        if x<3 || x>maxX-2 ; %leave bufferzone for coral checking
           x = xOld; %trap on boundary
        end
                     
        %NO PBC FOR CEILING EITHER (herpderp)
        if (y>yStart)   
            y = yStart; % just set back to max if is gets too high
        end
        
        %IF CORAL IS ENCOUNTERED
        if (A(x,y)>0)
            x=xOld;
            y=yOld; %reject nutrient diffusion.
            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 neighor
        yD = y-1; % yDOWN  = y coordinate of my lower neighbor
        
                    
        
        %THEN NUTRIENT TO STOP DIFFUSING; CORAL CAPTURES NUTRIENT 
        if (A(xR,y) > 0 | A(xL,y) > 0 | A(x,yU) > 0 | A(x,yD) > 0)
            N(x,y) = N(x,y) + nutrientpower;
            if (y+yBuffer>yStart && y+yBuffer<maxY) 
                yStart = y+yBuffer;
            end
            break
        end        
    end
    
    
%...............................................................    
    
%CORAL GROWTH

    %Nutrient has stopped diffusing! Find the coral cause and determine how much
    %nutrients it has access to.
    if A(xR,y) > 0
        nutrientaccess = N(xR+1,y)+N(xR-1,y)+N(xR,y+1)+N(xR,y-1); %check all cells around it
        if nutrientaccess>calculatethreshold(A,xR,y);
            A(x,y) = 1;
            N(x,y) = 0; %This cell is no longer available for nutrients
        end
        
    else if A(xL,y) > 0
            nutrientaccess = N(xL+1,y)+N(xL-1,y)+N(xL,y+1)+N(xL,y-1); %check all cells around it
            if nutrientaccess>calculatethreshold(A,xL,y);
                A(x,y) = 1;
                N(x,y) = 0;
            end
            
        else if A(x,yU) > 0
                nutrientaccess = N(x+1,yU)+N(xR-1,yU)+N(xR,yU+1)+N(xR,yU-1); %check all cells around it
                if nutrientaccess>calculatethreshold(A,x,yU);
                    A(x,y) = 1;
                    N(x,y) = 0;
                end
                
            else if A(x,yD) > 0
                    nutrientaccess = N(xR+1,yD)+N(xR-1,yD)+N(xR,yD+1)+N(xR,yD-1); %check all cells around it
                    if nutrientaccess>calculatethreshold(A,x,yD);
                        A(x,y) = 1;
                        N(x,y) = 0;
                    end
                end
            end
        end
    end
    
%...............................................................

%DRAW IMAGE
    if (i==nParticles | mod(i,nNewParticlesPerFrame)==0 | yStart+1==maxY)
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(1-gray);
        axis equal;
        axis tight;
        drawnow
    end
    %saveas (gcf,['Coral',sprintf('%03d',f),'.png'])
    
%...............................................................
    
%STOP AT MAX HEIGHT
    if (yStart+1==maxY) 
        fprintf('Growth reached Y limit after only %d particles\n',i);
        break;
    end
end
        
