% L. Nielsen 
% Prof. Burkhard Militzer
% EPS 109 final movie project Fall 2009
% Ionic particle precipitation due to DLA through stagnant fluid

clear all
clf

nParticles = 7200;
maxX = 200;
maxY = 200;
maxZ = 200;

nNewParticlesPerFrame = 10;

% Initialize 3D matrix
A = zeros(maxX, maxY, maxZ);

% Introduce a sticky particle nucleus in the center of the grid
seed = 1;
A(99,100,99) = 1;
A(100,99,99) = 1;
A(99,99,99) = 2;
A(100,100,99) = 2;
A(99,100,100) = 2;
A(100,99,100) = 2;
A(99,99,100) = 1;
A(100,100,100) = 1;

% Create buffer region around particle from which to begin diffusion
Buffer = 10;
% To create a movie on a PC initialize: 
%aviobj = avifile('precip2.avi','quality',100,'compression','Cinepak');
shot = 0;

for i = 1:nParticles
    % Select type of particle (+ or -)
    if rand() < 0.5
        p = 1;
    else p = 2;
    end
    
    % Select starting point for diffusing ion
    checkx = 0;
    checky = 0;
    checkz = 0;
    while( (checkx ~= 1)||(checky ~= 1) || (checkz ~=1))
        xm = 1+floor(maxX*rand);
        ym = 1+floor(maxY*rand);
        zm = 1+floor(maxZ*rand);
        if (((xm <= ((maxX/2)+seed+Buffer+floor(i/50)))&&(xm > ((maxX/2)+seed+floor(i/50)))) || ((xm >= ((maxX/2)-seed-Buffer-floor(i/50))) && xm < ((maxX/2)-seed-floor(i/50))))
            x  = xm;
            checkx = 1;
        end
        if (((ym <= ((maxY/2)+seed+Buffer+floor(i/50)))&&(ym > ((maxY/2)+seed+floor(i/50)))) || ((ym >= ((maxY/2)-seed-Buffer-floor(i/50))) && ym < ((maxY/2)-seed-floor(i/50))))
            y  = ym;
            checky = 1;
        end
        if (((zm <= ((maxZ/2)+seed+Buffer+floor(i/50)))&&(zm > ((maxZ/2)+seed+floor(i/50)))) || ((zm >= ((maxZ/2)-seed-Buffer-floor(i/50))) && zm < ((maxZ/2)-seed-floor(i/50))))
            z  = zm;
            checkz = 1;
        end
    end
    
    while 1
        xOld = x; 
        yOld = y;
        zOld = z;

        % Select direction of movement
        r = rand();
        if r < .1667
            x = x-1; 
        elseif (r >= 0.1667 && r < 0.3333)
            x = x+1; 
        elseif (r >= 0.3333 && r < 0.5)
            y = y+1; 
        elseif (r >= 0.5 && r < 0.6667)
            y = y-1; 
        elseif (r >= 0.6667 && r < 0.8333)
            z = z-1; 
        elseif (r >= 0.8333)
            z = z+1; 
        end

        % Apply periodic boundary conditions in buffer box
        xUlim = (maxX/2)+seed+Buffer+floor(i/50);
        xLlim = (maxX/2)-seed-Buffer-floor(i/50);
        yUlim =(maxY/2)+seed+Buffer+floor(i/50);
        yLlim = (maxY/2)-seed-Buffer-floor(i/50);
        zUlim = (maxZ/2)+seed+Buffer+floor(i/50);
        zLlim = (maxZ/2)-seed-Buffer-floor(i/50);
        
        if (x > xUlim)
            x = xLlim; % if x exits to the right, let it enter on the left
        end
        if (x < xLlim)
            x = xUlim; % if x exits to the left, let it enter on the right
        end

        if (y > yUlim)   
            y = yLlim;
        end
        if (y < yLlim)   
            y = yUlim; 
        end
        
        if (z > zUlim)   
            z = zLlim; 
        end
        if (z < zLlim)   
            z = zUlim; 
        end
        
        if (A(x,y,z)>0) 
            x=xOld;
            y=yOld;
            z=zOld;
            continue;
        end

        % Apply periodic boundary conditions to movements
         xR = x+1;     % xRIGHT = x coordinate of my right neighbor
        if (xR > xUlim)
            xR = xLlim;
        end
         xL = x-1;     % xLEFT  = x coordinate of my left  neighbor
        if (xL < xLlim)
            xL = xUlim;
        end
         yU = y+1; 
        if (yU > yUlim)
            yU = yLlim;
        end 
         yD = y-1; 
        if (yD < yLlim)
            yD = yUlim;
        end
         zU = z+1; 
        if (zU > zUlim)
            zU = zLlim;
        end 
         zD = z-1; 
        if (zD < zLlim)
            zD = zUlim;
        end
         
        
        % If neighbor is occupied by oppositely-charged particle, form a
        % "bond"
        if p == 1           
            if  (A(xR,y,z) ==2 || A(xL,y,z) ==2 || A(x,yU,z) ==2 || A(x,yD,z) ==2 || A(x,y,zU) ==2 || A(x,y,zD) ==2)
                A(x,y,z) = p; % mark site as taken by positive ion
                fprintf('i= %5d x=%4d y=%4d z=%4d\n',i,x,y,z);
                break;
            end
        else if p == 2
                if  (A(xR,y,z) == 1 || A(xL,y,z) == 1 || A(x,yU,z) == 1 || A(x,yD,z) == 1 || A(x,y,zU) == 1 || A(x,y,zD) == 1)
                    A(x,y,z) = p; % mark site as taken by negative ion
                    fprintf('i= %5d x=%4d y=%4d z=%4d\n',i,x,y,z);
                    break;
                end
            end
        end
    end
     
    %alfa = 30; % For a stationary particle image
    %beta = 50; 
    
    % Plot the particle: 
    if (i==nParticles || mod(i,nNewParticlesPerFrame)==0)
        shot = shot + 1;  
        %To rotate the growing particle image
        if shot < 360
            alfa = shot;
        end
        if 360 <= shot < 720
            alfa = shot - 360;
        end
        beta = 50;
        % Intial values
        [l,m,n]=size(A);
        M=max(size(A));
        M1=round(sqrt(2)*M)+2;
        b=zeros(M1,M1);
        d=-M1*ones(M1,M1);
        alfa=alfa*pi/180;
        beta=beta*pi/180;
        T=zeros(3,3); % Rotation matrix
        T(1,1)=cos(alfa);
        T(1,2)=sin(alfa);
        T(1,3)=0;
        T(2,1)=-cos(beta)*sin(alfa);
        T(2,2)=cos(beta)*cos(alfa);
        T(2,3)=sin(beta);
        T(3,1)=sin(beta)*sin(alfa);
        T(3,2)=-sin(beta)*cos(alfa);
        T(3,3)=cos(beta);
        xyz=[1,1,1]';
        for p=1:l,
            for q=1:m,
                for r=1:n,
                    xyz=round(T*[p-l/2,q-m/2,r-n/2]');
                    x1=round(M1/2+xyz(1)+1);
                    y1=round(M1/2+xyz(2)+1);
                    z1=round(M1/2+xyz(3)+1);
                    if A(p,q,r)>0,
                        if d(x1,y1)<z1,
                            b(x1,y1)=A(p,q,r);
                            d(x1,y1)=z1;
                        end
                    end
                    
                    
                end
            end
        end
        imagesc(b)
        h = figure(shot);
        [fx,fy]=gradient(d);
        lm=(b.*(d.^1.5+mean(mean(d/2))*(cos((atan(fy)+atan(fx))/2)).^2));
        pcolor(flipud(lm));
        shading 'interp' ;
        saveas(h,strcat('jpgs\test_movie_',num2str(999+shot),'.jpg'));
        %aviobj = addframe(aviobj,F(shot));
        close all
    end
end
%aviobj = close(aviobj);