%flocking via agent-based
clear all
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);

N = 100; % number of flyers  %Remark: no flocking below N = 45
M = 100;  % number of time steps
L = 10; %Length of periodic domain
v = 1.0; %advective velocity
d = 0.01 % diffusion coefficient
dt = 0.05;  % time step
delt = sqrt(2*d*dt);
w = L/100;  %neighboring window size,
t = 0;

x = rand(1,N)*L; % initial position for flyers
 Xs = x; % save the positions
typ =  (1-2*(rand(1,N)>0.5));  % random distribution of initial directions
  tp = find(typ==1);
 tm = find(typ==-1);
 xp(1) = length(tp); % number of positive v flyers
 xm(1) = length(tm); % number of negative v flyers
 T(1) = t;
%make a plot
 figure(1)
 tp = find(typ==1);
 tm = find(typ==-1);
 plot(x(tp),typ(tp),'r*',x(tm),typ(tm),'b*')
 axis([0 L -1 1])
 pause(0.1)
 kstep = 20;
for k = 2:M
    for kk = 1:kstep
 t = t +dt;   
for j = 1:N
    xnb = find((abs(x-x(j))<w)|(abs(x+L-x(j))<w)|(abs(x-L-x(j))<w));  %find number of near neighbors
    W(j)=sum(typ(xnb))-typ(j);
end
y =  x+v*typ*dt +delt*randn(1,N);  %advect the points
yf = floor(y/L);
x = y-L*yf; % the domain is periodic;  map back to [0,L]
 Xs = [Xs;x]; % save the positions
 f = tanh(W)/2;
 al  = 1/2-typ.*f; %switching rate
 sw= (rand(1,N)<al*dt);
 typ = typ.*(1-2*sw); %switch directions sometimes
 
    end
    
 %make a plot
 figure(1)
 tp = find(typ==1);
 tm = find(typ==-1);
 xp(k) = length(tp);
 xm(k) = length(tm);
 T(k) = t;
 
 plot(x(tp),typ(tp),'r*',x(tm),typ(tm),'b*')
  axis([0 L -1 1])
  xlabel('Position','fontsize',20)
  ylabel('Direction of flight','fontsize',20)
 pause(0.1)
end


figure(2)
plot(T,xp/N,'r',T,xm/N,'b','linewidth',2)
axis([0 T(end) 0 1])
ylabel('proportion of flyers','fontsize',20)
xlabel('Time','fontsize',20)

 figure(3)
 %sample path
 th = 2*pi*Xs(:,1)/L;
 rad = L/(2*pi);
 tm = [1:size(Xs,1)]*dt;
 plot3(rad*cos(th),rad*sin(th),tm, 'linewidth',2)
xlabel('x','fontsize',20)
ylabel('y','fontsize',20)
zlabel('Time','fontsize',20)
title('Sample Path','fontsize',20)
 
