% code to simulate flocking behavior via Method of lines
function flock_sim
 global v N dx  D
 
 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7);



 % there are two populations Up and Um
 
 %parameters
 v = 0.5; % advective velocity
 L = 20; % length of domain  
 N = 100; %number of cells for method of lines
 D = 0.01 ; %diffusion coefficient
 
 
 dx = L/(N);
 x = [0:N-1]*dx;
 Ut = 1.65; %total population density; Total population is Ut*L
Ut = 1.95  % threshold is near Ut = 1
% setup for pde solve
% initial data
 Up = 0.52+0.1*sin(2*pi*x/L)/2;% slightly biased to favor positive over negative
 Um = 0.48-0.1*cos(4*pi*x/L)/2;
 s0 = Ut*[Up,Um]; 
tstep = 1;
t_end = 115;
 
%specify the output points
tspan = [0:tstep:t_end];
 
[T,S] = ode23s(@deRHS,tspan, s0);  

   figure(1)
% formatSpecF = '%6.2f\n';
 % title(strcat('Velocity = ',sprintf(formatSpecF,vel)),'fontsize',18)
  
 n = length(T);
 for j = 1:n
     sig = S(j,1:N)+S(j,N+1:2*N);
     W = S(j,1:N)-S(j,N+1:2*N);
     figure(1)
     plot(x,S(j,1:N),x,S(j,N+1:2*N),'linewidth',2)
     %axis([0 L 0 Ut])
     xlabel('x','fontsize',20)
     legend('U_+', 'U_-','fontsize',20)
     axis([0 L 0 3*Ut])
     pause(0.01)
     
 end


 figure(2)
 contour(S(1:n,1:N),[1,1],'linewidth',2)
 %  estimate the speed of the wave
 Tst = S(:,N/2);
 nsp = find((Tst(1:n-1)<=2)&(Tst(2:n)>2));
 figure(3)
 plot(T(nsp),nsp,'*')
 
 L/(T(nsp(end))-T(nsp(end-1)))
  
%the right hand side for ode simulation:
function s_prime=deRHS(t,y)
global v N dx D  
 Up = y(1:N);
 Um = y(N+1:2*N);
 vp = 0;
 W = Up-Um;
 f = tanh(W)/2;
 %switching rates
 al = 1/2-f;
 bt = 1/2+f;
  
%advection using periodic boundary conditions 
% use upwinding for the advective part
Fp = bt.*Um-al.*Up-v*(Up-[Up(end);Up(1:end-1)])/dx...
    +D*(-2*Up+[Up(end);Up(1:end-1)]+[Up(2:end);Up(1)])/dx^2;
Fm = al.*Up-bt.*Um+v*([Um(2:end);Um(1)]-Um)/dx...
    +D*(-2*Um+[Um(end);Um(1:end-1)]+[Um(2:end);Um(1)])/dx^2;

s_prime = [Fp;Fm];

