% David Barth
% EPS 109 Final Project
% Simulation of porous silicon formation by electrochemical etching

% Modified from code written by B. Militzer, University of California, Berkeley
%% Set Initial Conditions
clear; clf;
rng(1);
nParticles = 50000;
maxX = 100;
maxY = 100;
maxZ = 100;
q=0.005; %Probability of Dissolution of SiO2
f=1;
nNewParticlesPerFrame = 1000;

A = zeros(maxX, maxY, maxZ);

% Seed domain
A(:,:,1) = 2;

%% Random Walk
zBuffer = 5;
zStart = 1+zBuffer;

for i = 1:nParticles
    %i;
    % Add new particle in silicon region
    x  = 1+floor(maxX*rand);
    y=1+floor(maxY*rand);
    z  = zStart; %always start at upper limit
    
    
    % Add new particle in
    while 1
        xOld = x; %store current x and y values
        yOld = y;
        zOld = z;
        
        % Arbitrarily select direction to turn to
        r = rand();
        
        if (r<1/6)
            x=x+1;
        elseif (r<2/6)
            x=x-1;
        elseif (r<3/6)
            y=y+1;
        elseif (r<4/6)
            y=y-1;
        elseif (r<5/6)
            z=z+1;
        elseif (r<6/6)
            z=z-1;
        end
        
        
        %%% Periodic Boundary Conditions
        
        
        if (x<=0)
            x=maxX;
        elseif (x>maxX)
            x=1;
        end
        if (y<=0)
            y=maxY;
        elseif (y>maxY)
            y=1;
        end
        
        
        if (z>zStart)   % did the particle diffuse through the ceiling?
            z = zStart; % just set back to max if is gets too high
        end
        
        if (A(x,y,z)>0) % ok, this site has been taken already, reject the move
            x=xOld;
            y=yOld;
            z=zOld;
            continue;
        end
        
        xR = x+1; % xRIGHT = x coordinate of my right neighbor
        xL = x-1; % xLEFT  = x coordinate of my left  neighbor
        yR = y+1; % yUP    = y coordinate of my upper neighor
        yL = y-1; % yDOWN  = y coordinate of my lower neighbo
        zU = z+1;
        zL = z-1;
        %%% Now apply periodic boundary conditions to xL and xR
        
        if (xL==0)
            xL=maxX;
        elseif (xR>maxX)
            xR=1;
        end
        if (yL==0)
            yL=maxY;
        elseif (yR>maxY)
            yR=1;
        end
        if (zU>maxZ)
            zU=z;
        end
        
        % Any of my 4 neighor filled already? If yes, get stuck here
        neighbor = [A(xR,y,z),A(xL,y,z),A(x,yL,z),A(x,yR,z),A(x,y,zU),A(x,y,zL)];
        if (sum(neighbor(neighbor==2))>0)
            A(x,y,z) = 1; % mark site as oxidized
            if (z+zBuffer>zStart && z+zBuffer<maxZ)
                zStart = z+zBuffer;
            end
            break;
        end
        
    end
    
    % Plot 3D figure. Markers show etched area.
    if (i==nParticles || mod(i,nNewParticlesPerFrame)==0 || zStart+1==maxZ)
        i
        clf;
        C=A;
        C(:,:,1)=0;
        C(C==1)=0; %Do not plot oxidized areas
        C(C==2)=1;
        [xx,yy,zz]=ind2sub(size(C), find(C));
        plot3(xx,yy,zz,'o','markeredgecolor',[0.5,0.5,0.5],'markersize',3)
        axis equal;
        axis tight;
        drawnow
        axis([0 maxX 0 maxY 0 maxZ])
        %saveas(gcf, ['project3d',sprintf('%03d',f),'.png']);
        f=f+1;
    end
    B=A;
    B(B==2)=0; %Find all oxidized areas
    ind=find(B);
    for kk=1:length(ind)  %At every step, each oxidized area might dissolve
        r2=rand;
        if (r2<q)
            A(ind(kk))=2; 
        end
    end
    if (zStart+1==maxZ)
        fprintf('Growth reached Z limit after only %d particles\n',i);
        break;
    end
end

%% Plot Cross Sections
C=A;
C(:,:,1)=0;
C(C==1)=0;
C(C==2)=1;
cssize=4;  %Plot region 4 units thick
for k = 1:length(C)-cssize
    D=C(:,k:k+cssize,:);
    [xcs,ycs,zcs]=ind2sub(size(D), find(D));
    plot3(xcs,ycs,zcs,'s','markeredgecolor',[0.5 0.5 0.5],'markerfacecolor',[0.5 0.5 0.5],'markersize',7);
    view(0,0)
    axis([0 maxX 0 maxY 0 maxZ])
    pause(0.3)
    %saveas(gcf, ['projectxs',sprintf('%03d',k),'.png']);
end