clear
clf

% Set up grid size
N = 100;
L = 20;
maxX = N + 1;
maxY = N + 1;

% Set up initial DLA
nParticles = 2500;
nIterPerImage = 100;
A = zeros(N+1, N+1);
A(:,1) = 1;

% Set up initial T
sourceT = 50;
T = zeros(maxX, maxY);
for i = 1:maxX
    for j = 1:maxY
        if ((i-maxX)^2 + (j-maxY)^2 < N/3)
            T(i, j) = sourceT;
        end
    end
end 
T_new = T;

% For plotting
xx = 0:L/N:L;
yy = 0:L/N:L;
[X,Y] = meshgrid(xx,yy);

% Set up DLA stuff
yBuffer = 5;
yStart = 1+yBuffer;
framesInUpdate = 20;

% Set up the colormap
cmap = colormap;
lastMap = size(cmap,1)
cmap(lastMap,:) = [0 1 0];
colormap(cmap);

img = 0;

for iter = 1:nParticles
    disp(iter);
    
    x  = 1+floor(maxX*rand); 
	y  = yStart; %always start at upper limit
    
    while 1 
        xOld = x;
        yOld = y;
        
        % Random walk
        r = rand();
        if (r < .25)
            x = x - 1;
        elseif (r < .5)
            x = x + 1;
        elseif (r < .75)
            y = y + 1;     
        elseif (r < 1)
            y = y - 1;    
        end

        % Periodic boundary conditions
        if (x < 1)
            x = maxX;
        end
        if (x > maxX)
            x = 1;
        end

        if (y > yStart)
            y = yStart;
        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 < 1)
            yD = maxY;
        end
        
        if (A(x,y) > 0)
            x = xOld;
            y = yOld;
            continue;
        end
        
        if (A(xR,y) > 0 | A(xL,y) > 0 | A(x,yU) > 0 | A(x,yD) > 0)
            tRand = rand();
            p = tRand * exp(T(x, y)) / 2;
            % Grow if not too hot
            if (p < 1)
                A(x,y) = 1;
            end
            if (y+yBuffer>yStart && y+yBuffer<maxY) 
                yStart = y+yBuffer;
            end
            break;
        end
    end
    
    if (yStart+1==maxY) 
        fprintf('Growth reached Y limit after %d particles\n',i);
        break;
    end
    
    % Update temperature
    tempUpdateRate = 4;
    for tempIter = 1:tempUpdateRate
    for i = 1:maxX
        for j = 1:maxY
            if ((i-maxX)^2 + (j-maxY)^2 < N/3)
                T(i, j) = sourceT;
            end
        end
    end 
    
    for ii = 1:maxX
        T(ii,1) = T(ii,2);
        T(ii,maxY) = T(ii,maxY-1);
    end
    
    for jj = 1:maxY
        T(1,jj) = T(2,jj);
        T(maxX,jj) = T(maxX-1,jj);
    end
    
    for i = 2:N
        for j = 2:N
            T_new(i,j) = .25 * (T(i+1,j) + T(i-1,j) + T(i,j+1) + T(i,j-1));
        end
    end
    
    for ii = 1:maxX
        T_new(ii,1) = T_new(ii,2);
        T_new(ii,maxY) = T_new(ii,maxY-1);
    end
    
    for jj = 1:maxY
        T_new(1,jj) = T_new(2,jj);
        T_new(maxX,jj) = T_new(maxX-1,jj);
    end
    
    T = T_new;
    end

    if (mod(iter, framesInUpdate) == 0)
        Z = T;
        for i = 1:maxX
            for j = 1:maxY
                if (A(i, j) == 1)
                    Z(i, j) = sourceT + 1;
                end
            end
        end
        imagesc(Z');
        set(gca,'YDir','normal');
        hold off
        if (mod(iter, nIterPerImage) == 0)
            img = img + 1;
            saveas(gcf, ['finalproj',sprintf('%03d',img),'.png'])
        end
        pause(0.5);
    end
end
