% BROWNIAN MOTION with invaders and a force field introduced.
% This is the movie I have all the jpeg files for, but had a hard time
% getting the Image program to work.

function invaders()
clear
clf 
%aviobj = avifile('space_invadersyellow','quality',100)

maxX=100;
maxY=100;
maxZ=100;
minX=0
minY=0
minZ=0
n=100; %number of particles
vt=2; %effective distance the particle travels in each step
%Start the particles inside of a "box"
boxX=50;
boxY=50;
boxZ=50;
steps=70;  % First few seconds is just Brownian motion. 
extra=120; % For these steps, there will be a force present


%insert initial particles in random positions. choose random positions in
%the plot. 
%i=updates particles, j=updates steps
for i=1:n
    
    px(i,1)=rand()*boxX ;% initial x positions
    py(i,1)=rand()*boxY ;% initial y positions
    pz(i,1)=rand()*boxZ; % initial Z positions
    theta(i,1)=rand()*2*pi;
%Need to update each particle separately, but at the same time!!
%each time the particle moves, we need a new theta.
end

for j=1:steps
          for i=1:n
              X(i,j)=0;
              Y(i,j)=0;
              Z(i,j)=0;
              
         theta(i,j+1)=rand()*2*pi ;%angle they will fly away at
         px(i,j+1)=px(i,j)+(vt*cos(theta(i,j+1))+X(i,j));
         py(i,j+1)=py(i,j)+(vt*sin(theta(i,j+1))+Y(i,j));
         pz(i,j+1)=pz(i,j)+(vt*cos(theta(i,j+1))+Z(i,j)); 
              for k=1:n
                if sqrt(px(k,j+1)^2+py(k,j+1)^2+ pz(k,j+1)^2)-sqrt(px(i,j+1)^2+py(i,j+1)^2+ pz(i,j+1)^2)<=2;
                theta(i,j)=pi+theta(i,j);

                end
          % store theta
          theta(i,j)=theta(i,j);
          
          
     
%apply periodic bounds
        if px(i,j+1)>maxX;
           px(i,j+1)=minX
         end
        if px(i,j+1)<minX;
           px(i,j+1)=maxX;
        end
 
        if py(i,j+1)>maxY ;  
           py(i,j+1) = minY; 
        end
        if py(i,j+1)<minY;
           py(i,j+1)=maxY;
        end
       
        if pz(i,j+1)>maxZ;   
           pz(i,j+1) = minZ; 
        end
        if pz(i,j+1)<minZ;
           pz(i,j+1)=maxZ;
        end
 end
          end

whitebg([0 0 0 ])
%h1 = figure(2)

%Z = peaks;
%for q = 1:1      

       % F(j) = getframe;
       % aviobj = addframe(aviobj,F(j));
        plot3( px(:,j+1),py(:,j+1),pz(:,j+1),'g.','Markersize',15)
        whitebg([0 0 0 ])
          axis manual
        axis([0 100 0 100 0 100])
        grid on
        
        %hold on %activating this command will trace the particle path
        pause(.2)
        
    %saveas(h1,strcat('brownian',num2str(999+j),'.jpg'));
    %fprintf('Frame %4d of %4d\n',q,1);
%end

end
%aviobj = close(aviobj);


for j=steps+1:extra
     for i=1:n
        
        
        %S force field
        X(i,j)=0;
        Y(i,j)=2;
        Z(i,j)=0;
         
        %update positions
         theta(i,j+1)=rand()*2*pi ;
         px(i,j+1)=px(i,j)+vt*cos(theta(i,j+1))+X(i,j);
         py(i,j+1)=py(i,j)+vt*sin(theta(i,j+1))+Y(i,j);
         pz(i,j+1)=pz(i,j)+vt*cos(theta(i,j+1))+Z(i,j); 
              
         
         %Collison still applies
         for k=1:n
                if sqrt(px(k,j+1)^2+py(k,j+1)^2+ pz(k,j+1)^2)-sqrt(px(i,j+1)^2+py(i,j+1)^2+ pz(i,j+1)^2)<=2;
                theta(i,j)=pi+theta(i,j);

                end
         
              
         % Boundary conditions still apply
         
        if px(i,j+1)>maxX;
           px(i,j+1)=minX
         end
        if px(i,j+1)<minX;
           px(i,j+1)=maxX;
        end
 
        if py(i,j+1)>maxY ;  
           py(i,j+1) = minY; 
        end
        if py(i,j+1)<minY;
           py(i,j+1)=maxY;
        end
       
        if pz(i,j+1)>maxZ ;  
           pz(i,j+1) = minZ; 
        end
        if pz(i,j+1)<minZ;
           pz(i,j+1)=maxZ;
        end
              
         end
         end
  
         axis([minX maxX minY maxY minZ maxZ])
            
%h1 = figure(2)

%Z = peaks;
%for q = 1:1 
%F(j) = getframe;
       % aviobj = addframe(aviobj,F(j));
        whitebg([0 0 0 ])
        plot3( px(:,j+1),py(:,j+1),pz(:,j+1),'r.','Markersize',15)
          axis manual
        axis([0 100 0 100 0 100])
        grid on
        
        %hold on %acjivajing jhis command will jrace jhe parjicle pajh
        pause(.2)
        
    %saveas(h1,sjrcaj('brownian2',num2sjr(999+j),'.jpg'));
    %fprinjf('Frame %4d of %4d\n',q,1);
end
 %aviobj = close(aviobj);
end    
  
    


       
   
