% SNOWING

% Original Author of DLA : B. Militzer, University of California, Berkeley 
% Date   : Sep 21, 2009

% Modified by: S. Liu, Dec 01 2013

clf;
clear all;

nParticles = 3000;
maxX = 150;
maxY = 50;

nNewParticlesPerFrame = 150;

% 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);
Z = zeros(maxX, maxY);
f=0;

% Introduce a solid wall of already aggreated particles at y=1 
A(:,1) = 1;

% To save computer time, we want to inject the new particle not too far
% above growing aggreate. We inject at on a line 'yStart', which
% keeps being increased so that it is always 'yBuffer' lines above the
% highest frozen particle
yBuffer = 5;
yStart = 1+yBuffer;

for i = 1:nParticles
	% Compute new starting point on line y=yStart
	x  = 1+floor(maxX*rand); 
	y  = yStart; %always start at upper limit

    while 1
        xOld = x; %store current x and y values
        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
        %%% 25% go RIGHT
        %%% 25% go UP
        %%% 25% go DOWN
        %%% --> need to enter about 10 lines of code
        
        if (r <= 0.25)
            x = x-1; %go left
        end
        
        if (r > 0.25 && r <= 0.5)
            x = x+1;
        end
        
        if (r>0.5 && r <= 0.75)
            y = y+1;
        end
        
        if (r>0.75)
            y=y-1;
        end
             
        %%% Now apply periodic boundary conditions in X direction
        %%% x must stay with 1 <= x <= maxX
        %%% if x exits to the left, let it enter on the right
        %%% if x exits to the right, let it enter on the left
        %%% --> enter 6 lines of code (or less)

        if (x<1)
            x = x+maxX;
        end
        
        if (x>maxX)
            x = 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 (A(x,y)>0) % ok, 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 neighor
        yD = y-1; % yDOWN  = y coordinate of my lower neighbo
        
        %%% Now apply periodic boundary conditions to xL and xR
 
        if (xR > maxX)
            xR = xR-maxX;
        end
        
        if (xL < 1)
            xL = xL+maxX;
        end
        
 		% Any of my 4 neighor filled already? If yes, get stuck here
        n=0;
        p = rand();
        p2 = 0.001;
        if (A(xR,y) > 0)
           n=n+1; 
        end
        
        if (A(xL,y) > 0)
            n=n+1;
        end
        
        if (A(x,yU) > 0)
            n=n+1;
        end
        
        if (A(x,yD) > 0)
            n=n+1;
        end
        
        if (n>=2)
            p2 = (n-1)*10*p2;
        end
        
        if ((A(xR,y) > 0 | A(xL,y) > 0 | A(x,yU) > 0 | A(x,yD) > 0) &&  p<=p2)
            A(x,y) = 1; % mark site as taken
            if (y+yBuffer>yStart && y+yBuffer<maxY) 
                yStart = y+yBuffer;
            end
            fprintf('i= %5d x=%4d y=%4d yStart=%4d\n',i,x,y,yStart);
            break;
        end

    end
    %{
    Z=-5.*A;
    xx= 0:maxX-1;
    yy=0:maxY-1;
    [X,Y]=meshgrid(yy,xx);
    %}
    if (i==nParticles | mod(i,nNewParticlesPerFrame)==0 | yStart+1==maxY)
        %subplot(1,2,1);
        imagesc(A');
        set(gca,'YDir','normal');
        colormap([0.05 0.05 0.55; 1 1 1]);
        axis equal;
        axis tight;
        title('Accumulation of Snow');
        
        %{
        subplot(1,2,2);
        hold off
        pcolor(Y,X,Z)
        shading interp
        hold on
        contour(X,Y,Z,30,'k');  
        %}
        
        drawnow
        saveas(gcf, ['snow_',sprintf('%03d',f),'.png'])
        f=f+1;
    end

    if (yStart+1==maxY) 
        fprintf('Growth reached Y limit after only %d particles\n',i);
        break;
    end
end


% heating
T=A.*(-5);
T(:,maxY)=20;
%map = [0 0.4 1;1 0 1; 0 0 0.5; 1 1 1];
nColors=5;
map=ones(nColors,3);
map(1,:)=[0 0.4 1];
map(nColors,:)=[1 1 1];
for i = 2:nColors-1
    for j = 1:3
        map(i,j)=rand()*map(nColors,j);
    end
end


subplot(2,1,1);

imagesc(A');
set(gca,'YDir','normal');
colormap(map);
axis equal;
axis tight;
title('The Melting of Snow (top) & Temperature Distribution (bottom)')
drawnow

subplot(2,1,2);


[X,Y]=meshgrid(1:maxY,1:maxX);

hold off
pcolor(Y,X,T)

shading interp
hold on
contour(Y,X,T,30,'k');
colorbar

saveas(gcf, ['snow_',sprintf('%03d',f),'.png'])
f=f+1;
saveas(gcf, ['snow_',sprintf('%03d',f),'.png'])
f=f+1;

iter = 1500;
for k=[1:iter]
    Tnew = T;
    for i=2:maxX-1
        for j = 2:maxY-1
            Tnew(i,j)=(T(i,j+1)+T(i,j-1)+T(i-1,j)+T(i+1,j))/4;
            %fprintf('Tnew(%d,%d) = %d\n', i,j,Tnew(i,j));
            if (A(i,j)==1 && Tnew(i,j)>=0)
               A(i,j)=0; 
            end
        end
    end
    for i=2:maxX-1
        for j = 1
            Tnew(i,j)=(T(i,j+1)+T(i-1,j)+T(i+1,j))/3;
            %fprintf('Tnew(%d,%d) = %d\n', i,j,Tnew(i,j));
            if (A(i,j)==1 && Tnew(i,j)>=0)
               A(i,j)=0; 
            end
        end
    end
    for i = 1
        for j=2:maxY-1
            Tnew(i,j)=(T(i,j+1)+T(i,j-1)+T(i+1,j))/3;
            %fprintf('Tnew(%d,%d) = %d\n', i,j,Tnew(i,j));
            if (A(i,j)==1 && Tnew(i,j)>=0)
               A(i,j)=0; 
            end
        end
    end
    for i = maxX
        for j=2:maxY-1
            Tnew(i,j)=(T(i,j+1)+T(i,j-1)+T(i-1,j))/3;
            %fprintf('Tnew(%d,%d) = %d\n', i,j,Tnew(i,j));
            if (A(i,j)==1 && Tnew(i,j)>=0)
               A(i,j)=0; 
            end
        end
    end
    Tnew(1,1)=(T(1,2)+T(2,1))/2;
    if (A(1,1)==1 && Tnew(1,1)>=0)
       A(1,1)=0; 
    end
    Tnew(maxX,1)=(T(maxX-1,1)+T(maxX,2))/2;
    if(A(maxX,1)==1&&Tnew(maxX,1)>=0)
       A(maxX,1)=0; 
    end
    
    
    fprintf('iteration = %d\n', k);
    T=Tnew;
    
    if (mod(k,30)==1)

        subplot(2,1,1);
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(map);
        axis equal;
        axis tight;
        title('The Melting of Snow (top) & Temperature Distribution (bottom)')
        drawnow
        
        subplot(2,1,2);
        
        hold off
        pcolor(Y,X,T)
        
        shading interp
        
        hold on
        contour(Y,X,T,30,'k');
        colorbar
        
        saveas(gcf, ['snow_',sprintf('%03d',f),'.png'])
        f=f+1;
    end
end
%{
imagesc(A');
set(gca,'YDir','normal');
colormap(map);
axis equal;
axis tight;
drawnow
%}