function flock

global alpha beta radat radrep delta_rep

format long

% system parameters
radius=1;
L=20;
noise=0.;


% numerical parameters
% initialize the random number generator in a well-defined state (for
% repeating the exactly same initial condition with different dt
state=137;
random_number_generator='state';
rand('state',state);
dt=0.2;
tmax=200;
numsteps=tmax/dt+1;
method=1;

% plotting options
log_plot_movie=1;

% initial conditions
Nboids=50;
x0=L*rand(Nboids,1);
y0=L*rand(Nboids,1);
vx0=2*(rand(Nboids,1)-.5);
vy0=2*(rand(Nboids,1)-.5);


% preallocate the memory for the positions to speed up the code
x=zeros(Nboids,numsteps);
y=zeros(Nboids,numsteps);
vx=zeros(Nboids,numsteps);
vy=zeros(Nboids,numsteps);
angle=zeros(Nboids,numsteps);
vxaverage=zeros(numsteps,1);
vyaverage=zeros(numsteps,1);

% start integration
x(:,1)=x0;
y(:,1)=y0;
vx(:,1)=vx0;
vy(:,1)=vy0;
xp=x0;
yp=y0;
time=0;

if (log_plot_movie==1) % set-up for movie
    HandleFig=figure(1);
    subplot(1,4,[1:3])
    HandleM=plot(xp,yp,'o','MarkerSize',5);
    set(gcf,'DoubleBuffer','on');
    axis([0 L 0 L]);
end

% integration loop 
dti=dt;
k=1;
while (time < tmax)
    if (time+dt>tmax)
        dti=tmax-time;
    end
    if (method==1) % forward Euler
        [Fx,Fy,Fvx,Fvy] = F(x(:,k),y(:,k),vx(:,k),vy(:,k),radius,Nboids,L);
        x(:,k+1) = x(:,k) + dti*Fx;
        y(:,k+1) = y(:,k) + dti*Fy;
        vx(:,k+1) = vx(:,k) + dti*Fvx;
        vy(:,k+1) = vy(:,k) + dti*Fvy;
    end
    
% determine order parameter    
    vxaverage(k+1,1)=mean(vx(:,k+1));
    vyaverage(k+1,1)=mean(vy(:,k+1));
% add a bit of noise to velocities
    if (noise>0)
        vx(:,k+1)=vx(:,k+1)+noise*2*(rand(Nboids,1)-.5);    
        vy(:,k+1)=vy(:,k+1)+noise*2*(rand(Nboids,1)-.5);
        vnorm=sqrt(vx(:,k+1).^2+vy(:,k+1).^2);
        vx(:,k+1)=vx(:,k+1)./vnorm;    
        vy(:,k+1)=vy(:,k+1)./vnorm;
    end
    k=k+1;
% periodic boundary conditions
    xp=mod(x(:,k),L);
    yp=mod(y(:,k),L);
    angle(:,k)=atan2(vy(:,k),vx(:,k));
    time=time+dti;
    if (log_plot_movie==1) % update plot    
        set(HandleM,'XData',xp,'MarkerFaceColor','r');
        set(HandleM,'YData',yp);
        subplot(1,4,4)
        imagesc(angle(:,k))
        caxis([-pi pi])
        colorbar
        drawnow
    end
end
kmax=k;

% output trajectories

figure(2);
plot(x',y')
title(' Trajectories')
xlabel(' x-position')
ylabel(' y-position')

figure(3)
plot(angle','o')

xlabel(' Time Steps')
ylabel(' Angle')

figure(4)
subplot(1,2,1)
vaverage=sqrt(vxaverage.^2+vyaverage.^2);
plot(vaverage)
xlabel(' Time Steps')
ylabel(' Order Parameter')
subplot(1,2,2)
angleaverage=atan2(vyaverage,vxaverage);
plot(angleaverage)
xlabel(' Time Steps')
ylabel(' Average Angle')

% output positon for convergence test
mod(x(1,numsteps),L)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [Fx,Fy,Fvx,Fvy] = F(x,y,vx,vy,radius,N,L)

global alpha beta radat radrep delta_rep

Fx=vx;
Fy=vy;
Fvx=zeros(N,1);
Fvy=zeros(N,1);

xp=mod(x,L);
yp=mod(y,L);

for i=1:N
    distancex=xp-xp(i)*ones(N,1);
    distancey=yp-yp(i)*ones(N,1);
    distancex(find(distancex>L/2))=distancex(find(distancex>L/2))-L;    
    distancex(find(distancex<-L/2))=distancex(find(distancex<-L/2))+L;
    distancey(find(distancey>L/2))=distancey(find(distancey>L/2))-L;    
    distancey(find(distancey<-L/2))=distancey(find(distancey<-L/2))+L;
    % periodic boundary conditions
    distancesq=max(distancex.^2+distancey.^2,.1e-10);
    weights=exp(-distancesq/radius^2);
    vxaverage=vx'*weights/N;
    vyaverage=vy'*weights/N;
    xaverage=x'*weights/N;
    yaverage=y'*weights/N;
    xaverage=mean(x);
    yaverage=mean(y);
    vavnorm=sqrt(vxaverage^2+vyaverage^2);
    Fvx(i,1)=vxaverage/vavnorm-vx(i);
    Fvy(i,1)=vyaverage/vavnorm-vy(i);
end

end % F.m
