function M = Animator(L,N,Na,T,c,K) %Function to animate the 2D local attraction model with periodic %boundary conditions with Na asynchronously updating particles and %N-Na synchronous updating particles. rng('shuffle') %generate new random initial configuration % PARAMETERS % L =side length of square w periodic boundary conditions % N = Total number of particles % Na = number of particles using the Asynchronous heading/position update (so N-Na are using Synchronous update) % T = Simulation time % c= relative strength of attraction to the local center of mass % K = animation R=4; %interaction radius delta=0.5; %displacement per timestep %INITIATE PARTICLE POSITIONS AND HEADINGS P=zeros(N,3); % N particle matrix for i=1:N %For each particle P(i,1)=rand*L; %Initiate x-coordinate. P(i,2)=rand*L; %Initiate y-coordinate. P(i,3)=rand*2*pi-pi; %Initiate heading [-pi,pi]. end for k=1:T %Run simulation for T time steps %----Plotting------------------------------------------------ if mod(k,K)==0 plot(P(:,1),P(:,2),'k.','markersize',10); hold on for r=1:N plot([P(r,1),P(r,1)+cos(P(r,3))],[P(r,2),P(r,2)+sin(P(r,3))],'r-'); end hold off axis([0 L 0 L]); xlabel('X position') ylabel('Y position') axis manual U(k)=getframe; %makes a movie fram from the plot end %------------------------------------------------------------- % UPDATE THE Na ASYNCHRONOUSLY UPDATING PARTICLES for i=1:Na %FIND THE NEIGHBORS AND CALCULATE THE LOCAL ATTRACTION VECTOR (CM) %for particle i using the nhS.m function. %nhS.m calculates (for particle i) %1. the direction from it to the local center of mass of its neighbors (CM) %2. the number of non-self neighbors (n) CMAL=nhS(i,P,R,L); CM=CMAL(1,:); % Direction to the LCM n=CMAL(2,1); % Number of non-self neighbors %PARTICLE i's CURRENT HEADING (D) D=[cos(P(i,3)),sin(P(i,3))]; %PARTICLE i's NEW HEADING (Dir) if n==0 %If no neighbours no local interactions Dir=D; else %if at least one neighbour local interactions occur Dir=c*CM+D; end NormDir=(1/(Dir(1,1)^2+Dir(1,2)^2)^0.5)*Dir; %Normalized new heading P(i,3)=atan2(NormDir(1,2),NormDir(1,1)); %New heading angle %Update position P(i,1)=mod(P(i,1)+delta*NormDir(1,1),L); %New x-coordinate P(i,2)=mod(P(i,2)+delta*NormDir(1,2),L); %New y-coordinate end % UPDATE THE N-Na SYNCHRONOUSLY UPDATING PARTICLES PP=zeros(N,3); %Initiate temporary matrix to store updated positions and headings for Synchronously updating particles for i=Na+1:N % %FIND THE NEIGHBORS AND CALCULATE THE LOCAL ATTRACTION VECTOR (CM) %for particle i using the nhS.m function. %nhS.m calculates (for particle i) %1. the direction from it to the local center of mass of its neighbors (CM) %2. the number of non-self neighbors (n) CMAL=nhS(i,P,R,L); CM=CMAL(1,:); % Direction to the LCM n=CMAL(2,1); % Number of non-self neighbors %PARTICLE i's CURRENT HEADING (D) D=[cos(P(i,3)),sin(P(i,3))]; %PARTICLE i's NEW HEADING (Dir) if n==0 %If no neighbours no local interactions Dir=D; else %if at least one neighbour local interactions Dir=c*CM+D; end NormDir=(1/(Dir(1,1)^2+Dir(1,2)^2)^0.5)*Dir; %Normalized new heading PP(i,3)=atan2(NormDir(1,2),NormDir(1,1)); %Store new heading angle in temp matrix PP % Store new position in temp matrix PP PP(i,1)=mod(P(i,1)+delta*NormDir(1,1),L); %New x-coordinate PP(i,2)=mod(P(i,2)+delta*NormDir(1,2),L); %New y-coordinate end %UPDATE P MATRIX P(Na+1:N,:)=PP(Na+1:N,:); %Insert all Synchronously updated positions/angles (in PP) into the particle matrix P. %PERMUTE THE ROWS IN THE P MATRIX (To avoid particles being updated in a specific order in every time step) Q=P; dd=randperm(N); % d is a random permutation of the vector [1,2,3,..,N] for i=1:N % This loop randomly change the order of the rows, and hence the individual particles, in P. P(i,:)=Q(dd(i),:); end end function CMAL=nhS(i,P,R,L) %Neighborhood function that determines the neighbors of particle i taking into %account the periodic boundary conditions (9 cases) and calculates %the direction from particle i to the local center of mass of the neighbors. N=size(P,1); %Number of potential neighbours AG=zeros(N,3); %To store the coordinates and heading angles of the neighbors. count=0; %Initiate neighbor counter for j=1:N %For all particles if j~=i %not itself if (P(i,1)-P(j,1))^2+(P(i,2)-P(j,2))^20 %If have any neighbors CMx=mean(AG(:,1)); %x-coordinate of local center of mass CMy=mean(AG(:,2)); %y-coordinate of local center of mass CM=[CMx-P(i,1),CMy-P(i,2)]; %Direction toward the center of mass of the neighbors CM=(1/((CM(1,1)^2+CM(1,2)^2)^(1/2)))*CM; %Normalized direction toward the local center of mass CMAL=[CM;size(AG,1),0]; %Output else %If no neighbours CMAL=[0,0;0,0]; %Output end