% Douglas Fir generation
%based on Diffusion Limited Aggregation by Burkhard Militzer

clear;
clf;

nParticles = 12000;
%make square for easier branch design
maxX = 300;
maxY = 300;

nNewParticlesPerFrame = 100;

% empty "air"
A = zeros(maxX, maxY);

% sticky ground
A(:,1) = 1;
A(150,1:180)=1;
x0 = 1;
h = 240; %tree height and trunk
for y = 1:40:h %branches: 4 of decreasing size
    A(x0+60:h-x0,y) = 1;
    x0=x0+20;
end


	

yBuffer = 2; %smaller buffer size for more growth in lower branches
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;
        
        % random walk
		r = rand();

        if r < .25;
            x=x-1;
            elseif r < .50;
            x=x+1;
            elseif r < .75
            y=y+1;
            elseif r < 1.0;
            y=y-1;
        end
       
        %xboundary: go through wall

        if x == maxX+1
            x = 1;
        end
        if x == 0
            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 gets too high, reset at bottom and begin again
        if yStart>=200  
            yStart = 20;
            y = yStart;
        end
        
        if yStart>maxY
            yStart =  20;
            y = yStart;
        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 cord of my right friend
        xL = x-1; % xLEFT  = x cord of left pal
        yU = y+1; % yUP    = y cord of my upper chum
        yD = y-1; % yDOWN y cord of lower buddy
        
        %%% x neighbor boundary = allow particles to flow through walls if
        %%% neighbor is wall
        if xR == maxX+1
            xR = 1;
        end
        if xL == 0
            xL = maxX;
        end
            rr=rand(1); %number to beat
            s=.1; %stick chance
        
 		 
        %neighbor counter
        n=0;
        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
        
        % I want a low chance of initial stick, but a better chance to
        % stick once surrounded by 3 or 4 neighbors, better simulating
        % branch and leave growth
        if  n==1 && rr<s*.1 || n==2 && rr<.05*s || n==3 && rr<4*s || n==4 && rr<2.5*s
                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

    if (i==nParticles || mod(i,nNewParticlesPerFrame)==0 || yStart+1==maxY) 
        imagesc(A');
        set(gca,'YDir','normal');
        colormap(summer);
        axis equal;
        axis tight;
        drawnow
    end
%MOVIE CODE!!!! every 100th frame save for later compression
     if mod(i,100) ==1
        h1 = A(2);
        saveas(h1,strcat('dougfir',num2str(999+i),'.jpg'));
        fprintf('Frame %4d of %4d\n',i,nParticles);
     end

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

