%Generating an Eden cluster in three spatial dimensions
N=50;
lattice=ones(N,N,N);
s=N/2;
lattice(s,s,s)=2;
list(:,1)=[s s-1 s];
list(:,2)=[s s+1 s];
list(:,3)=[s-1 s s];
list(:,4)=[s+1 s s];
list(:,5)=[s s s+1];
list(:,6)=[s s s-1];
Nperm=6;


tag=0;

  while tag==0
    
   for i=1:N
       for j=1:N
            if (lattice(i,j,1)==2) 
                tag=1;
            end
            if (lattice(i,j,N)==2) 
                tag=1;
            end
            if (lattice(1,i,j)==2) 
                tag=1;
            end
            if (lattice(N,i,j)==2) 
                tag=1;
            end
            if (lattice(i,1,j)==2) 
                tag=1;
            end
            if (lattice(i,N,j)==2) 
                tag=1;
            end
       end
   end
    if tag==1 
        break
    end
    
    r=ceil(rand()*Nperm); %pick perimeter points at random
    x=list(1,r); %get x,y for perimeter point picked
    y=list(2,r);
    z=list(3,r);
    lattice(x,y,z)=2;
    %take point out of perimeter list
    list=[list(:,1:r-1) list(:,r+1:end)];
    Nperm=Nperm-1;
    
    %add new perimeter points to perimeter list if aren't in cluster
    %already
    xnew=x-1;
    ynew=y;
    znew=z;
    if lattice(xnew,ynew,znew)==1
        lattice(xnew,ynew,znew)=2;
        list=[list [xnew ynew znew]'];
        Nperm=Nperm+1;
    end
    xnew=x+1;
    ynew=y;
    znew=z;
    if lattice(xnew,ynew,znew)==1
        lattice(xnew,ynew,znew)=2;
        list=[list [xnew ynew znew]'];
        Nperm=Nperm+1;
    end
    xnew=x;
    ynew=y-1;
    znew=z;
    if lattice(xnew,ynew,znew)==1
        lattice(xnew,ynew,znew)=2;
        list=[list [xnew ynew znew]'];
        Nperm=Nperm+1;
    end
    xnew=x;
    ynew=y+1;
    znew=z;
    if lattice(xnew,ynew,znew)==1
        lattice(xnew,ynew,znew)=2;
        list=[list [xnew ynew znew]'];
        Nperm=Nperm+1;
    end
    xnew=x;
    ynew=y;
    znew=z+1;
    if lattice(xnew,ynew,znew)==1
        lattice(xnew,ynew,znew)=2;
        list=[list [xnew ynew znew]'];
        Nperm=Nperm+1;
    end
    xnew=x;
    ynew=y;
    znew=z-1;
    if lattice(xnew,ynew,znew)==1
        lattice(xnew,ynew,znew)=2;
        list=[list [xnew ynew znew]'];
        Nperm=Nperm+1;
    end
end

%plotting the cluster in 3d
l=1;
for i=1:N
    for j=1:N
        for k=1:N
            if (lattice(i,j,k)==2) 
                xx(l)=i;
                yy(l)=j;
                zz(l)=k;
                l=l+1;
            end
        end
    end
end
plot3(xx,yy,zz,'b+')

  % calculate the number of cluster sites (mass) within a given radius to be added below