clear all
clf
close all

nParticles = 20000;
nNewParticlesPerFrame=100;
map =[1, 1, 1
    .6, .1, 0
    .9, 0, 0];
maxX = 800;
height=400;
h = zeros(maxX, height); %array that has 2 indices
height1=zeros(maxX,1);
h(:,1)=1; %base for volcano
height1(:)=1; %everywhere there is one particle
y=2;
f=0;

c=1;
for i = 1:nParticles
    x=maxX/2-1;
    x1=maxX/2;
    r=rand();
    w=1; %increase this, decrease the volcano's width
    while 1
        dh=height1(x)-height1(x-w); %this is the elevation fall, increase height(x-___), decrease volcano's slope
        dh1=height1(x1)-height1(x1+w);
        
        if dh>1  || (dh==1 && r>.1 ) %larger than 1, go left, or 2 chance of 50%
            x=x-1;
        else
            break
        end
        
        if dh>1 || (dh==1 && r>.1 )
            x1=x1+1;
        else
            break
        end
    end
    height1(x)=height1(x)+1; %how many particles are in this column
    height1(x1)=height1(x1)+1;
    h(x,height1(x))=c; %sticks particles, where the particle gets stuck
    %want to get particles colored
    h(x1,height1(x1))=c;
    if mod(i,1000)==0
        c=3-c;
    end
    if height1(x)>height %stay in bounds
        break
    elseif height1(x1)>height
        break
    end
    
    
    
    if (i==nParticles || mod(i,nNewParticlesPerFrame)==0 )
        imagesc(h');
        set(gca,'YDir','normal');
        colormap(map);
        axis equal;
        axis tight;
        drawnow
        %         f=f+1;
        %         saveas(gcf, ['VOLCANO_', sprintf('%03d',f), '.png'])
        %
        %
        
    end
title('w=1')   
  
end


